Modelling
and Simulation
of Turbulence
and Heat Transfer
in Wall-Bounded Flows
kada sam imao petnaest godina naˇsao sam se u uˇzasu rata uslijed agresije na moju domovinu od strane njenih susjeda. proˇsavˇsi to teˇsko iskustvo ja sam nastavio moj put prema mom velikom ˇzivotnom cilju, a ˇciji rezultat se upravo nalazi pred vama. naˇzalost, mnogi moji sunarodnjaci nikad nisu dobili priliku da ostvare svoj ˇzivotni cilj jer su brutalnom silom zaustaviljeni. zbog toga ovaj moj rad posve´cujem svim neduˇznim ljudima koje je na njihovom ˇzivotnom putu zaustavila agresija na moju domovinu, bosnu i hercegovinu....
Modelling and Simulation of Turbulence and
Heat Transfer in Wall-Bounded Flows
PROEFSCHRIFT
ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,
op gezag van de Rector Magnificus Prof. dr. ir. J.T. Fokkema, voorzitter van het College voor Promoties,
in het openbaar te verdedigen op maandag 16 october 2006 te 10:00 uur
door
Mirza POPOVAC
Samestelling promotiecommissie:
Rector Magnificus voorzitter
Prof. dr. dipl.ing. K. Hanjali´c Technische Universiteit Delft, promotor
Prof. dr. D. Laurence University of Manchester, Groot Brittanni¨e; EDF, Frankrijk
Prof. dr. D.J.E.M Roekaerts Technische Universiteit Delft
Prof. dr. ir. Th.H. v.d. Meer Universiteit Twente
Prof. dr. ir. C.R. Kleijn Technische Universiteit Delft
dr. ir. B.J. Boersma Technische Universiteit Delft
dr. ir. M. Tummers Technische Universiteit Delft
Copyright by Mirza Popovac c° 2006.
All rights reserved. ISBN-10: 90-9021117-9 ISBN-13: 978-90-9021117-6
Design and photo by Marilena Motta c° 2006.
Contents
Summary vii
Samenvatting xi
Nomenclature xiii
1 Introduction 1
1.1 Background and motivation . . . 1
1.2 Objectives of the research . . . 2
1.3 Outline of the thesis . . . 4
2 Theoretical background and literature overview 7 2.1 Introduction . . . 7
2.2 Reynolds averaged Navier-Stokes equations . . . 7
2.2.1 Near-wall behaviour . . . 10
2.3 The k − ε model . . . . 11
2.3.1 Transport equation for k . . . 11
2.3.2 Transport equation for ε . . . 12
2.3.3 Boundary conditions . . . 12
2.4 Wall functions . . . 12
2.4.1 Wall function assumptions . . . 13
2.4.2 Scaling of flow variables . . . 13
2.4.3 Standard wall function expressions . . . 13
2.5 The υ2− f model . . . . 14
2.5.1 Transport equation for υ2 . . . . 14
2.5.2 Elliptic relaxation function f . . . 14
2.5.3 Boundary conditions . . . 15
2.6 Large-eddy simulation . . . 15
2.6.1 Filtering . . . 16
2.6.2 SGS modelling . . . 16
2.7 Validation test cases . . . 17
2.7.1 Channel and pipe flows . . . 18
2.7.2 Backward facing step flow . . . 22
2.7.3 Two-dimensional wavy wall . . . 22
2.7.4 Axisymmetric impinging jet . . . 23
3 The ζ − f model 27
3.1 Introduction . . . 27
3.2 Reynolds stress model . . . 28
3.3 Elliptic relaxation model . . . 31
3.4 Derivation of ζ equation . . . 32
3.5 Differences between υ2 and ζ equations . . . . 34
3.6 Modified elliptic relaxation function . . . 36
3.7 Model coefficients . . . 37
3.8 Free-stream sensitivity . . . 39
3.9 Temperature field . . . 40
3.10 Remarks on the ζ-f model . . . 41
3.11 UMIST ϕ-f model . . . 42
3.12 Validation on test cases . . . 43
3.13 Discussion and conclusions . . . 48
4 Generalised Wall Function 55 4.1 Introduction . . . 55
4.2 Derivation of the generalised wall function . . . 56
4.3 Analysis of non-equilibrium terms . . . 60
4.4 Thermal wall function . . . 63
4.5 Notes on implementation . . . 65
4.6 Inclusion of additional effects . . . 65
4.7 Validation on test cases . . . 67
4.8 Discussion and conclusions . . . 70
5 Compound Wall Treatment 73 5.1 Introduction . . . 73
5.2 Rationale for compound wall treatment . . . 74
5.3 Compound wall treatment strategies . . . 79
5.4 Wall shear stress . . . 82
5.5 Turbulent kinetic energy dissipation rate . . . 83
5.6 Production of the turbulent kinetic energy . . . 85
5.7 Elliptic relaxation function . . . 86
5.8 Validation on test cases . . . 87
5.9 Discussion and conclusions . . . 93
6 Impinging jet in a cross-flow cooling of a heated cube 97 6.1 Introduction . . . 97
6.2 Description of the flow . . . 100
6.3 Computational aspects . . . 103
6.3.1 Computational mesh . . . 104
6.3.2 Computational costs . . . 108
6.3.3 Initial and boundary conditions . . . 110
6.3.4 Numerical and model details . . . 113
6.3.5 Validation of LES results . . . 114
CONTENTS v
6.4.1 Flow around the impinging jet . . . 120
6.4.2 Flow around the central cube . . . 128
6.4.3 Budgets of the turbulent stresses . . . 136
6.5 Temperature field . . . 145
6.5.1 Surface temperature . . . 145
6.5.2 Turbulent heat fluxes . . . 146
6.5.3 Budgets of the turbulent heat fluxes . . . 154
6.6 The ζ-f predictions . . . 158
6.6.1 Flow field . . . 158
6.6.2 Turbulent quantities . . . 158
6.6.3 URANS with the ζ-f model . . . 166
6.6.4 Temperature field . . . 166
6.7 Discussion and conclusions . . . 168
7 Conclusions and recommendations 171 7.1 Conclusions . . . 171
7.2 Recommendations for the future research . . . 174
A The ζ-f model and compound wall treatment 177 A.1 Flow equations . . . 177
A.2 Near-wall treatment . . . 178
A.3 Set of coefficients . . . 180
B Generic test flows results 181 B.1 Plane channel flow . . . 181
B.2 Backward-facing step . . . 185
B.3 Axisymmetric impinging jet . . . 185
B.4 Periodic wavy hill . . . 185
C Numerical method 189 C.1 Discretisation . . . 189 C.2 Convective terms . . . 190 C.3 Diffusive terms . . . 190 C.4 Source terms . . . 191 C.5 Time integration . . . 191 C.6 Pressure-velocity coupling . . . 191 C.7 Multiple materials . . . 191 Bibliography I Acknowledgement XI
Summary
Modelling and Simulation of Turbulence and
Heat Transfer in Wall-Bounded Flows
At present it is widely accepted that there is no universal turbulence model, i.e. no turbulence model can give acceptably good predictions for all turbulent flows that are found in nature or engineering. Every turbulence model is based on certain assumptions, and hence it is aimed at certain type of flows for which it recovers the effects that are occurring in that flow. Therefore the basic requirement for the development of a new turbulence model is the understanding of the fundamental fluid flow principles. The work presented in this thesis aims at improving the turbulence modelling applicable to the complex wall-bounded engineering flows and heat transfer on one side, and on the other it gives a detailed analysis of a complex wall-bounded flow and heat transfer with multiple flow phenomena.
The elliptic relaxation model of Durbin (1991) brings significant improvement in cal-culation of wall-bounded flows, in comparison to the standard k-ε model and its low-Re
variants. The υ2-f model solves the transport equation for the wall-normal Reynolds
stress component υ2, which brings in a limited way some features of the Reynolds stress
anisotropy, and the elliptic equation for the relaxation function f , which captures non-viscous wall effects. However, despite the fact that in many complex turbulent flows with
heat transfer the υ2-f model gives much better results than the k-ε model, it is still not
widely used for the computation of complex engineering flows because of its restrictions: the near-wall behaviour, the boundary condition for f and the mesh sensitivity.
In order to derive the turbulence model which is to be implemented into a com-mercial engineering CFD code, and which is suitable for the calculation of complex non-equilibrium wall flows using segregated solvers, Durbin’s elliptic relaxation approach with the integration to the wall was chosen for the starting point. As a way to overcome the
shortcomings of the υ2-f model, Popovac and Hanjali´c (2006a) and Hanjali´c and Popovac
(2004) presented a modification of Durbin’s eddy-viscosity model, called the ζ-f model,
in which the transport equation for the velocity scale ratio ζ = υ2/k is solved instead of
the equation for υ2. Another improvement that was elaborated is the application of the
SSG quasi-linear pressure-strain model in the equation for the elliptic relaxation func-tion, instead of using the basic IP model. As a consequence of this transformation the transport equation for the new velocity scale ratio, in comparison to its parent transport equation, appears to have better numerical stability.
Capturing and resolving near-wall effects is crucial for an accurate prediction of the characteristics of a wall-bounded turbulent flow. If the governing equations are solved all the way to the wall, then the model which is used has to include the near-wall effects: non-viscous wall-blocking and viscous damping. In this case exact boundary conditions for flow variables can be imposed on the solid boundary, and the accuracy in reproducing the local flow effects (separation, impingement etc.) depends on how accurate these effects are described by the model. The price one has to pay for this approach is a very fine mesh resolution near the wall, which can be a limiting factor in case of complex flow calculations.
However, when a coarse mesh is used so that the first near-wall cell is in the fully tur-bulent region, the majority of the local near-wall effects are not resolved. In this case wall functions, which relate the values of the variables in the cell centres with those at the wall through pre-integrated simplified expressions, are used to define the boundary conditions at the wall. The drawback of the wall function approach is that the simplifications used for their deriving exclude accurate capturing of general flow effects.
Popovac and Hanjali´c (2006a) and Popovac and Hanjali´c (2005a), following the argu-mentation of Craft et al. (2002), derived the generalised wall function in a more compact form. The idea was to integrate analytically the parabolised momentum equation, as-suming the linear variation of the turbulent viscosity throughout the near-wall region. Using the generalised wall function for complex flow computations on coarse meshes, non-equilibrium flow effects (pressure gradient, convection and transient effects) are recovered to some extend. But it is unlikely that any assumptions can be sufficiently accurate and reliable to fully recover all the near-wall effects of a non-equilibrium flow.
The most frequent flows of engineering relevance are bounded with solid walls of a complex topology that involve protrusions, cavities or other geometrical irregularities. In that case it is very difficult to control the mesh quality near the wall, hence it is hard to know in advance whether a near-wall cell in such a flow will be in the viscous sublayer, in the fully turbulent region or somewhere in between. Or viewed from another perspective, it is hard to know beforehand whether the integration to the wall or the pre-integrated wall functions can be applied for the calculation of such a flow. This problem is even more pronounced in a computation of flows which involve a strong acceleration and deceleration or a streamline curvature, where by the virtue of the flow the same cell size can fall either in the viscous or the fully turbulent zone. Therefore there is a need for the wall treatment which combines the wall limiting and the farfield definition of the boundary condition for the variable under consideration, given the computational mesh and the flow conditions.
SUMMARY ix The major quality of this wall treatment is that it significantly reduces the influence of the size or the distribution of the near-wall cells on the quality of the numerical predictions. Together with the ζ-f model this wall treatment is implemented into the commercial CFD package FIRE of AVL List, and set as the default turbulence model and the near wall treatment (Basara 2006; Brohmer et al. 2005; Tatschl et al. 2006).
Originating from the electronics industry applications, the LES study of the vortical structures and heat transfer in the case of the cooling of one in a series of wall-mounted cubes placed in an array on the wall of a channel flow by a round impinging jet is pre-sented. The geometry in this numerical analysis is defined by an in-line array of five cubes mounted on the wall of a plane channel. The central cube consists of the copper core and a layer of epoxy around the core, and only this cube is heated by specifying the constant temperature of the copper core. This cube is cooled by two mutually perpen-dicular streams: a channel flow parallel to the mounting wall, and a round impinging jet issued from the orifice which connects the channel with cubes and the chamber placed above it.
This LES analysis, presented by Popovac and Hanjali´c (2006b), followed the experi-ment reported by Tummers et al. (2005) and Flikweert (2005), in which detailed PIV,
LDA and liquid crystal measurements were given. The computational mesh had 8.5 ×106
cells, the subgrid scale stresses were calculated by the dynamic Smagorinsky model of Germano et al. (1991), for the convective scheme the central differencing was used, and the heat conduction trough the solid was calculated together with the convective heat transfer on the solid-fluid interface.
In addition to the analysis of the flow structures and their correlation with the heat transfer, the computation of fluid flow and the conjugate heat transfer for a jet in the cross-flow in such a complex geometry was also a challenging test for the ζ-f model with the near-wall treatment. It is very difficult to predict accurately multiple flow phenomena that are occurring in a single flow (impingement, separation, reattachment, cross-flow etc.), especially when the geometry is complicated. It is difficult even to generate the mesh which will meet either the requirements for the integration to the wall or for the wall functions.
Samenvatting
Modellering en Simulatie van Turbulentie
en Warmte Overdracht in Wandstromingen
Er is brede consensus in de literatuur dat er geen universeel turbulentiemodel is dat goede voorspellingen geeft voor alle type turbulente stromingen. Ieder turbulentiemodel is gebaseerd op bepaalde veronderstellingen waardoor het model in het algemeen goed bruikbaar is voor de stromingen waarvoor deze veronderstellingen voor gelden. Voor het ontwikkelen van een nieuw turbulentiemodel voor wandstromingen is het dus nodig om de belangrijke processen die optreden in een dergelijke stroming goed te begrijpen. Het werk gepresenteerd in dit proefschrift richt zich op het verbeteren van turbulentiemodellering voor stroming en warmteoverdracht in de buurt van een wand. Verder wordt de stro-ming en warmteoverdracht in de buurt van de wand geanalyseerd. Dit proefschrift kan worden opgedeeld in twee delen. In het eerste gedeelte (Hoofdstukken 3, 4 en 5) wordt gekeken naar de turbulentiemodellering, waarin de wandmodellering wordt beschreven. In het tweede gedeelte (Hoofdstuk 6) wordt gekeken naar de simulatie van stromingen die begrensd zijn door een wand.
Het elliptische relaxtiemodel van Durbin (1991) geeft een aanmerkelijke verbetering in de berekening van wandstromingen in vergelijking tot standaardmodellen zoals het
k-²-model en de laag-Reynolds varianten van dit model. Het υ2 − f model lost een
transportvergelijking op voor de wand-normale component van de Reynoldsspanning υ2.
Deze transportvergelijking weerspiegelt enigzins de anisotropie van de Reynoldsspanning. De elliptische vergelijking voor de relaxatiefunctie f beschrijft de niet-visceuze effecten bij
de wand. Ondanks dat het υ2− f model veel betere voorspellingen geeft voor turbulente
stroming en warmteoverdracht wordt het nog steeds niet veel gebruikt. De reden hiervoor is voornamelijk dat de correcte formulering van randvoorwaarden voor f erg lastig is. Als beginpunt voor het ontwikkelen van een turbulentiemodel dat goed te implementeren is in een commercieel software pakket, gebruiken we de elliptische relaxatie aanpak van Durbin. De nadelen van deze aanpak kunnen gedeeltelijk te niet worden gedaan door gebruik te maken van het ζ − f model, Popovac and Hanjali´c (2006a) and Hanjali´c
and Popovac (2004) waarbij een vergelijking voor ζ = υ2/k wordt opgelost in plaat van
een vergelijking voor υ2. Een andere verbetering is het gebruik van het SSG
quasi-lineare druk-snelheidsmodel in de vergelijking voor de elliptische relaxatie. Dit laatste is een alternatief voor het gangbare IP-model. Het gevolg hiervan is dat de nieuwe transportvergelijking betere numerieke eigenschappen heeft.
Het oplossen van de wandeffecten is van cruciaal belang voor een accurate voorspelling van de karakteristieken van een wand-gebonden turbulente stroming. Voor het geval dat de vergelijking wordt opgelost tot aan de wand, moet het juiste wand gedrag in het model worden verwerkt, zoals visceuze demping en blokkage. In dit geval kunnen exacte randvoorwaarden worden opgelegd op de vaste wand. Dit impliceert wel dat het rekenrooster erg fijn moet zijn en dit brengt de nodige kosten met zich mee. Als er een grof rekenrooster wordt gebruikt en het eerste rekenpunt in het turbulente gebied ligt, zal het merendeel van de wandeffecten niet worden opgelost. In dit geval is het nodig om gebruik te maken van wandfuncties. Popovac and Hanjali´c (2006a) en Popovac and Hanjali´c (2005a) gebruiken de argumentatie van Craft et al. (2002), om een gegeneraliseerde wandfunctie af te leiden. Het basisidee is hier om de parabolische impulsvergelijking te integreren waarbij verondersteld wordt dat turbulente viscositeit een linear gedrag vertoont in het wandgebied. Deze wandfuncties beschrijven ten dele niet-evenwichts effecten. Het is niet aannemelijk dat deze aanpak voldoende nauwkeurig is om alle effecten te beschrijven.
De meest voorkomende technische stromingen zijn stromingen in pijpen en kanalen, met complexe geometrien. Het is in dit geval erg lastig om de kwaliteit van het reken-rooster in de buurt van de wand te voorspellen. Het is dus niet duidelijk in welk gebied wel en in welk gebied geen wandfunctie nodig is. Dit probleem wordt nog lastiger wanneer er sprake is van versnellende of vertragende stromingen. Er is dus behoefte aan een benader-ing van het wandgebied die breed toepasbaar is. Een eerste aanzet voor een dergelijke randvoorwaarde staat beschreven in Popovac and Hanjali´c (2006a) en Popovac and Han-jali´c (2005b). Hierin wordt de randvoorwaarde gegeven ongeacht de stromingscondities en het rekenrooster bij de wand. De methode vereenvoudigt tot een standaard inte-gratieprocedure bij de wand, in het geval dat de eerste rekencel in de visceuze sublaag ligt. Als dit niet het geval is wordt er een wandfunctie gebruikt. Het samenstellen van de twee voorwaarden is gedaan met behulp van het criterium van Kader (1981). Ondanks enige verschillen kan deze samenstelling gezien worden als een vereenvoudiging van de “algebraic elliptic blending” zoals voorgesteld door Manceau and Hanjali´c (2002). Het grote voordeel van de wandbehandeling is dat de oplossing minder gevoelig wordt voor de kwaliteit van het rekenrooster bij de wand. Dit wandmodel is samen met het ζ − f model gemplementeerd in het commercile software pakket FIRE van AVL List en is nu het standaardmodel van dit pakket (Basara 2006; Brohmer et al. 2005).
Een probleem dat zijn oorsprong in de electronicaindustrie heeft, is de LES van vor-tices en warmteoverdracht rond een rij van kubusvormige obstakels, die gekoeld worden door middel van een straal. De rekengeometrie bestaat uit vijf kubusjes die bevestigd zijn op de wand van een kanaal. De middelste kubus is gemaakt van koper waaromheen een mantel van epoxy zit. Alleen deze kubus wordt verwarmd door de temperatuur van het koper voor te schrijven. De kubus wordt gekoeld door twee stromingen: een kanaal-stroming parallel aan de wand en door een rondje straal die tegen de bovenkant van de kubus stroomt. De LES-berekeningen zijn gepresenteerd in Popovac and Hanjali´c (2006b) en vergeleken met de experimentele resultaten van Tummers et al. (2005) en Flikweert (2005), waarin gedetailleerde PIV en LDA metingen, alsmede metingen met vloeibare
kristallen zijn gegeven. Het rekenrooster bevatte 8.5 · 106 cellen, voor het subgrid model
Nomenclature
Mean flow quantities
Symbol Description Unit
CU = ρDUDt +∂x∂Pi parabolized momentum convection kg/m2s2
CΘ = ρcpDΘDt − QΘ parabolized temperature convection kgK/m3s
E power density −
fs natural frequency of structures 1/s
˙
m mass flux kg/s
P pressure kg/m s2
qw wall heat flux kg/s2
Sij = 12 ³∂U j ∂xi + ∂Ui ∂xj ´ rate-of-strain tensor s−1 t time s
U , V , W streamwise/wall-normal/spanwise velocity components m/s
Ui, Uj, Uk cartesian velocity components m/s
Uz, Ur, Uθ cylindrical velocity components m/s
|U| velocity magnitude m/s
αm mass flux ratio −
αv velocity ratio −
ατ wall shear stress ratio −
Θ temperature K
τw wall shear stress kg/m s2
Ωij = 12 ³∂U j ∂xi − ∂Ui ∂xj ´ rate-of-rotation tensor s−1 |ω| vorticity magnitude s−1
ψU = 1 − U κuCUyτ non-equilibrium momentum correction factor −
ψΘ = 1 − ρκuP rτtCcpΘ∆Θy non-equilibrium temperature correction factor −
Physical properties
Symbol Description Unit
a speed of sound m/s
cp specific heat capacity m2/s2K
K thermal activity −
P r Prandtl number −
Sc Schmidt number −
α thermal diffusivity m2/s
µ dynamic molecular viscosity kg/ms
ν kinematic molecular viscosity m2/s
ρ fluid density kg/m3
λ thermal conductivity kg m/s3K
Turbulent quantities
Symbol Description Unit
A2 second stress-anisotropy invariant −
aij = uikuj − 23δij Reynolds stress anisotropy tensor −
Dij Reynolds stress diffusion tensor m2/s3
f elliptic relaxation function m2/s3
I = u
U turbulent intensity −
k turbulent kinetic energy m2/s2
L = k3/2
ε turbulent length scale m
L characteristic length scale m
P production of turbulent kinetic energy m2/s3
Pij Reynolds stress production tensor m2/s3
P Jayatilleke function −
p fluctuating pressure kg/m s2
P rt turbulent Prandtl number −
Rij Reynolds stress redistribution tensor m2/s3
T = kε turbulent time scale s
T characteristic time scale s
U+ = U
uτ normalised velocity −
ui fluctuating velocity component m/s
uiuj Reynolds stress component m2/s2
uk= c1/4µ k1/2 turbulent kinetic energy velocity scale m/s
uτ =
q τw
ρ friction velocity m/s
υ2 wall-normal velocity scale m2/s2
X = 2νef k ∂k ∂xk ∂ζ ∂xk cross-diffusion of ζ m 2/s2 y+ = yuτ
ν normalised distance to the wall −
αt= P rνtt = ρcλtp turbulent thermal diffusivity m2/s
Γ blending factor exponent −
ε dissipation rate of turbulent kinetic energy m2/s3
εij Reynolds stress dissipation tensor m2/s3
φ generic quantity −
λt turbulent thermal conductivity kg m/s3K
µt = ρνt turbulent dynamic viscosity kg/ms
NOMENCLATURE xv
Θ+ = ∆Θρcpuτ
qw normalised temperature −
θ fluctuating temperature K
θui turbulent heat flux Km/s
ζ = υ2
k normalized wall-normal velocity scale −
Geometrical quantities
Symbol Description Unit
D nozzle/orifice diameter m
h channel half-height, height of the step/hill/cube m
H channel height, nozzle-to-plate distance, backstep height m
n normal unit vector m
S surface area m2
V volume m3
x, y, z streamwise/wall-normal/spanwise coord. components m
xi, xj, xk cartesian coordinate components m
z, r, θ axial/radial/azimuthal coordinate components m; m; rad
αr radius ratio −
∆ filter width m
Constants
Symbol Description
aT = 0.6 time scale realisability constant
C1 = 1.4, C2 = 0.3 redistribution model constant
CL = 0.36, Cη = 85 length scale constant
CT = 6.0 time scale constant
CS = 0.1 Smagorinsky constant
Cε1 = 1.4, Cε2 = 1.9 dissipation equation constants
Cµ = 0.22 Durbin’s turbulent viscosity constant
cµ = 0.09 turbulent viscosity constant
EU = 8.3 log-law additive constant
P rt = 0.9 turbulent Prandtl number
κ = 0.41 von Karman constant
σk= 1.0 turbulent Prandtl number for k
σε= 1.3 turbulent Prandtl number for ε
Dimensionless numbers
Symbol Description
CF L = ∆tU∆x Courant-Friedrich-Levy number
cf = ρU2τw2 friction factor
M = Ua Mach number N u = qwD ∆Θλ Nusselt number Re = U L ν Reynolds number St = qw ∆ΘρcpU Stanton number Sts= fsUD Strouhal number
Abbreviations and acronyms
Symbol Description
1-D one-dimensional
2-D two-dimensional
3-D three-dimensional
BiCG bi-conjugate gradient
CDS central differencing scheme
CFD computational fluid dynamics
CG conjugate gradient
CGS conjugate gradient squared
CWT compound wall treatment
DNS direct numerical simulation
ERM elliptic relaxation model
IP isotropisation of production
ItW integration to the wall
LRR Launder, Reece, Rodi
LES large-eddy simulation
MinNOx minimisation of NOx emissions
RANS Reynolds averaged Navier-Stokes
RI return to isotropy
RSM Reynolds stress model
SGS subgrid scale
SMC second moment closure
SSG Speziale, Sarkar, Gatski
Chapter 1
Introduction
The first chapter gives an introduction to this thesis. It presents the motivation for this work, explains the source of the problems whose solution is sought and the solutions which are offered. At the end of this chapter the outline of the thesis is given.
1.1
Background and motivation
The developments within the European Union project MinNOx, from which this re-search has been sponsored, explain clearly the motivation for the work presented in this thesis. The project gathered some of the biggest names of the European automotive industry (Daimler-Chrysler, Ford and Volvo) on one side, and Europe’s leading univer-sities (Chalmers University of Technology G¨oteborg, King’s College London, Universit¨at Stuttgart and of course Technische Universiteit Delft) on the other. The aim of the project was the development of a physical model for accurate predictions of heat transfer, ve-locity and temperature field for the flow inside the cylinder of an internal combustion engine. This model was to be validated by comparing numerical predictions with the measurement results provided by the experimentalist partners in the project.
In the first half of the project the TU Delft team presented a new turbulence model with improved accuracy, stability and robustness, which has been verified on several generic non-equilibrium test cases. By that time the experimentalists did not produce any conclusive result because they were still struggling to make their experimental setup run. Soon after the development at TU Delft another partner in the project, AVL List Graz, implemented this model into their commercial code FIRE and obtained the first results on a real engine numerical setup – long before there was any experimental reference data to compare with.
Apart from the time table, a very important detail in the MinNOx story concerns the money. The local computer that I was using was an eight years old defected Pentium II with a monitor the colours of which were fading, whereas most of the individual com-ponents used by the experimentalists for their engine measurements costed much more than that single computer. Accordingly, the budget of the experimentalist partners in the project was twice as big as ours.
This example gives a clear indication how important computational fluid dynamics (CFD) is today in different fields of engineering. The ability to produce and analyse
many results within a reasonable time frame and without a large financial investments is a very strong argument in favour of this engineering approach. It introduces a signif-icant speed-up in the design and optimisation process, thus improving the research and development in actual engineering. The numerical approach is particularly attractive for the automotive industry and alike where experimental research is both very expensive and time consuming, and where it is not certain that a full range of information will be achieved. For engineering applications the most widely used approach is the Reynolds averaged Navier-Stokes (RANS) solution of the fluid flow problem.
Apart from engineering applications, CFD also has a great potential in the research of the physics of fluid flows. For this purpose, however, advanced numerical techniques such as large-eddy simulation (LES) have to be used. Over the years LES has matured into a very powerful numerical method that can tackle even very complex fluid flow problems, but the user has to have sufficient knowledge and experience to set up the computation properly and to analyse the results correctly. Still, the major obstacle for total dominance of the numerical over the experimental approach is the reliability of the numerical results – experiments are still needed to verify numerical predictions. The only way to make further advances in numerical fluid flow and heat transfer predictions is to follow fully the physical rationale, and to search for the relationship between the flow structures and the effects that they cause, in order to be able to describe those phenomena mathematically. The work presented in this thesis follows these guidelines.
1.2
Objectives of the research
The first part of the thesis presents the new turbulence model developed to meet the requirements of the MinNOx project. The basic prerequisite for the accurate predictions of the NOx emissions of internal combustion engines is the accurate prediction of the temperature field in the cylinder of the engine. Naturally, the accurate and reliable velocity field is necessary for the accurate temperature field predictions. Therefore, in order to be able to design the engines towards the reduction of the NOx emissions, which was the main idea of the MinNOx project, it is important to have a physical model that gives an accurate solution for the fluid flow inside the cylinder. This was the task of the TU Delft team within the MinNOx project.
The proposed model was to be simple so that it can be implemented into a commercial CFD code, robust (the antonym of fragile) so it can handle the computation of such a complex flow, and yet be able capture the most important effects occurring in the in-cylinder flow: the near-wall and transient effects, strong pressure gradients, impingement, separation, reattachment, swirl and similar. For all but the last requirement the standard k-ε model performs satisfactorily: it is a simple, fast and robust turbulence description. The biggest disadvantage of k-ε is its isotropic eddy viscosity for all shear components, which locally produces erroneous prediction of the fluid flow characteristics especially in near wall flows .
1.2. Objectives of the research 3 is of the Poisson type), yet it includes most of the effects that occur in wall-bounded flows, in particular both the viscous and the non-viscous wall-blocking effect. The main drawback of the original elliptic relaxation model that was considered is its very stiff boundary condition, which makes it unstable for the complex flow computations. In addition to that, in its original form that model requires a computational mesh suited for integration to the wall, which is a prohibitive requirement in complex engineering applications.
The objective of the first part of the work, therefore, is to bring improvements in numerical predictions by using a higher level turbulence modelling, as compared to the predictions obtained with standard k-ε model. The ultimate requirement for the new model is that all the improvements should be neither at the expense of the computational stability and robustness, nor the quality and size of the mesh. For that, the chosen higher level turbulence model has to be kept as simple as possible, and the appropriate near-wall treatment has to accommodate the correct wall boundary conditions irrespective of the mesh near the wall. The final outcome of the work performed within the MinNOx project (a)
(c)
(b)
(d)
(the new turbulence model and the near-wall treatment) has been successfully used for the purpose which was the original motivation of the project: Fig. 1.1 shows the engine calculations by AVL’s code FIRE with the turbulence model and wall treatment proposed by this work.
The second part of the thesis presents a physical insight into a wall-bounded turbulent flow with heat transfer. Contrary to the work within the MinNOx project, which was oriented towards engineering applications (that implied searching for pragmatic solutions in order to obtain reasonably good result), in the last part of the thesis the focus is on the investigation of the physics of the particular wall-bounded flow, its structures, dynamics and the heat transfer effects. In particular, the correlation between the local flow structures and heat transfer was analysed. For this analysis large-eddy simulation has been performed in a selected complex flow, since LES is well established advanced numerical research tool. The flow that was examined is a round jet in a cross-flow impinging on a wall-mounted cube, and it features several flow phenomena, such as impingement, cross-flow, impingement and separation. Because of its complexity, this was also a challenging test case for the newly developed turbulence model.
1.3
Outline of the thesis
This thesis is organised into seven chapters and three appendices. The chapters are ordered chronologically, according to developments in the research. Each chapter begins with a short abstract written in cursive, and all chapters except the first and the last (there it was simply redundant) finish with conclusions and discussion.
According to the time schedule of the MinNOx project, the first period was reserved for a survey of related literature and testing of existing fluid flow modelling strategies. Therefore after this introductory Chapter 1, the overview of existing published work, which was used as a basis for this research, is briefly presented in Chapter 2. This includes the basics of turbulence modelling, an overview of some well known turbulence models and wall functions, and the validation test cases together with their flow features. The turbulence models presented include the most widely used k-ε model with wall functions
and Durbin’s υ2-f model, as well as the principles of large-eddy simulation. Instead of
the joint conclusion section, in this chapter at the end of each thematic entity there is an unindented paragraph of conclusions or final remarks.
1.3. Outline of the thesis 5 As an extension of the investigation of the near-wall modelling strategies, the devel-opment of a new wall function and near-wall treatment was the task planned for the third and the last year of MinNOx project. Extending the idea of integration of the parabolised momentum equation in the near-wall region, Chapter 4 gives the derivation and the validation of the simplified generalised wall function which takes into account the non-equilibrium effects. The improvement presented in this chapter is the inclu-sion of transient effects and the derivation of a compact expresinclu-sion which enables simple implementation in the standard log-law manner.
For turbulence models with integration to the wall the mesh resolution in the near-wall region has to be very fine, while for complex flow calculations (i.e. in the vast majority of engineering flow cases) the computational mesh is mostly suited for the wall function approach. Therefore, having discussed the turbulence model for the integration to the wall in the chapters above, Chapter 5 explores the possibilities of using that model on the meshes that could locally be both coarse as for the wall function, but also sufficiently fine for the integration to the wall. For that reason an effort was made to treat the buffer region properly. The line that was followed is to find a method which will appropriately combine the integration to the wall and the wall function approach, depending on the local flow parameters and the given mesh, having in mind the simplicity and robustness of the method.
In the last year of the research the focus was on the investigation of the physical aspects of some challenging case of fluid flow with conjugate heat transfer. The selected case includes a jet in cross-flow cooling the central of five in-line wall-mounted cubes by impinging onto its top surface. The computation of this demanding flow case was performed using LES, so that the skill and proficiency is demonstrated in performing a complex flow simulation using this powerful advanced technique and analysing the ob-tained data. Namely, it is now widely regarded that a researcher or engineer in CFD should sufficiently master both LES and RANS, in order to be able to select the ap-propriate tool for the apap-propriate flow problem. This has been recognised also by the CFD vendors, who are now offering LES as a standard option in their CFD packages. However, LES brings some concerns about accuracy and reliability related to the grid resolution (much more than in RANS), and proper interpretation of LES data. Therefore LES is very potential engineering tool that offers much more information about the flow structure and its unsteady dynamics, but only if it is correctly used.
Having these aspects in mind, Chapter 6 gives the analysis of the flow phenomena for an impinging jet in a cross-flow and conjugate heat transfer from a wall-mounted cube obtained with LES. It also discusses the results, the mesh requirements and the numerical details of LES. Additionally, the calculation of this complex flow with the new turbulence model and the near-wall treatment was presented in order to demonstrate the capability of the model to capture multiple complex fluid flow phenomena.
the numerical method implemented into the code that was used for all the computations presented in this thesis.
Chapter 2
Theoretical background and
literature overview
This is an overview chapter which briefly presents the published work that was used as the basis for the present research. That includes the basic elements of turbulence modelling, the Reynolds averaged Navier-Stokes equations, the standard
k-ε and υ2-f turbulence models, the large eddy simulation and the generic test cases
featuring different flow effects. This chapter also gives references to the literature where a given work has been published.
2.1
Introduction
The research presented in the first part of this thesis, aiming at the improvement of the existing near-wall RANS modelling strategies applied to complex turbulent flows, originates from the work of a number of researchers who published their findings on wall-bounded turbulent flow modelling. This chapter gives an overview of the work relevant to this research, including the basics of the Reynolds averaged Navier-Stokes approach and large-eddy simulation method. At the end, the generic test cases and their flow features used for validating the model will be discussed.
Being not of the prime interest of this research on one hand, and being well docu-mented in the literature on the other, the description given here brings only the aspects of each work which is of the importance for this research. More details can be found in Pope (2000), Durbin and Pettersson Reif (2000) and other turbulence textbooks.
2.2
Reynolds averaged Navier-Stokes equations
The instantaneous turbulent flow and heat transfer in an incompressible Newtonian fluid, considered in this thesis, is governed by the set of equations commonly referred to as the Navier-Stokes equations. Those are the momentum equation:
∂uei ∂t +uej ∂uei ∂xj = − 1 ρ ∂pe ∂xi + ν ∂ 2ue i ∂xi∂xj (2.1) the thermal energy (enthalpy) equation:
x x3 1 P1 s2 s3 s1 P3 P2 n n P s 0 r 0 V x v z y x2 P
Figure 2.1: Control volume of arbitrary shape and the notation used.
∂ eθ ∂t +uei ∂ eθ ∂xi = α ∂ 2eθ ∂xi∂xi (2.2) and the continuity equation:
∂uei ∂xi
= 0 (2.3)
where uei is the instantaneous velocity, p is the instantaneous pressure, ee θ is the
instan-taneous temperature, ρ is the fluid density, ν is the molecular viscosity, α is the thermal
diffusivity of the fluid, t is the time and xi are the space coordinates. The thermal
dif-fusivity α is related to the thermal conductivity λ as α = λ/ρcp, with cp being fluid’s
specific heat capacity.
Here the convention, to be used throughout this thesis, is introduced: when the index
notation is used, xi and correspondingly Ui(i = 1, 2, 3) are referring to the components of
the Cartesian coordinate system. However, when the coordinate components are denoted as x, y and z, with corresponding velocity components U , V and W , they are referring to the coordinate system originating at the wall, defined by the streamwise, the wall-normal and the spanwise flow direction respectively. This convention is illustrated in Fig. (2.1). Having in mind important characteristics of turbulent flows, such as complexity, un-steadiness, irregularity, three-dimensionality and a wide range of scales (Tennekes and Lumley 1976), it is clear that a turbulent flow is fully described only if Eq. (2.1) and Eq. (2.2) resolve the smallest scales occurring in the flow, i.e. the Kolmogorov time and length scale (Kolmogorov 1941):
TKol = ν 1/2 ε1/2 LKol = ν 3/4 ε1/4 (2.4)
2.2. Reynolds averaged Navier-Stokes equations 9 kinetic energy ε. The method where Navier-Stokes equations are directly applied to the fluid flow problem is called the direct numerical simulation (DNS).
For most of the practical problems, however, it is too expensive to perform DNS even on the fastest supercomputers. This is why a much simpler statistical approach is sought. Reynolds (1894) proposed the averaging of the flow quantities (denoted with the overline):
the instantaneous flow variables, denoted with the tilde (uei, eθ and p) are decomposede
into the long-time mean, denoted with the capital letters (ui = Ui, θ = Θ and p = P ),
and the fluctuating part, denoted with small letters (ui, θ and p). The averaged value φ
is obtained from the instantaneous field of the given variable eφ as follows:
φ = 1 tav Z t+tav t e φ(t)dt (2.5)
where tav is the averaging time interval, longer than any time scale occurring in the flow.
The fluctuating part is then the difference between the instantaneous and the averaged
value φ = eφ − φ.
Applying the averaging operator given by Eq. (2.5) to the decomposed Navier-Stokes equations leads to the Reynolds averaged Navier-Stokes equations (RANS):
∂Ui ∂t + Uj ∂Ui ∂xj = − 1 ρ ∂P ∂xi + ν ∂ 2U i ∂xi∂xj − ∂uiuj ∂xj (2.6)
where the second moment tensor uiuj is called the Reynolds stress or the turbulent stress,
and it arises from averaging of the non-linear terms in Eq. (2.1). The Reynolds stress tensor represents the effects of turbulence, and it is a new unknown variable.
By repeating the same averaging procedure with Eq. (2.2), the RANS equation for temperature is obtained: ∂Θ ∂t + Ui ∂Θ ∂xi = ν P r ∂2Θ ∂xi∂xi − ∂θui ∂xi (2.7) where Prandtl’s analogy is used to express the thermal diffusivity α through the molec-ular viscosity ν as α = ν/P r, with P r being the molecmolec-ular Prandtl number of the fluid.
The tensor θui is a new unknown variable called the turbulent heat flux, which is the heat
transfer counterpart of the turbulent stress.
In order to close the set of equations given by Eq. (2.6), Eq. (2.7) and Eq. (2.3),
tensor uiuj and vector θui have to be expressed in terms of known variables through a
certain physical rationale – they have to be modelled. By analogy to the viscous stress, the common way to model the turbulent stress and the turbulent heat flux (Boussinesq 1877)
is to introduce the turbulent or the eddy viscosity νt, hence the name for this approach
the eddy viscosity model:
−uiuj = νt µ ∂Ui ∂xj +∂Uj ∂xi ¶ − 23kδij (2.8)
and in the case of the turbulent heat flux the eddy diffusivity model reads:
where, again from the Prandtl’s analogy, the turbulent thermal diffusivity αt (which is
related to the turbulent thermal conductivity λtas αt= λt/ρcp), is expressed through the
turbulent molecular viscosity νt and the turbulent Prandtl number P rt as αt= λt/P rt.
Analogous to the definition of the molecular viscosity in the kinetic theory of gases, the turbulent viscosity is defined by the characteristic velocity scale U and the length scale L. For the molecular viscosity, however, the velocity and the length scale are properties of the fluid, related to the mean free path and the average molecular speed (Prandtl 1925; Batchelor 1967). On the other hand, for the turbulent viscosity these two scales are properties of the flow, and therefore they have to be calculated for every specific flow. Obviously the Reynolds averaging cannot close the set of equations given by the
Navier-Stokes and the continuity equation, since by averaging new unknowns , uiuj and
θui, appear. But in order to have an idea on how to model Reynolds stresses, it is useful
to derive the transport equations for uiuj (the analysis for θui is analogous, but it will not
be given here). The final equation is obtained by subtracting Eq. (2.6) from Eq. (2.1),
the resulting equation for ui is multiplied by uj, the indices interchanged and equations
summed: ∂uiuj ∂t | {z } time variation Vij + Uk ∂uiuj ∂xk | {z } convection Cij = −∂ukuiuj ∂xk | {z } turb. transport Dt ij −ujuk ∂Ui ∂xk − u iuk ∂Uj ∂xk | {z } production Pij + ν∂ 2u iuj ∂xk | {z }
mol. dif f usoin Dν ij −1 ρ µ uj ∂p ∂xi + ui ∂p ∂xj ¶ | {z } redistribution Πij − 2ν∂ui ∂xk ∂uj ∂xk | {z } dissipation εij (2.10) and under each term its meaning is given in short.
2.2.1
Near-wall behaviour
In the immediate vicinity of a solid boundary the variation of the fluctuating velocities in the wall-normal direction can be analysed from its expansion in Taylor series in the wall-normal direction. For fixed x, z and t it reads:
u(y) = a1+ b1y + c1y2+ . . .
v(y) = a2 + b2y + c2y2 + . . .
w(y) = a3 + b3y + c3y2 + . . .
(2.11) where a’s, b’s, c’s and so on are the coefficients statistically independent on x, z and t.
From the no-slip boundary condition (u = v = w = 0 at y = 0) it follows that
a1 = a2 = a3 = 0. At the wall the derivatives (∂u/∂x)y=0 and (∂w/∂z)y=0 are also zero,
since u and w are zero for all x and z. Hence, from the continuity equation (∂v/∂y)y=0
must also be zero. This implies zero value for the b2 constant, b2 = 0. In order words, the
2.3. The k − ε model 11 From Eq. (2.11) it is easy to obtain the near-wall behaviour of the Reynolds stress tensor components: uu(y) = b2 1y2 + c21y4 + . . . vv(y) = c2 2y4 + . . . ww(y) = b2 3y2 + c23y4 + . . . uv(y) = d4y3 + . . . (2.12)
i.e. uu and ww vary with the distance to the wall as y2, whereas uv and vv vary faster
(y3 and y4 respectively). The importance of this finding will be clear later on.
Presented here is the starting point for the analysis of RANS turbulence models (more can be found in many turbulence textbooks, Durbin and Pettersson Reif (2000) and Pope (2000) among others). The meaning of the terms in the Reynolds stress transport equation and the way to model them is given in Chapter 3.
2.3
The
k − ε model
This is the most widely used high Reynolds number turbulence model, first postulated
by Hanjali´c (1970)1, whose basic idea is to define the time and the velocity scale through
the turbulent kinetic energy k and its rate of dissipation ε. On dimensional grounds one
can see that the time scale is T = k/ε, while the velocity scale is U = k1/2. Thus the
definition of the turbulent viscosity reads:
νt = cµ
k2
ε (2.13)
where the most common value for the turbulent viscosity constant is cµ= 0.09.
2.3.1
Transport equation for
k
Having in mind the definition of the turbulent kinetic energy, k = 12uiui, the transport
equation for k can be derived from the Reynolds stress equation by taking half the trace of Eq. (2.10) for i = j: ∂k ∂t + Uj ∂k ∂xj = P − ε + ∂ ∂xj · (ν + νt) ∂k ∂xj ¸ (2.14) where the turbulent transport (triple correlation) and the pressure redistribution in Eq. (2.10) are jointly modelled by a simple gradient diffusion hypothesis.
1
2.3.2
Transport equation for
ε
Unlike k, the transport equation for ε is not systematically derived, but it is taken to have the same form as that for k (Eq. 2.14), only scaled with the turbulent time scale T = k/ε. The argumentation is that the dynamics of ε follows the dynamics of k.
∂ε ∂t + Uj ∂ε ∂xj = Cε1P − Cε2ε T + ∂ ∂xj ·µ ν + νt σε ¶ ∂ε ∂xj ¸ (2.15)
with the empirical coefficients Cε1 = 1.44 and Cε2 = 1.92 added, together with the
turbulent Prandtl number σε = 1.3. In case of integration to the wall, Cε1 is reformulated
in order to account for the near-wall anisotropy (see Table 3.2).
2.3.3
Boundary conditions
From the no-slip boundary condition, given by Eq. (2.12) it follows that the value for
turbulent kinetic energy at the wall is kw = 0. The boundary condition for ε follows from
the near-wall analysis of the transport equation for k, given by Eq. (2.14):
εw = 2ν lim y→0 µ k y2 ¶ (2.16) However, the transport equations in the standard k-ε model are not solved up to the wall. Instead, the boundary conditions are imposed in the near-wall cell centres by certain pre-defined expressions, as explained in the following section.
The k-ε model is very simple and therefore is widely used for engineering flow calcu-lations. In case of complex turbulent flows, however, two flow variables, k and ε, are simply not sufficient to describe accurately the flow characteristics. The major source of k-ε model inaccuracy is the isotropic definition of the turbulent viscosity and the deficiency of the modelled transport equation for ε.
2.4
Wall functions
In order to integrate Eq. (2.14) and Eq. (2.15) all the way up to the solid boundary, they have to be modified to account for the near-wall effects. But these modifications ruin the simplicity of the model, and what is more serious, these modifications are realized through empirical damping functions. Therefore Patankar and Spalding (1970) and later Launder and Spalding (1974) elaborated on a method in which the transport equations are not solved in the near-wall region, but the pre-integrated values of the flow variables are imposed in the near-wall cell centres.
Of course, this requires the near wall cell to be located outside the viscous layer (above
≈ 30 in wall units defined below) yet still below y ≈ 2δ99, where δ99 is the height where
2.4. Wall functions 13
2.4.1
Wall function assumptions
In order to derive the algebraic wall function expressions, certain assumptions on the flow characteristics have to be made. The first assumption for standard wall functions is that in the near-wall region the turbulent shear stress is equal to the wall shear stress
−ρuv = τw. Furthermore, it is also assumed that the turbulent shear stress is proportional
to k, i.e. that relation |uv|/k = c1/2µ holds. It is also assumed that the flow in the
near-wall region is in the local energy equilibrium, for which the equality of the production and dissipation P = ε is valid.
2.4.2
Scaling of flow variables
The wall function concept relies strongly on the scaling of the flow variables, i.e. on the normalisation with the appropriate velocity scale U and the length scale L = ν/U which yields a universal behaviour for different flow conditions. For the wall bounded flows it is common to use the velocity scale that follows from the assumption of the
shear stress balance in the wall vicinity U = uτ = |uv|1/2 =
p
τw/ρ, which is called the
friction velocity. The normalisation with uτ gives the wall or plus units, denoted with
the superscript “+”, but it can lead to singularity when τw naturally goes to zero (e.g. at
flow reattachment). In order to avoid this problem, the velocity scale definition through k is introduced, following the assumption of constant normalised shear Reynolds stress
U = uk = uv1/2 = c1/4µ k1/2. The normalisation with uk gives the star units, denoted with
the superscript “∗”. Depending on the specific flow conditions under consideration, other velocity scales expressed through additional flow parameters can be defined (Kaltenbach 2003; Shih et al. 2003).
2.4.3
Standard wall function expressions
Using the assumptions introduced above, the following wall function expressions are derived: U∗ = 1 κln (EUy ∗) dU dy = uk κy (2.17)
for the velocity, and for the dissipation and production:
εp =
u3 k
κyp = Pp
(2.18)
where κ = 0.41 is the Von K`arm`an’s constant, EU = 8.9 is the log-law additive constant,
and y∗ is the distance from the wall normalised with u
k = c1/4µ k1/2.
analysis of the wall functions together with the improvement of their performance in non-equilibrium flow conditions is given in Chapter 4. Chapter 5 examines the possibility to combine exact boundary conditions with wall function expressions.
2.5
The
υ
2− f model
Durbin (1991) noted that a more appropriate velocity scale for the definition of the turbulent viscosity is the wall-normal component of the Reynolds stress tensor vv instead of k. This, in fact, follows from the modelled equation for the full stress tensor when it is truncated to the single “source–sink” form close to the wall. In order to make use of this finding for an arbitrary geometry, this variable is treated as a scalar which resembles
the stress component normal to the streamlines. Denoting this new scalar as υ2, the
definition of the turbulent viscosity reads:
νt= Cµ
υ2k
ε (2.19)
with the adopted value for the turbulent viscosity constant Cµ = 0.22.
2.5.1
Transport equation for
υ
2By putting i = j = 2 in Eq. (2.10) and replacing vv = υ2, the transport equation for this
new velocity scale is obtained:
∂υ2 ∂t + Uj ∂υ2 ∂xj = kf − ε υ2 k + ∂ ∂xk " (ν + νt) ∂υ2 ∂xk # (2.20) where the pressure-velocity coupling term and the difference between the inhomogeneous and the homogeneous part of ε are treated jointly as a product of k and the elliptic relaxation function f . This term is crucial for correct representation of the Reynolds stress anisotropy effects.
2.5.2
Elliptic relaxation function
f
Arguing that the non-viscous wall effects are of the elliptic nature (Durbin 1993), an elliptic relaxation equation is constructed for the pressure-strain correlation to model non-viscous wall blocking effects:
L2∇2f − f = C1− 1 T Ã υ2 k − 2 3 ! − C2P ε (2.21)
2.6. Large-eddy simulation 15 T = max " k ε, CT µ ν3 ε ¶1/2# L = CLmax " k3/2 ε , Cη µ ν3 ε ¶1/4# (2.22)
and the coefficients used C1 = 1.4, C2 = 0.45, CT = 6, CL= 85 and Cη = 0.36.
2.5.3
Boundary conditions
The boundary conditions for k and ε are given in Eq. (2.16), and Eq. (2.12) shows that
the at the wall υ2
w = 0. The boundary condition for f , obtained from the near-wall
analysis of Eq. (2.20), reads:
fw = −20ν2lim y→0 Ã υ2 εy4 ! (2.23) This boundary condition is very stiff, and it can cause numerical instabilities, especially for the complex flow calculations.
Although for wall-bounded flows the υ2-f model, compared to the standard k-ε,
per-forms significantly better at relatively low price (only two extra equations, where f is not a full transport equation but an elliptic one), k-ε remained the most widely used
turbulence model for engineering applications. The reason why υ2-f did not overshadow
k-ε is the requirement for a high mesh resolution near the wall, and even more its high sensitivity to the near-wall mesh quality due to its boundary conditions. More details on
the υ2-f model (the rationale and coefficients), together with more robust modification
of this model is presented in Chapter 3.
2.6
Large-eddy simulation
Unlike in RANS, where the Navier-Stokes equations are averaged over a certain time interval, in large-eddy simulation (LES) the equations are filtered over a small region of space by applying a predefined spatial filter. The filtered equations are resolving the flow variables of the scale larger than the filter width, and modelling those on the smaller scale. When filtering is applied to the Navier-Stokes equations, the following equation is obtained: ∂ui ∂t + uiuj ∂xj = − 1 ρ ∂p ∂xi + ν ∂ 2u i ∂xi∂xj − ∂τSGS ij ∂xj . (2.24)
and as in the Reynolds averaging, the nonlinear terms produce new unknowns τSGS
ij ,
which are in LES context referred to as the subgrid scale Reynolds stresses (SGS):
These subgrid stresses play a similar role in LES as the Reynolds stresses do in RANS. However, they represent a much smaller portion of the total turbulent stress and kinetic energy than in RANS, and thus the model accuracy in LES is not as important as in RANS. Furthermore, the small scales which are modelled in LES tend to be more isotropic, therefore the SGS model can be much simpler than the RANS model.
For the LES of wall bounded flows the estimated number of cells is related to the
Reynolds number as Re9/4. This relationship is dictated by the small structures of the
turbulent flow in the wall vicinity. Thus, to accurately resolve the structures in the near wall region with LES, the first grid point in the wall-normal direction must be
located at y+ ≤ 1, and the grid spacing in streamwise and spanwise direction must be
∆x+ ≈ 50 − 100 and ∆z+ ≈ 15 − 40 respectively (Piomelli 1997; Reynolds 1990). For
numerical reasons it is advised to have cells of the aspect ratio close to unity, thus the grid resolution near the wall dictates its resolution in the rest of the domain.
2.6.1
Filtering
The first step in LES is the filtering in order to separate the quantities which are to be computed explicitly from those that will be modelled. For any variable φ its filtered quantity φ (in analogy with RANS also denoted by the overline) is calculated from:
φ =
Z x+∆
x
G(x, x0)φ(x0)dx0 (2.26)
where ∆ is a filter width and G(x, x0) is a filter kernel which determines the size and the
structures of the small scales.
The filter used in this research is the box filter, which assumes that the value computed in the cell centre is also the filtered value in that cell. In the physical space the box filter at point x is defined as:
G(x, x0) =
½
1/∆ if |x − x0| < ∆/2
0 otherwise (2.27)
where x0 is the point around x over which the box filter is being defined, and the filter
width ∆ is the local grid spacing. For an unstructured code that handles arbitrary cell
shapes the only appropriate definition is ∆ = (V)13, with V being the cell volume.
2.6.2
SGS modelling
Most of the presently used SGS models are based on the eddy viscosity concept of Boussi-nesq (1877), and are of the form:
τijSGS −1
3τkkδij = −2νtSij, (2.28)
where 1
3τkkδij is the kinetic energy of the unresolved velocity field, νt is the turbulent
viscosity and Sij = 12(∂uj/∂xi+ ∂ui/∂xj) is the mean strain rate tensor.
Usually the eddy viscosity νt is obtained algebraically in order to avoid the
2.7. Validation test cases 17 algebraic models should describe their physics accurately. The most widely used SGS model is the one of Smagorinsky (1963):
νt = (CS0∆)2|S| (2.29)
where |S| = (2SijSij)1/2 and CS0 is the Smagorinsky constant. Unfortunately, different
flows seem to require different values of CS0.
Deardoff (1970) proposed CS0 = 0.1, while some other authors suggest different values:
CS0 = 0.17 by Lilly (1992), or for the channel flow CS0 = 0.065 by Ferziger (1995).
Furthermore, in order to account for viscous damping near the wall the Smagorinsky constant should be damped, and the most commonly used in LES calculations is the Van Driest damping (also used in RANS):
CS = CS0
³
1 − e−y+/A+´
2
(2.30)
where A+ = 25 is the damping constant, and y+ is the normalised distance from the
wall.
To avoid the arbitrariness in the definition of the Smagorinsky constant, Germano et al. (1991) proposed it to be calculated in a dynamic manner. The idea behind the dy-namic procedure is that at the cutoff wave number the resolved and the unresolved scales are similar, so that the stresses from the original fine grid and a somewhat coarser test grid (whose values are denoted by the hat) can be compared. Applying the Boussinesq definition both for the subgrid stresses with the filter width ∆, and the subtest stresses
with the filter width b∆ (which is usually taken to be b∆ = 2∆), after some transformations
the following CS definition can be obtained:
CS = − 1 2 LijMij MijMij (2.31)
where Lij = duiuj− buiubj represents the resolved turbulent stress of the scales between ∆
and b∆, and Mij = b∆2|bS|cSij− d∆2|S|Sij represents the contribution of the modelled stress
of those scales.
The dynamic adjustment of the Smagorinsky constant reproduces the correct
near-wall behaviour, and can account for backscattering by taking a negative value for CS.
LES is a very powerful tool for the analysis of fluid flow and heat transfer as it re-solves a wide spectrum of scales (only the smallest scales remain unresolved). As a result of the constant increase in computer power, this method is becoming an important tool for engineering applications, although still limited due to its computational requirements. Chapter 6 presents different considerations required by LES, together with some results for the case of fluid flow and heat transfer in complex geometries.
2.7
Validation test cases
it is to calculate the simplest possible generic flow with one dominating flow effect. But very often it is difficult to find such a simple test case, because there are no sufficient reliable data. And even if there are, sometimes it is difficult to get that data. The author remembers vividly the wet basement without light, where one of the reference articles
was stored. Getting that article was a bit odd experience,2 but the final outcome was
worth this small sacrifice.
This section presents the test cases and their flow characteristics used to validate the new model developments. The following cases are well documented in the literature, and are widely used as reference test cases. The results provided for each test case are to demonstrate the correct implementation of the numerical method and turbulence models in T-FlowS (originally developed by Niˇceno (2001) under the name T-Rex), the code that was used for all computations presented in this thesis. The numerical set-up for each of these test cases is explained in Appendix B, and more information about numerical details of T-FlowS is given in Appendix C.
2.7.1
Channel and pipe flows
The simplest wall-bounded non-homogeneous turbulent flows are fully developed steady circular pipe (Fig. 2.3) and plane channel flow (Fig. 2.2). This is why they are the basic test cases for turbulence modelling. Furthermore, the pipe flow has very high practical importance, as it is widely used for numerous engineering applications. In these simple cases the turbulence non-homogeneity is in the wall-normal direction only (or for the pipe in the radial direction), whereas the gradients and temporal variation of the statistical quantities are equal to zero in the streamwise and the spanwise direction (that is the axial and the azimuthal direction for the pipe). The force that drives the flow is the streamwise (axial) constant non-zero mean pressure gradient.
x z y p =
{
f(t) const. ∆ region region buffer viscous sublayer 2h logarithmic 1 10 100 1000 y+ 0 5 10 15 20 25 U+, DNS, Reτ=590 U+, v2-f, Reτ=590 U+, k-ε, Reτ=590Figure 2.2: Sketch of a plane channel (left) and the profiles of the velocity obtained with
the standard k-ε model and Durbin’s υ2-f model (right).
2
2.7. Validation test cases 19 If the mean pressure gradient is varying periodically, unsteady effects are observed. By this periodical variation, which can be either pulsating (if the temporal mean is non-zero) or oscillating (if the temporal mean equals non-zero), the shear generated at the wall is changed. The pulsating flow between two parallel planes or in a pipe is the simplest case for testing the unsteady behaviour of turbulence model.
Steady channel flow
The mean flow equation in the case of a steady flow between two flat planes shows that the total shear stress (the sum of the viscous and the turbulent shear stress) varies linearly with the distance to the wall (Tennekes and Lumley 1976; Durbin and Pettersson Reif 2000): dU+ dy+ − uv + = 1 − y + Reτ (2.32)
where y+ and U+ are the normalised distance to the wall and velocity respectively.
From Eq. (2.12) it follows that in the immediate vicinity of the wall the turbulent stress
is negligible (uv+ ∼ y+3), and in Eq. (2.32) only the viscous stress remains (dU+/dy+ ≈
1), which yields the linear velocity profile U+ = y+. The region in which the viscous
stress is dominant is usually referred to as the viscous sublayer, and the DNS results are
showing that it extends up to y+ ≈ 5.
For high Re far away from the wall (y+>≈ 40) the viscous term in Eq. (2.32) becomes
negligible, thus uv+≈ 1. Using the mixing length assumption ν
t = κuτy (Prandtl 1925)
gives the relation −uv = u2
τ = νtdU/dy = κuτydU/dy, which integrates to the log-law:
U+ = 1 κln ¡ y+¢+ B = 1 κln ¡ EUy+ ¢ (2.33) where the constant κ = 0.41 is named after Von K´arm´an who was the first to derive the logarithmic velocity profile from similarity arguments (Von K´arm´an 1930), and B = 5.1 is the integration constant obtained by equating the log-law and the linear profile for the given effective height of the viscous sublayer. For obvious reason, this region is referred to as the logarithmic region.
The part of the flow between the viscous and the logarithmic region (usually taken as
5 < y+ < 30) is called the buffer region, and there both stresses are of equal importance.
It is in this region that the linear and the logarithmic velocity profile are intersecting, hence the maximum of P and k occurs, and the gradient of U is still high. This is why for a good prediction of wall-bounded flow it is crucial to accurately resolve the buffer region.
A lot of reference data for fully developed turbulent channel flow can be found in the literature, but here DNS of Tanahashi et al. (2004) and Moser et al. (1999) were used
for channel flow without heat transfer (Reτ = 800 and Reτ = 590 respectively), Kasagi
et al. (1989) for Reτ = 150 channel flow with heat transfer, and experiments of Wei and
