LES, RANS and combined simulation of impinging flows and heat transfer

(1)

LES, RANS and Combined Simulation

of

(2)

(3)

LES, RANS and Combined Simulation of

Impinging Flows and Heat Transfer

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 9 januari 2006 om 10:30 uur

door

Muhamed HADˇZIABDI ´C

(4)

Prof. dr. K. Hanjali´c

Samestelling promotiecommissie:

Rector Magnificus voorzitter

Prof. dr. K. Hanjali´c Technische Universiteit Delft, promotor

Prof. dr. M.A. Leschziner Imperial College, London, UK

Prof. dr. ir. A.A. van Steenhoven Technische Universiteit Eindhoven

Prof. dr. ir. Th. van der Meer Twente Universiteit

Prof. dr. R.V.A. Oliemans Technische Universiteit Delft

Prof. dr. D.J.E.M. Roekaerts Technische Universiteit Delft

dr. H.J.J. Jonker Technische Universiteit Delft

(5)

(6)

(7)

Nomenclature

Roman

a heat diffusivity

CF L Courant-Friedrich-Lewis number u∆ t

d Cf skin friction coefficient

C convection term

Cs Smagorinsky coefficient

Cdes constant in the Detached Eddy Simulation model

Cα constant in the hybrid LES/RANS ζ − f model

cp heat capacity of a substance at constant pressure D pipe diameter or diameter of jet nozzle

D diffusion term

h channel half height

ht heat transfer coefficient

h enthalpy

f frequency

fb body force per unit volume

f elliptic relaxation parameter

L characteristic geometry of a flow

ls Smagorinsky lenght scale

l lenght scale proportional to the filter width ∆

k turbulence energy

k thermal conductivity

N u Nusselt number ht D

k

n unit vector orthogonal to S and directed outwards

p pressure

Pk production of turbulent kinetic energy

Pr Prandtl number

ν a

Prt turbulent Prandtl number

νt a

S surface enclosing a control volume (CV)

Sij strain rate tensor

t time

T characteristic time period

(12)

T temperature

Tτ friction temperature

qw ρ cp Uτ

t fluctuating component of the temperature T

U velocity of a flow Uτ friction velocity s τw ρ ui velocity in i direction

uiuj Reynolds stress tensor

Ub bulk velocity for channel flow

v fluid velocity (vector)

υ2 _{velocity scale}

y+ _{dimensionless distance} x, y, z Cartesian coordinates

r, θ, z Cylindrical-polar coordinates

Wb bulk velocity for pipe flow and impinging jet

qh heat flux vector

Q second invariant of the velocity-gradient tensor ₋1

2(Sij· Sij− Ωij· Ωij)

Greek

α length scale parameter max(1, L

∆)

L characteristic length scale of the flow

T characteristic time scale of the flow

γ blending factor

δ distance from the wall to the grid cell centre

∆ characteristic mesh size, filter width for top-hat filter (∆ V )13

ε dissipation

ζ flow variable υ

2 k

τ shear stress

τw wall shear stress

η Kolmogorov length scale

λ thermal conductivity

µ molecular viscosity

ν kinematic viscosity

νt eddy viscosity

νt,RAN S eddy viscosity obtained by RANS model νt,LES eddy viscosity obtained by LES model t uj turbulent heat-transfer vector

ρ fluid density

φ flow variable

< φ > mean component of the flow variable φ

(13)

CONTENTS vii ˆ

φ test-filtered flow variable φ

ψ conserved quantity per unit mass

∆x cell size in the x-direction

∆y cell size in the y-direction

∆z cell size in the z-direction

∆r cell size in the radial direction

r∆θ cell size in the azimuthal direction

∆ V control volume

Ωij vorticity tensor

Subscripts

+ indicates that the related variable is normalized by friction velocity and kinematic viscosity

Abbreviations

DES Detached Eddy Simulation

DNS Direct Numerical Simulation

LNS Limited Numerical Simulation

LES Large Eddy Simulation

LDA Laser Doppler Anemometry

PITM Partially Integrated Transport Model

PIV Particle Image Velocimetry

RANS Reynolds-Averaged Navier-Stokes

rms root mean square

SGS Subgrid-Scale Stresses

(14)

(15)

Chapter 1 Introduction

1.1 Motivation

Computational Fluid Dynamics (CFD) has become a standard tool for both fundamental research and for the industrial and environmental computations of turbulent flows and associated transport phenomena. For research purposes - understanding and discovering physics of turbulence mechanisms and interactions, structural and other details - the main methods are direct numerical simulations - DNS, and, to a lesser extent, large eddy simulation - LES. For practical solution of complex realistic problems, the dominant method is the Reynolds-averaged Navier-Stokes - the RANS approach, using statistical turbulence models, though, more recently - LES has also been gaining in popularity, at least for some flow types.

RANS models have dominated industrial CFD for over three decades. However, RANS models can hardly provide any information about the turbulence spectral dynamics and, thus, they cannot be regarded as reliable for investigating the physics of turbulence and turbulent transport processes. Moreover, the intuitive and empirical nature of the RANS models, which are used for closing the Reynolds equations, imply various approximations, which limit their role as a predictive tool. For these purposes, the only reliable compu-tational methods are the properly resolved direct numerical simulations (DNS), which provide full resolution in time and three-dimensional space of all eddy sizes including the smallest dissipative ones. However, the DNS is very expensive and limited to low Reynolds numbers and simple geometries.

The developments in LES and the continuous increase in computer power have bred expectation that LES will soon emerge as the industrial standard and, perhaps, eliminate RANS. Even though LES proved to be a powerful method, the very demanding grid-resolution requirements, especially in numerical simulation of near-wall turbulence at large Reynolds numbers, still limits LES to low-to-moderate Reynolds numbers and to relatively simple geometries. The perspectives outlined by Spalart et al. (1997), Pope (2000) and Hanjali´c (2005) and others clearly indicate that the situation will not improve in the foreseeable future.

It is known, however, that not all flow types pose the same problems to the LES. Separated flows with vortex shedding, those with strong large-scale coherent structures, such as in thermal convection over horizontal hot surfaces, and other flows with strong forcing, are usually less demanding in terms of grid resolution than the attached flows with

(16)

a broad turbulence spectrum. Especially demanding are non-equilibrium attached wall flows and in particular when wall phenomena, friction and heat transfer are in the focus, simply because there are no large eddies in the near-wall regions in such flows. Hence, despite impressive developments, there are still many open issues in LES of wall-bounded non-equilibrium attached or weakly separated flows. In most cases of practical relevance, designing a rational grid that will provide the desirable resolution and yet be affordable with the available computers, is the main challenge, especially for unknown complex flows. Unstructured grids are regarded to be more suitable than the structured ones, allowing better grid rationing, but complex cell shapes introduce additional problems with filtering and discretization, thus bringing additional uncertainties and diminishing the LES’s predictive power. All these uncertainties call for additional research and testing of LES in carefully selected demanding test flows, where the LES’s shortcomings can best be demonstrated and researched.

These issues provided the motivation for the first task elaborated in this thesis, which focuses on detailed LES studies of a round jet impinging normally on a flat surface of different temperature, thus also involving heat transfer. This geometry is selected because a round impinging jet has several interesting physical features that pose a challenge to LES. Of course, other generic test flows are also considered, such as a plane channel and a pipe flow in which the subgrid-scale models, numerical discretization and computations, as well as boundary and inflow conditions, have been tested.

Recognizing the limitations of RANS and LES, and in search for more efficient solution methods, the CFD community has recently turned its attention to hybrid LES/RANS modelling as an alternative strategy for the wall-bounded flows at high Reynolds numbers. The rationale behind the LES/RANS hybrid approach is to improve the accuracy of simulations while keeping the mesh size typical of that used in the off-wall LES or in the RANS computation. Whilst RANS also requires a fine-grid clustering in the near-wall layer, it allows much coarser spacing in spanwise and streamwise directions than is permitted in conventional LES. By applying a RANS model in the near-wall region, the mesh size can be kept in proportion to ln(Re)0.5_{, instead of Re}9/4_{, as required by LES} (Chapman 1979), whereas in the outer region the mesh can satisfy the away-from-wall LES constraint proportional to Re0.5_{. This makes the mesh criteria much less stiff than} in the classical LES approach.

While the concept of the hybrid approach sounds appealing, many issues need to be resolved before such methods can be used for reliable computations. Despite the fact that the governing equations in RANS and LES have the same form, it is noted that the RANS equations are obtained after applying Reynolds time- or ensemble averaging, whereas the LES equations are obtained after performing a filtering operation, giving thus different physical meaning to the same variables in the LES and RANS frameworks. Mixing the two approaches to solve a flow field requires some modifications of one or each of the methods and a careful matching at a prescribed or fictitious interface plane or region. These are only some of many issues that need to be researched.

(17)

1.1. MOTIVATION 3 RANS-type computational grid. The target test case was again the round impinging jet for which comprehensive data for validation have been provided by the fine-resolved LES. Alongside the LES and hybrid RANS/LES studies, some RANS models have also been studied and used for solving the same and some other test flows. One of the motivations was to illustrate the potential and advantages, but also the shortcomings of some popular RANS approaches as compared with LES and hybrid simulations. The other subtask in this research segment was the test of the performance of selected RANS models, which have been regarded as potentially suitable for use in combination with LES in a hybrid mode. Because in such combinations a RANS model is usually applied only near solid walls with LES prevailing in the outer region, the selected models here tested are those which have been developed aiming at capturing primarily the near-wall phenomena, thus not necessarily being very sophisticated.

The final task is more of a physical nature. Because very detailed information have been gathered by well-resolved and verified LES of a round impinging jet, these data have been used to study various physical features and phenomena encountered in impinging flows and associated heat transfer. It is recalled that next to their industrial relevance in enhancing heating and cooling in various applications, impinging flows contain several generic features of interest for studying physics of turbulence and heat transfer, such as stagnation, jet deflection (associated with strong acceleration), formation of wall jets and their radial spreading (in round jets) subjected to a strong adverse pressure gradient. Understanding in which way the turbulent structures and heat transfer in an impinging jet interact has been regarded as a prerequisite for effective control and enhancement of heat-transfer efficiency. Because in an impinging jet flow one can distinguish several different phenomena and flow regimes, this flow has been recognized and frequently served as a generic test case for various turbulence models and other modelling concepts.

(18)

1.2 Objectives

Based on the above motivations and arguments, the following thesis objectives have been formulated:

1. To perform LES of a round jet impinging on a flat surface of different temperatures using a finite-volume in-house computational code with unstructured computational grid.

Because the physical picture of an impinging jet flow structure and the heat-transfer mechanism is not unique, but depends on the jet configuration, i.e. orifice shape and profile, inflow turbulence, distance from the target surface, Reynolds number, the jet configuration was to be selected as to agree as close as possible with the experimental set-up with which the data will be compared. The two parameters, the jet-to-plate distance and the Reynolds number have been found to be espe-cially influential. Different combination of these parameters lead to different flow structures and heat-transfer distributions.

The configuration selected was a jet issuing from a fully-developed pipe flow at a moderate Reynolds number (≈ 20.000), with a jet-exit-to-plate distance H/D = 2, which corresponds to the configuration investigated experimentally by Baughn and Shimizu (1989) and Geers et al. (2005), and the LES performance was to be validated using these and other experimental data available in the literature and in our section at TU Delft.

(19)

1.3. OUTLINE 5 5. Using the verified and validated well-resolved LES,

• to extract information about the turbulence eddy structure, its generation prior to impingement, its subsequent evolution during radial flow spreading, in particular in the near-wall region, and the correlation of the eddy structure with the instantaneous wall temperature (thermal imprints);

• to investigate and explain the causes of mean heat transfer non-uniformities, especially the origin of the double-peak in the Nusselt number;

• to analyze the turbulent-kinetic-energy and turbulent-stresses budgets in the stagnation region, to verify (or disprove) some controversial experimental find-ings such a negative kinetic energy production in the stagnation zone, local instantaneous wall-jet separation, and others.

6. To draw conclusions arising from the research, especially on all interesting and novel findings.

1.3 Outline

This thesis consists of seven chapters.

In Chapter 2 an outline of the numerical methods used in the present research is presented. Section 2.2 presents the form of the governing equations for heat and fluid flow and their integration over control cells as used in the adopted computational method. A brief description of the discretization practice is given in Section 2.3. The chapter ends with some details about the implementation of the advanced convective schemes into the unstructured code.

Chapter 3 introduces the principle ideas behind Large Eddy Simulation. Section 3.2 describes two sub-grid-scale models, used in this research, the Smagorinsky and dynamic model. Some specifics of LES on unstructured grids are discussed in Section 3.3. Sec-tions 3.4 and 3.5 present the results of the two test cases, the plane channel and the pipe flow respectively. The flow structures are examined and discussed for the plane channel flow.

Chapter 4 begins with a description of the Reynolds-averaged-Navier Stokes (RANS) approach. The eddy-viscosity hypothesis is presented in Section 4.3. Within this section, three turbulence models, used in this research, are described and commented upon. As an example of the models performances, plane-channel-flow results are given in Section 4.3.4. Chapter 5 begins with the motivation behind the hybrid LES/RANS modelling. It continues with the literature overview. Section 5.3 presents the Detached Eddy Simula-tion. Plane channel-flow results are presented and commented upon. The zonal hybrid approach is discussed in Section 5.4. Sections 5.5 and 5.6 present two seamless hybrid models. The results of the plane channel flow are examined in detail. The consequences of the hybrid modelling on the flow structures are described. The chapter ends with conclusions on the feasibility of the hybrid modelling and different hybrid approaches.

(20)

review of major physical features and of challenges for simulation and modelling in this flow, is followed by Section 6.2 in which the computations with several engineering-type RANS models are presented. The models selected have been assessed primarily in re-gard to their potential as the near-wall models for hybridization with LES simulations. Section 6.3 presents the LES simulations. Effects of grid density, boundary conditions, discretization scheme and subgrid-scale models are discussed in detail. This is followed by illustration of the instantaneous fields, flow and eddy structure and their imprints on the wall-deducted by different structure-identification method, which lead to some interesting conclusions on the flow physics. Statistically averaged results for the first and second-moments of the velocity field, their budgets and integral properties of temperature and velocity fields (Nusselt number) are also discussed in comparison with experimen-tal data of Geers et al. (2005) and Cooper et al. (1993). The main issues regarding the heat-transfer computations are presented in Section 6.3.7. Section 6.4 presents the computations of the same flow using a hybrid RANS/LES approach. This section leans much on the previous one and provides comparison of some major data with the LES simulations.

(21)

Chapter 2 Numerical method

This chapter gives an overview of the numerical methods used in the present research. It provides a brief overview of the finite-volume method applied on unstructured grids. The implementation of some advanced convective schemes in the unstructured code is described in some details.

2.1 Introduction

The choice of a numerical method for unsteady simulation of turbulent flow is a very important issue. The chosen numerical method has to be accurate, but at the same time, it has to be robust and applicable to complex geometries. These two requirements make a choice of the numerical method a difficult task. In LES, the most widely used methods are various spectral approaches and the orthogonal finite-volume (or finite-difference) methods of fourth or higher order of accuracy. Unfortunately most of the high-order accuracy methods are restricted to simple geometries such as a channel. The target of the present research were high-Reynolds number flows in complex geometries, and the choice was the second-order-accurate finite-volume method for unstructured and hybrid grids (Demirdˇzić et al. 1997) with collocated grid arrangement. This type of solver is not widely used in LES simulations because of several reasons. It has certain limitations on the order of accuracy and generally demands significantly more allocated memory. The data structure is more complicated than in structured code and parallelization is conse-quently often a difficult task. But, on the other hand, for high-Reynolds-number flows in complex geometries the unstructured code is sometimes the only possible solution. Unstructured meshes offer the necessary flexibility in economizing meshing of the com-putational domain, which generally decreases the computing costs. For this research an in-house fluid-flow solver called T-FlowS was used. The code was developed by Dr. Bojan Niˇceno and has been continuously improved by this author and a number of users. The code was tested in a number of laminar and turbulent test cases for which it produced satisfactory results (Niˇceno and Hanjalić 1999; Niˇceno and Hanjalić 2001; Niˇceno 2001; Niˇceno et al. 2002; Niˇceno et al. 2005).

(22)

x y z P s V S r

Figure 2.1: General, pebble-shaped control volume (from Niˇceno (2001)).

2.2 Governing equations for heat and fluid flow

The governing equations which describe a turbulent flow and heat transfer are derived from the conservation laws of mass, momentum and energy. They form a set of par-tial differenpar-tial equations which can appear in different forms. The choice of the form of the governing equations is determined by the discretization and numerical methods adopted. The subsequent discretization of the governing equations and coding of com-puter programs have to be as straightforward as possible. The finite-volume method is based on an a priory integration of the equations over the grid cells, and solves the equations in the integral form. It is valid for an arbitrary control volume V bounded by surface S(Fig. 2.1). Throughout this Thesis, the fluid is assumed to be incompressible, Newtonian with constant physical properties. The conservation of mass is expressed by

the following equation: _Z

Sρu · dSn = 0, (2.1)

where ρ is the fluid density, u is the velocity vector (u = (u, v, w)T_{), and Sn is the} outward pointing surface vector.

The law of momentum conservation applied to a control volume results in the following equation, also know as Cauchy’s law of motion:

d dt Z V ρudV + Z Sρuu · dS = Z ST · ds + Z V ∇pdV + Z V fbdV, (2.2)

where T is the Cauchy stress tensor, ∇p is the pressure gradient and fb is the body force per unit volume. The stress tensor for an incompressible Newtonian fluid is expressed as:

T = 2µD, (2.3)

where D is the rate of deformation:

(23)

2.3. DISCRETIZATION 9 and: ∇uT ₌         ∂u ∂x ∂v ∂x ∂w ∂x ∂u ∂y ∂v ∂y ∂w ∂y ∂u ∂z ∂v ∂z ∂w ∂z         . (2.6)

The conservation of the enthalpy for a control volume reads: d dt Z V ρhdV + Z Sρhu · dS = Z Sq h· ds + Z Vs hdV, (2.7)

where h is the specific enthalpy, qh is the heat flux vector and sh is the heat source or sink. In Eq. (2.7) it is assumed that the fluid is thermally perfect, i.e. that enthalpy does not depend on pressure but only on the temperature:

h = cpT, (2.8)

where cp is the specific heat, and T is the temperature. The heat flux qh is related to the temperature gradient by the Fourier law:

qh = k∇T, (2.9)

where k is the thermal conductivity.

Equations (2.3) and (2.9) represent the response of the particular fluid to external effects (deformation of fluid or heat flux), and when they are introduced into the equations which govern conservation of momentum and enthalpy, the following forms are obtained, respectively: d dt Z V ρudV + Z Sρuu · dS = Z Sµ(∇u + (∇u) T ) · ds + Z V ∇pdV + Z V f bdV, (2.10) d dt Z V ρcpT dV + Z SρcpT u · dS = Z Sk∇T · ds + Z Vs hdV, (2.11)

which subsequently will be discretized by the finite-volume method.

2.3 Discretization

A conservation law for a generic variable can be written as follow: d dt Z V ρBφ ∂φ ∂tdV | {z } Rate of Change + Z SρBφφu · dS | {z } Convection = Z SΓφ∇φ · dS | {z } Diffusion + Z SqφSdS + Z V qφVdV | {z } Sources , (2.12)

(24)

general: φ Bφ Γφ qφS qφV

mass: 1 0 0 0 0

momentum: ui 1 µeff µeff(∇u)T fi+_∂x∂p_i

enthalpy: T cp keff 0 0

Table 2.1: Meaning of various terms in general transport equation 2.12.

s 2 s_i s₁ s 3 r z y x Ci C3 C1 C0 C2 0

Figure 2.2: Possible control volume (from Niˇceno (2001)).

2.3.1 Discretization in space

The finite-volume method on an unstructured grid (see Demirdˇzi´c et al. (1997), Niˇceno and Hanjali´c (2001)) is based on spatial discretization of the governing equations. The cells of any shape can be used. A possible control volume is shown in Fig. 2.2. The computational point C0 is located at the center of the control volume such that:

Z

V(r − r0)dV = 0. (2.13)

All dependent variables, velocity components and pressure, are stored in the cell centers, i.e. they are co-located. Staggered grid arrangements have been found to be more robust than the collocated ones, but their implementation on unstructured grids is very complex (Rosenfeld et al. 1991; Rosenfeld 1996; Wesseling et al. 1992).

All dependent variables are assumed to vary linearly over the cell, according to:

φ = φ0+ ∇φ · (r − r0). (2.14)

The governing equations consist of several volume and surface integrals. Taking into account Eq. (2.13) and Eq. (2.14) it follows that:

Z

V φ(x)dV ≈ φ0∆V, (2.15)

where ∆V is the volume of the computational cell. Hereafter, the subscript 0 is omitted from the cell centered values for the sake of simplicity. In a similar way, for the integration on the cell face it can be deduced that:

Z

(25)

2.3. DISCRETIZATION 11 where j denotes a cell-face value, defined at the cell face centroid, Sjnj is the discrete surface vector of the cell face. The form of equation (2.16) shows that a mid-point rule can be applied in discretization of all volume and surface integrals.

The dependent variables at the cell face are calculated by a linear interpolation: φj =

φ1,j + (∇φ)1,j · (rj− r1,j) + φ2,j+ (∇φ)2,j· (rj − r2,j)

2 , (2.17)

where subscripts 1, j and 2, j denote the cells surrounding the cell face j. Hereforth the subscripts 1, j and 2, j are replaced by 1 and 2 for the sake of simplicity. The gradients of dependent variables are calculated at cell centers, and are assumed to be constant in each cell. The values of the gradients at cell faces are calculated from:

(∇φ)j = (∇φ)1+ (∇φ)2

2 , (2.18)

where (∇φ)j stands for the gradient at the cell face and (∇φ)1 and (∇φ)2 are gradients at cells surrounding the face. The gradients at the cell centers are calculated by applying least-square fit.

∇(φ) ≈ G2(φ), (2.19)

where G2 is the second-order accurate discrete gradient operator. Convective terms

Special attention was paid to the treatment of the convective terms. These nonlinear terms are usually the main source of numerical instabilities when the Central Differenc-ing Scheme (CDS) is used for discretization. Central differencDifferenc-ing for convection causes instabilities when the local Peclet number reaches values of order of 103, no matter what the value of the local CFL number is (Niˇceno 2001). Instabilities do not necessarily lead to a divergence of the solution, but can be in the form of nonphysical oscillations. These kinds of instabilities are especially dangerous because they are difficult to detect. In order to prevent the occurrence of nonphysical oscillations several numerical remedies can be used, such as the mesh refinement (very undesirable because of the increase in the cost of computation) or the use of more diffusive convection schemes, usually high-order upwind-biased schemes (not recommended for LES because of dissipative nature of these schemes) or to blend CDS scheme with some percentage (usually very small one) of high-order upwind-biased schemes.

The code uses the mid-point rule for discretization of the convective terms, i.e. by multiplying the cell face value with the surface of the cell face:

Z

Sρφu · Sdn ≈ C2(φ) =

X

j

φjFj, (2.20)

where C2 is the discrete, 2nd order accurate convection operator discretised by central differencing scheme, and Fj are the mass fluxes defined on the cell faces by:

Fj = [(ρu + ∆tG2(p))j+ ∆t

(p1 − p2) |dj|

(26)

where |dj| is the distance between two neighboring cell centers. Eq. (2.21), represents the implementation of Rhie and Chow (1983) technique on arbitrary grids. The cell-centered pressure gradients (G2(p)) are extracted from cell-centered momentum equations, and the staggered terms (p1− p2)/|dj| are added.

An additional term C1 is introduced if the blending of CDS scheme with high-order upwind-biased schemes is to be conducted. C1 is discretized by one of the upwind-biased high-order schemes(e.g. QUICK, LUDS, SMART).

The complete convective term reads:

C(φ) = CL(φ) + CR(φ), (2.22)

where CL is the implicit and CR is the explicit 1 part of the convective flux.

The implicit part (L) is added to the system matrix to increase the stability of the linear system, whereas the explicit part (R) is added to the right hand side. The implicit part is defined as follows:

CL(φ) = C1(φ), (2.23)

and explicit as:

CR(φ) = γ(C2(φ) − C1(φ)), (2.24)

where γ is the blending factor ranging from 0 to 1. When γ is 1 the full central differencing scheme is obtained, which is second-order accurate and most often used in LES. In this work the local CF L number was defined on the cell faces as:

CF L = uj· dj∆t |dj|

(2.25) Diffusive terms

The diffusive fluxes are discretised by approximating the integrals with sums, yielding the following expression:

Z

SΓφ∇φ · Sdn ≈

X

j

ΓφG2(φ)jSjnj, (2.26)

where G2(φ)_j is the value of the gradient defined on the cell face. The straightforward dis-cretization of Eq. (2.26) leads to a numerical stencil containing first and second neighbors of each cell. In order to make the parallelization by a domain decomposition technique easier and more simple the gradient at the cell face is decomposed according to the following expression: (∇φ)?j = (∇φ)j+ φ2− φ1 |dj| nj+ (∇φ)j· dj |dj| nj. (2.27)

The first term on the right-hand side is the second-order accurate gradient of φ, the second term on the right hand side is the first order accurate gradient of φ which is calculated from the values surrounding the cell face j. The third term on the right hand side

1_{In principle, terms implicit and explicit as introduced here do not have any connection with time}

(27)

2.3. DISCRETIZATION 13

(a) (b)

Figure 2.3: (a) Numerical stencil for discretization of explicit viscous stresses DE on unstructured grid and (b) necessary communication pattern after parallelization (from Niˇceno (2001)).

was introduced to recover second-order accuracy of the discretization. This additional term represents the difference between the second- and first-order approximation of the derivative at the cell face. It contains values from the wider numerical stencil (Fig. 2.3(a)), but it is treated explicitly using the values either from the previous iteration, or old time step, leaving the number of entries in the system matrix low.

The final form of the diffusive terms is: DL(φ) = X j Γφ[ φ2− φ1 |dj| nj]Sjnj, (2.28) DR(φ) = X j Γφ[(∇φ)j +(∇φ)j · dj |dj| nj]Sjnj. (2.29)

where DL is the implicit and DR is the explicit part of the diffusive flux. Surface sources

The surface sources are discretised using the mid-point rule:

Z

SqφSdS ≈ QS(φ) = qφSj∆Sjnj. (2.30)

Volume sources

The midpoint rule is used to discretize the volume sources, i.e. :

Z

(28)

(a) (b)

Figure 2.4: (a) Numerical stencil for discretization of implicit viscous stresses DI on unstructured grid and (b) necessary communication pattern after parallelization (from Niˇceno (2001)).

The final form of the governing equations

The final form of the discretised system of equations is written as follows:

C(1) = 0, (2.32)

which is the discretised continuity equation and:

I(u) + C(u) = D(u) + QS(u) + QV(u), (2.33)

which is the discretised Navier-Stokes equation, and

I(T ) + C(T ) = D(T ). (2.34)

which is the discretised enthalpy equation.

It has to be noted that convective and diffusive terms were decomposed into implicit and explicit part:

C(u) = CL(u) + CR(u), (2.35)

D(u) = DL(u) + DR(u), (2.36)

and, depending on the time stepping scheme, each part can be treated differently in time.

2.3.2 Integration in time

Unsteady terms are integrated by the midpoint rule, i.e. the value of the dependent variable in the cell center is assumed to be representative for the entire cell. A linear variation between the current and previous time step lead to the following form of the discretised inertial terms:

(29)

2.4. DATA STRUCTURE 15 C₁ C₂ S C₁ S C₂ C 1 C₂ S C₁ C2 S

Figure 2.5: Cell face data structure and some of the possible shapes of control volumes (from Niˇceno (2001)).

where I is the discrete inertial operator, n is the new time, n − 1 is the old time and ∆t is the discrete time step. Three different time advancing schemes have been implemented for time integration: the semi-implicit fractional-step method (Chorin 1968; Gresho 1990; Kim and Moin 1985), the semi-implicit Crank-Nicolson and the fully-implicit Gear scheme (Fletcher 1991). The fractional-step method uses the second-order Adams-Bashforth method for the convective terms and Crank-Nicolson method for diffusive terms. This scheme poses a too restrictive limit on the numerical time step and requires high-quality grids. The semi-implicit Crank-Nicolson scheme is more stable than the fractional-step method (Niˇceno 2001). By this scheme all terms are discretised with the Crank-Nicolson scheme. Finally, the fully-implicit Gear scheme discretises the inertial term using the quadratic span over three time steps. This scheme appeared to be the most economic in terms of computer memory needed. All three schemes are of second-order accuracy in time.

2.4 Data structure

(30)

2.5 Auxiliary numerical details

2.5.1 Implementation of bounded convective schemes for

un-structured grids

The implemented bounded convective schemes have been described in Przulj and Basara (2001). Transport of a scalar φ is considered along the local coordinate ψ, which is defined by the upstream (U), central (C) and downstream (D) computational nodes, Fig.2.6(a). The velocity direction determines the actual labelling of the nodes. The cell face value φ has to fulfill the following two conditions:

• it is bounded by the neighboring cell values φC and φD;

• it ensures monotonicity of the local function φ(ψ) which passes through the points φU and φC.

It is convenient to introduce the normalized variable ˜φ defined as: ˜

φ = φ − φU φD− φU

(2.38) giving ˜φD = 1 and ˜φU = 0.

Fig.2.6(b), where the case with φC > φU is illustrated, shows the meaning of these conditions. The first condition imposed interpolative boundedness and can be written as: 0 < φj− φC φD − φC < 1 or 0 < φ˜j − ˜φC 1 − ˜φC < 1 (2.39)

The second, monotonicity condition, reads: φj− φC φC− φU > 0 or φ˜j − ˜φC ˜ φC > 0 (2.40)

Both conditions are satisfied simultaneously if the following relations hold: ˜ φj < 1 and φ˜j > ˜φC if φ˜C < 1 (2.41) ˜ φj = ˜φC if φ˜C < 0 or φ˜C > 1 (2.42)

These conditions, called CBC constraints, are derived by Gaskell and Lau (1988). If the monotonic profile is considered, defined by φD − φC

φC − φU

> 1, one can bound the cell face value by the extrapolated downstream value φDe= 2φC− φU (Fig.2.6(b)) rather than by the actual value φD. This constraint can be expressed as:

φj− φC φC− φU < 1 or φ˜j − ˜φC ˜ φC < 1 (2.43)

(31)

2.5. AUXILIARY NUMERICAL DETAILS 17 P C U D Pj j n_j m. 0 ξ d_j (a)

φ

U

φ

D D D D C

φ

C D e e e D

U

C

j

D

φ

_j

ξ

(b)

Figure 2.6: (a): Upstream, central and downstream nodes; (b): An illustration of variable values at computational nodes.

˜

φj < 1 and ˜φj < 2 ˜φC and ˜φj > ˜φC if 0 < ˜φC < 1 ˜

φj = ˜φC if ˜φC < 0 or ˜φC > 1

(2.44) The only avaliable values of φ for arbitrary unstructured grids is in neighboring points C and D. A gradient of φ is used in point C to reconstruct the needed value of φU

∗ as: φU

∗

= φD− ∆φC · 2 ~dj (2.45)

The star (∗) denotes an imaginary upstream cell which is defined in such a way that the vector identity ~_{CU = − ~}CD is satisfied and that an imaginary face between C and U is placed at the same distance from C as the considered face j, Fig.2.6(a). The reconstructed value of φU

∗

has to be limited by the values of φ at neighboring cells which surround the upstream node. This is done by checking φU∗ against the minimum and maximum values of φ over the cells that share faces of the central cell.

φU = max[ φneighboursmin , min( φU,φneighbourmax ) ] (2.46) The TVD and CBC conditions are illustrated in Fig.2.7. The values are normalized as previously noted. The cell-face values ˜φj should lie within the shaded areas in the monotonic range 0 < ˜φC < 1, and on the line ˜φj = ˜φC outside the monotonic range. Some of the convective schemes with linear characteristic in the NVD (Normalized Variable Diagram) diagram may violate the boundedness criteria. It is obvious from the NVD diagram that CDS, LUDS and QUICK are that kind of schemes.

A general function ˜φj = f ( ˜φC) can be expressed in the form that includes the flux limiter:

˜

φj = ˜φC + ϕj( 1 − ˜φC) (2.47)

(32)

UDS LUDS CDS Q MINMOD C j C 2 f_j* 0.5 1.0 1.0 φ φ φ (a) C j 0.5 1.0 1.0 φ φ φC 2.25 φC 3.0 AVL smart quick smart Q o (b)

Figure 2.7: (a): TVD criteria and characteristics of several convective schemes; (b): CBC criteria and characteristics of several convective schemes

Eq. (2.47) can be expressed in terms of normalized variables: φj = φU DSC + ϕj

˙ mj | ˙mj|

( φP j − φP) (2.49)

The unbounded limiter for the situation shown in Fig.2.7 can be obtained as:

ϕj = (gD− αj) + (gU + αj)rj (2.50) gD = 1 2f ∗ j(1 + f ∗ j) , gU = 1 2f ∗ j(1 − f ∗ j) (2.51)

Families of the second- and third-order accurate schemes are defined by parameter α. The convective flux C Eq. 2.22 can be discretized as:

(33)

Chapter 3 Large-eddy Simulation

This chapter provides a short overview of basic principles of the large-eddy simulation (LES) technique, with specific reference to issues related to unstructured computational meshes. For testing the computational code in general, of the implementation and perfor-mance of the subgrid-scale models, and of the transformation of the code into cylindrical coordinates, LES results for a plane channel and a pipe flow at several Reynolds numbers in the range of Reτ = 180−590 are presented. The obtained results were used to examine the turbulence structures especially close to the wall in the plane channel flow in order to identify the streak structures. These structures were subsequently used as an indicator for validating the resolution performance of the hybrid LES/RANS models. Whereas only the hydrodynamic field was considered in channel flows, the pipe flow at Reτ = 180 was simulated together with the heat transfer. The budgets of the transport equations for the Reynolds stress components and the turbulent heat-flux-vector components are calculated for this case. The main purpose of this calculation was to investigate the performance of LES in cylindrical coordinates and the potential of LES to accurately calculate budgets of the transport equations for Reynolds stresses and components of the turbulent heat-transfer vector.

3.1 Introduction

The fluid motion is mathematically described by the Navier-Stokes equations of motion, complemented by the continuity and energy equations. Here, incompressible, constant-density fluids are considered for which these equations read:

(34)

where T is the temperature of the fluid and a is the heat diffusivity. Here, thermally perfect fluid is assumed, i.e. that enthalpy does not depend on temperature and pressure but on the temperature only.

The direct numerical solution of the Navier-Stokes equations is called Direct Numer-ical Simulation (DNS). If properly resolved, the DNS solution represents a real picture of turbulence containing all scales of motion, from the largest scales to the smallest dissipative scales. Since DNS resolves all sizes of the motion, the computational cost in-creases with the Reynolds number enormously, roughly with Re2.2_{. This limits the DNS} to low-Reynolds-number flows in relatively simple geometries. However, by performing a simulation which resolve only the large, energy-containing eddies and model the rest of the velocity spectrum, one can significantly reduce the cost of the simulation. This type of simulation is called Large Eddy Simulation (LES). LES is based on an assumption that large, energy-containing eddies determine the main flow dynamics. The small eddies that lie in the inertial and dissipative subranges of the energy spectra are assumed to be isotropic and not to influence significantly the main flow properties, but nevertheless they need to be modelled.

DNS, τmod = 0 LES, τmod = τLES RANS, τmod = τRAN S

E

log(κ)

Figure 3.1: Energy spectrum of isotropic turbulence (solid line) as a function of the wave number κ. The quantity τmoddesignates the part of the spectrum which is to be modelled by the three approaches, DNS, LES and RANS.

In order to separate the small scales of motion from the large scales, a filtering oper-ation is conducted over the Navier-Stokes equoper-ations. The velocity U (x, t) is decomposed into the sum of a filtered component U (x, t) and a subgrid scale component u(x, t):

U (x, t) = U (x, t) + u(x, t). (3.4)

The filtered velocity U (x, t) represents the motion of the large eddies. Similarly, the temperature is decomposed into the sum of a filtered and subgrid scale components:

(35)

3.1. INTRODUCTION 21 Notation of the filtered temperature T will be replaced with T for the sake of sim-plicity.

Leonard (1974) proposed to define U (x, t) by: U (x, t) =

Z +∞

−∞ G(r, x)U (x − r, t)dr. (3.6)

This type of an integral is called a convolution. G(r, x) is a filter function that satisfies the normalization condition:

Z +∞

−∞ G(r, x)dr = 1. (3.7)

The most commonly used filter functions are:

• The sharp Fourier cutoff filter (Leonard 1974) defined in wave space as: G(κ) =      1 if κ < π/2 0 otherwise (3.8)

• The Gaussian filter:

G(r) = s 6 π∆2exp − 6r2 ∆2 ! (3.9) • The box filter:

G(r) =      1/∆ _{if |r| < ∆/2} 0 otherwise (3.10)

Since a finite-volume unstructured code was used for this research, the box filter was used. With a finite-volume discretization, filter G does not appear explicitely at all. The way the implicit filter works is to discard any scale smaller than the grid size.

Applying the filter operator to the Navier-Stokes equations, the following equations for the filtered velocity components are obtained:

∂Ui ∂xi = 0, (3.11) ∂Ui ∂t + ∂UiUj ∂xj = ν ∂ 2_U i ∂xi∂xj − 1 ρ ∂p ∂xj , (3.12) ∂T ∂t + Uj ∂T ∂xj = ∂ ∂xj a∂T ∂xj − u jt ! . (3.13)

Eq. (3.12) differs from the Navier-Stokes equation (3.2) in the term UiUj. The dif-ference between the two terms represents the impact of the unresolved eddies on the resolved velocity field. It is defined by:

(36)

U ,U

u,u

Figure 3.2: Upper curve: a sample of the velocity field U (x) and the corresponding filtered field U (x) (bold line), using the Gaussian filter with ∆ ≈ 0.35. Lower curves: the residual field u(x) and the filtered residual field u(x) (bold line), (Pope (2000), reprinted with permission).

τR

ij is the residual-stress tensor which has to be modelled. It is analogous to the Reynolds-stress tensor:

huiuji = hUiUji − hUiihUji, (3.15)

where hi denotes Reynolds time- or ensemble- averaging.

The turbulent heat flux ujt, that appeared in Eq. 3.13, is modeled in a similar way: ujt = − νt Prt ∂T ∂xj , (3.16)

where Prt = 0.9 is the turbulent Prandtl number. The equation (3.16) is also known as simple gradient diffusion hypothesis (SGDH).

τR

ij comes from the nonlinear, convection term. In order to close the system of equa-tions, this term has to be modelled. It is obvious from the definition of the filtered velocity that the field of U (x, t) is random, three-dimensional and unsteady, see Fig. 3.2. It is noted that τR

(37)

3.2. SUBGRID-SCALE MODELLING 23

3.2 Subgrid-scale modelling

The effect of small structures, lost due to the filtering operation, has to be replaced by an appropriate subgrid-scale model. Compared to RANS modelling, where the whole spectrum is modelled, the subgrid scale modelling in LES does not influence the results on such a scale. As the cutoff wavelength is supposed to be in the inertial subrange of energy spectrum, the main burden lies in the resolved motion obtained directly from the resolved Navier-Stokes equations. The task of the SGS model is to dissipate the energy coming from the energy-containing range of the spectrum. Usually the dominant scales of the flow are much larger than the grid scales. This means that the coupling between these scales is weak. As the scales of motion become finer, they become more universal and lose their original orientation. This makes it easier to model them. In contrast, if the grid is coarse and the cut off wavenumber is close to the range of the energetic, anisotropic and inhomogeneous scales, the SGS model needs to be more sophisticated in order to correctly represent the influence of these scales.

3.2.1 Smagorinsky subgrid-scale model

The Smagorinsky model (Smagorinsky 1963) was the first SGS model proposed, and it is still widely used. It is based on the Boussinesq eddy-viscosity concept which relates the traceless part of the SGS stress components, τR

ij, to the strain of rate Sij of the resolved velocity field:

τ_ijR₋1

3τkkδij = −2νtSij. (3.17)

The eddy viscosity νt has to be determined. Dimensional analysis shows:

νt ∝ lsusgs, (3.18)

where ls is the length scale of the unresolved motion and usgs its velocity scale. There is still some discussion about what is the appropriate definition of the length scale of the unresolved motion. A natural choice is to determine the length scale by the mesh size. But if the grid is not uniform in all directions, it is not clear what is the appropriate length scale. The most widely used definition, which is also adopted here, is:

ls= ∆ = (∆V )

1

3, (3.19)

where ∆V is the volume of the computational cell.

The velocity scale is related to the gradients of Ui and it is defined as:

usgs = ls|S|, (3.20)

where |S| =q2Sij Sij is the magnitude of the strain rate Sij defined as: Sij = 0.5( ∂Ui ∂xj +∂Uj ∂xi ). (3.21)

(38)

νt = (Cs∆)2|S|. (3.22) Several drawbacks arise with the Smagorinsky model. The non-universality of the Smagorinsky constant Cs is one of them. Different types of flows have different optimal value of Cs. The values of Csvary from 0.065 (found to be optimal for the channel flow by Ferziger (1995)) to 0.2 for isotropic turbulence. Furthermore, Cshas to be reduced in the near wall region to account for the turbulence anisotropy. So-called damping functions are used in order to obtain the correct distribution of Cs in the near wall region. The most often used is the van Driest damping function:

Cs = Cs0

1 − e−_y+

/252_, _(3.23)

where y+ _{is the distance from the nearest wall in the viscous wall units y}+ ₌ yUτ ν (Uτ is the friction velocity defined as Uτ =

q

τ /ρ, where τ is the wall shear stress) and Cs0 is the Smagorinsky constant in the region away from solid boundaries. The damping functions are empirical expressions that add the additional uncertainty in flows with complex geometries. Another drawback is that the model is strictly dissipative (νt ≥ 0) and cannot predict the energy transfer in opposite direction, the so-called backscatter. However, the Smagorinsky model is still very popular, mainly because of its simplicity. The model is successful as long as the mesh is sufficiently fine and the transition processes are not dominant.

3.2.2 Dynamic subgrid-scale model

The dynamic model developed by Germano et al. (1991) overcomes the main drawbacks of the standard Smagorinsky model. The model uses the same Boussinesq definition of the unresolved stresses but the constant Cs is determined from a dynamic procedure. The dynamic procedure is based on the idea that the resolved scales close to the cutoff scale are similar to the modelled ones. The subgrid scale model is applied not only on the grid scale, defined by ∆, but also on the coarse scale defined by ∆. This is usuallyb

b

∆ = 2∆, and∆ is called the test filter. Now, stresses for the sub-grid and test levels canb be defined: τij = UiUj − Ui Uj ≈ τijmod C, ∆, U (3.24) Tij =UdiUj −Uci Ucj ≈ Tijmod C,∆,b Ub (3.25) The test-scale stresses Tijare obtained when the second, coarser filter is applied on the already filtered Navier-Stokes equations with filter width ∆, Eq. (3.12). From Eq. (3.25) the total stresses UiUj can be written as the sum of the resolved stresses on the grid ∆, Ui Uj, and the remainder τij. Inserting UiUj in the expression for Tij gives:

Tij = (Udi Uj−Uci Ucj)

| {z }

Lij

(39)

3.2. SUBGRID-SCALE MODELLING 25

Lij

resolved scales unresolved scales

π b ∆ π ∆ Tij τij logκ logE(κ)

Figure 3.3: Illustration of the sub-grid and test regions.

Expression (3.26) is known as the Germano identity. Lij represents the resolved stresses of the scales whose length is intermediate between the grid filter and the test filter widths. From the known resolved velocities Ui, the velocities Uci can be computed by applying the filter ∆ to Ub i associated with an appropriate function G. Similarly, Lb ij can be computed using already known velocity field Ui. On the other hand, the term Lij can be expressed by the modelled stresses τijmod and Tijmod using Eq. (3.26):

Lmod_ij _{= −2C}∆b2_|_S|b dSij | {z } Tmod ij +2C∆d2_|S|S ij | {z } −τmod ij . (3.27)

In an ideal case, Lij should be equal to Lmodij :

Lij− Lmodij = 0. (3.28)

From this condition Cs can be derived. As Eq. (3.28) is a tensor equation, it can be fulfilled only in some average sense. This is done by the least-squares minimalization proposed by Lilly (1992).

First, Eq. (3.27) can be written as:

Lmod_ij _{= −2C (}∆b2_|_S|b dSij−2C∆d2|S|Sij)

| {z }

Mij

. (3.29)

Multiplying Eq. (3.28) by Mij yields: Cs= − 1 2 Lij Mij Mij Mij . (3.30)

(40)

Figure 3.4: Illustration of the unstructured mesh used for the simulation of the round impinging jet.

feature of the model is the possibility of Cs to be negative, accounting for the effects of the backscatter. The price one has to pay for using the dynamic procedure to determine Cs is the increased complexity of the model and the higher cost of the simulation.

3.3 Unstructured-grid Large Eddy Simulation

The LES on unstructured grids is a relatively novel concept. The need for unstructured grids appeared at the moment when LES moved from the problems with relatively simple geometries and moderate Reynolds numbers to high-Reynolds numbers and complex-domain flows. It has been realized that the structured grids are less efficient in meshing complex domains then the unstructured ones. As the wall bounded flows require very fine mesh clustered in all three directions near the wall, the lack of flexibility to economize with the grids can lead to a prohibitively high cost of the simulation. Chapman (1979), and later Jim´enez and Moin (1991), argued that for specific complex, high-Reynolds-number flows, the structured grid should be abandoned to give way to the unstructured one. However, the unstructured-grid LES is not without problems. Some of the high-order convective schemes are not applicable for the unstructured codes which implies a higher numerical errors. Generally, implementation of turbulence models and different numerical schemes is more complex. The unstructured codes require more computer memory than the structured ones, which is unattractive. Specialized structured LES codes are usually faster then unstructured ones. Fig. 3.4 shows the unstructured mesh for the round impinging jet. It can be seen that the great number of cells are saved combining prismatic cells with hexahedral cells. Special attention has to be paid to the implementation of the dynamic procedure used in the Dynamic Smagorinsky model. The following filter operation is defined:

c φc_{= φ} 0fc+ nc X i=1 fiφci. (3.31) where φc

i is the value of the function at the center of the cell on the other side of the ith_{face, and f}

(41)

3.4. LES OF A PLANE CHANNEL FLOW 27

channel walls

2h

Flow direction

x, U

z, W

y, V

Figure 3.5: Schematic representation of the problem domain for the channel flow simu-lation.

3.4 LES of a plane channel flow

Fully-developed turbulent channel flow has long served as a standard benchmark for testing turbulence models both in the RANS and LES approaches. Detailed experimental research and well-resolved DNS have been performed by several authors. LES of a channel flow is performed in order to test the code and observe the channel-flow structures. Since a part of this research was the development of a hybrid RANS/LES turbulence model (for which the channel flow was the basic test case), it is considered that prior to these test a well-resolved LES should be conducted with the same computational code and numerical mesh.

3.4.1 Subgrid-scale model and computation details

A channel flow at Reτ = Uτh

(42)

differencing was used for both the diffusive and convective terms, without any under-relaxation procedure. The simulation was conducted on the supercomputer TERAS in Amsterdam. The domain was divided into 16 sub-domains. Statistics were collected over a period of 40 flow-through times.

3.4.2 Results

1 10 100 1000 y+ 0 5 10 15 20 25 30 U + (a) 0 100 200 300 400 500 600 0 2 4 6 8 10 uu + 0 100 200 300 400 500 600 0 0.5 1 1.5 2 vv + 0 100 200 300 400 500 600 y+ 0 0.5 1 1.5 2 ww + 0 100 200 300 400 500 600 y+ -0.9 -0.6 -0.3 0 uv + (b)

Figure 3.6: (a): Mean-velocity profile; (b): Resolved stresses; line: LES, symbols: DNS data of Kim, Moin and Moser (1989).

The LES results are compared with the DNS data of Kim et al. (1987) at Reτ = 590. The real Reτ, obtained by the LES simulation, was Reτ = 560. The mean-velocity profile, scaled with the inner-wall scales, Fig. 3.6(a), shows good agreement with the DNS data.

The error in the prediction of the skin-friction coefficient was defined as: Error = |Cf − Cf,Dean| Cf,Dean · 100%, (3.32) where Cf is calculated as Cf = τw 0.5 ρ U2 b , here τw = ν ∂U ∂y ! wall

is the wall shear stress and Ub is the bulk velocity.

The Cf,Dean was proposed by Dean (1978) and defined as: Cf,Dean = 0.073 Re

−1/4

b , (3.33)

where Reb = Ub h

ν , h is half of the channel height. Eq. (3.32) was also used for any error estimation in all other hybrid-scheme test cases.

(43)

3.4. LES OF A PLANE CHANNEL FLOW 29 play an important role in the production of turbulent kinetic energy. The agreement of the normalized shear stress with the DNS result is very good. It appeared that the contribution of the modelled energy compared to the resolved one is rather small. The modelled part of the total shear stress did not exceed 5 % of the resolved shear stress. The good agreement of the LES results compared with the DNS data allowed us to assume that the most important structures of the turbulent plane channel flow are truthfully resolved. In the next section these structures will be discussed in more detail.

3.4.3 Channel flow structures

Since the channel flow has been used as the basic test case for the hybrid LES/RANS models, which in principle are not expected to resolve streaks and other small-scale struc-tures in the near-wall region, the fine-resolved LES is used to identify the characteristic flow structures in order to establish the reference structure pattern, as well as to under-stand their role in the turbulence-production process. Channel-flow structures have been examined and reported in details, both experimentally and numerically, by a number of authors. Among others, experimental work has been carried out by Wei and Willmarth (1989), and numerical experiment has been conducted by Kim et al. (1987). A useful categorization of the flow structures has been made by Kline and Robinson (1990) and Robinson (1991). They identified eight categories of the coherent structures:

• Low-speed streaks in the region 0 < y+ _{< 10.} • Ejection of low-speed fluid outward from the wall. • Sweeps of high-speed fluid towards the wall. • Vortical structures of several proposed forms.

• Strong internal shear layers in the wall zone (y+_{< 80).}

• Near-wall pockets, small regions without marked fluid in certain types of flow vi-sualizations.

• Backs: surface across which the streamwise velocity changes abruptly. • Large-scale motions in the outer layers.

(44)

of turbulent kinetic energy Pk = −uv ∂U

∂y. It has been observed by different researchers, e.g. Wallace et al. (1972), that the process of ejection and sweeps is important for turbulence production. Fig. 3.8 shows instantaneous U velocity field at the plane parallel to the wall at y+ = 6. The contours of the streaks can be clearly seen. In the region y+_{< 100, the dominant structures are counter-rotating streamwise vortices or rolls. The} fluid, moving between the rolls towards the outer region, reduces axial velocity, sketched in Fig. 3.7(b). This can lead to the inflection point in the velocity profile, which means the flow is inviscidly unstable. It is believed that this is connected to the bursting process (Holmes et al. 1996). Fig. 3.9a shows a vector plot in the plane perpendicular to the flow direction. The near wall detail from the top wall, zoomed at Fig. 3.9(b), reveals a pair of the counter-rotating vortices. Further from the wall dominant structures are horseshoe or hairpin vortices. u 1 2 3 4 Ejections v Sweeps (a) x y z (b)

Figure 3.7: (a): The u−v sample space with the corresponding quadrants of ejections and sweeps (with permission from Pope (2000)); (b): A schematic picture of counter-rotating rolls in the near wall region (from Holmes et al. (1996)).

Contours of X Velocity (m/s)

FLUENT 6.1 (3d, segregated, lam) Apr 07, 2004 8.03e-01 7.69e-01 7.35e-01 7.01e-01 6.66e-01 6.32e-01 5.98e-01 5.64e-01 5.29e-01 4.95e-01 4.61e-01 4.27e-01 3.92e-01 3.58e-01 3.24e-01 2.90e-01 2.55e-01 2.21e-01 1.87e-01 1.53e-01 1.18e-01 Contours of X Velocity (m/s)

FLUENT 6.1 (3d, segregated, lam) Apr 07, 2004 8.03e-01 7.69e-01 7.35e-01 7.01e-01 6.66e-01 6.32e-01 5.98e-01 5.64e-01 5.29e-01 4.95e-01 4.61e-01 4.27e-01 3.92e-01 3.58e-01 3.24e-01 2.90e-01 2.55e-01 2.21e-01 1.87e-01 1.53e-01 1.18e-01

(45)

3.5. LES OF A TURBULENT PIPE FLOW 31

(a) (b)

Figure 3.9: (a): Velocity vector plot in the y − z plain; (b): Detail from the top wall region showing the pair of counter-rotating rolls.

3.5 LES of a turbulent pipe flow

Fully-developed pipe flow has also often served as a test case, both for the RANS and LES models. While the flow pattern is similar to the one in a channel flow, from the numerical point of view, the geometry is more complex. Detailed investigations of the pipe flow have been conducted among others numerically by Eggels (1994) and Friedrich and Nieuwstadt (1994), and experimentally by Durst and Wang (1989) and Toonder and Nieuwstadt (1997). Here, some results of the LES simulation of the pipe flow at two Reynolds numbers (Reτ = 590 and Reτ = 180) are presented. The aim of these simulations was to test the implementation of the dynamic SGS model. Additionally, the quality of the budgets of the transport equations for the Reynolds stresses and turbulent heat-flux components, calculated from the LES results, is investigated. This issue is relevant because the same procedure was later used to calculate the budgets of the transport equations in the LES results of the impinging jet flow.

z θ r U W V

(46)

3.5.1 Results for Re

τ

= 590

The dynamic Smagorinsky model was used for the simulation of the pipe flow at Reτ = 590. The model was implemented in the unstructured finite volume code. The negative values of Smagorinsky coefficient Cs, obtained from the dynamic procedure, were clipped to zero in order to prevent numerical instabilities. This eliminates the backscatter, but in simple equilibrium flows this is not a serious issue. The length of the pipe was 4R and periodic boundary conditions were imposed in the streamwise direction. Eggels (1994) has shown that the one-dimensional two-point correlations of the velocity and pressure fluctuations have very small values after L > 4R. The hybrid mesh, shown in Fig. 3.10, was used with hexahedral cells near the wall (r > 0.83R) and the tetrahedral prisms in the core region. The velocity components in the axial (z), radial (r) and azimuthal (θ) directions are denoted by W, U and V. The region covered by hexahedral cells had 380 cells in the azimuthal direction (θ), 70 cells were distributed in the axial (z) direction and the near-wall layer r+ _{< 35 had 20 cells clustered towards the wall. The total number} of cells was 801500. This produced ∆z+ _{= 33, and ∆r}+ _{= 0.7 and r∆θ}+ _{= 9.7 at} the wall (+ indicates normalization by the friction velocity Uτ and fluid viscosity). The constant initial velocity field was randomly perturbed every 50 time steps in the period of 1000 time steps. This helped to develop instabilities which led to a fully turbulent state. The time marching was performed using a fully-implicit three-level time scheme. The SIMPLE algorithm was used for coupling velocity and pressure. Central differencing was used for both the diffusive and convective terms, without any under-relaxation procedure. The flow was driven by the constant mass flux in the stramwise direction. Dimensionless time step, normalized by R and Wb, was 0.01 which produced CFL below 0.5. Statistics were collected over a period of 40 flow-through times.

The simulation results are shown in Fig. 3.11. U, V and W stand for the radial, azimuthal and axial components of the velocity vector. The mean-velocity profile agrees well with the experiment of Toonder and Nieuwstadt (1997). The appropriate log-law,

W+ ₌ 1

κln(Er

+_{), is found to yield κ = 0.4. The averaged C}

s at the central region of the pipe is found to be 0.07. This is close to the value found by Ferziger (1995) to be appropriate for the channel flow. The resolved Reynolds stresses, normalized with the inner-wall scales are in good agreement with the experiments. Since the experiment was at Re = 17000 and Re = 24000, the small differences in the stresses that appeared in the region away from the solid boundaries are due to the differences in the Reynolds numbers. The peak of the normal stress uu is well captured. The error in the prediction of the skin friction coefficient, calculated by Eq. 3.32, is 5%.

Fig. 3.12 shows the near-wall streaks in the plane parallel to the wall, at r+ _{= 6. The} flow structures are similar to one observed in the channel flow.

(47)

3.5. LES OF A TURBULENT PIPE FLOW 33 1 10 100 1000 r+ 0 5 10 15 20 25 30 W + (a) 0 50 100 150 200 r+ 0 0.5 1 1.5 2 2.5 3 w_rms(Exp.) u_rms (Exp.) w_rms (LES) u_rms (LES) (b) 0 0.1 0.2 0.3 0.4 0.5 r/D 0 0.2 0.4 0.6 0.8 1 uw + (c)

Figure 3.11: (a): Mean-velocity profile, solid line: LES, dash line: log-law, symbols: experiments of Toonder and Nieuwstadt (1997) Re = 24600; (b): Resolved-axial and -radial profiles of rms of the normal stresses, line: LES, symbols: experiments of Toonder and Nieuwstadt (1997) Re = 17800; (c): Resolved shear stress, solid line: LES, symbols: experiments of Toonder and Nieuwstadt (1997) Re = 24600.

(a) (b)

(48)

3.5.2 Results for Re

τ

= 180

A LES simulation of the pipe flow at Reτ = 180 was conducted in order to test the quality of the budgets of the transport equations for the Reynolds stresses and components of the turbulent heat-flux vector, calculated from the LES results. The small Reynolds number is chosen since DNS data for both the velocity and thermal field are available for this case. Since the impinging jets are homogeneous in the azimuthal direction, the transport equations are transformed from the Cartesian coordinate system into the cylindrical-polar coordinate system. The transformation rules and resulting equations are given in the Appendix 7.3. Furthermore, the influence of the mesh resolution on the budget accuracy is easier to test at smaller Reynolds number. The final goal was to use the LES results of the impinging jet flow to calculate the budgets of the transport equations.

The computation parameters were similar to the ones used in the previous LES sim-ulation of the pipe flow at Reτ = 590. Two different meshes were adopted. Both meshes had the same mesh resolution in the streamwise and spinwise directions (r∆ θ+ _{= 8} and ∆ z+ _{= 30) while the grid distribution in the wall-normal direction differed. The} coarse mesh had 18 points in the wall normal direction in the region r+ _{< 36 while the} fine mesh had 28 points in the same region. Two different boundary conditions for the energy equation were applied. A constant heat-flux condition was imposed on the walls in one case and a constant-temperature condition in another. The coordinates and flow variables are normalized by the channel width h, the friction velocity Uτ, the kinematic viscosity ν and the friction temperature Tτ (=

qw ρ cp Uτ

).

The flow-field results were compared with the DNS data of Fukagata and Kasagi (2002) while the results of the thermal-field simulation were compared with the DNS data of Kasagi et al. (1992). Kasagi et al. (1992) have simulated fully-developed thermal field in a turbulent channel flow for Reτ = 180. Since the basic physics does not differ significantly in a turbulent pipe and channel flows, the results of the Kasagi et al. (1992) were used for the comparison.

The mean-velocity and temperature profiles (for the latter very similar results are obtained for both types of the wall conditions), shown in Fig. 3.13, agree reasonably well with the DNS data. Only the results obtained with the finer mesh were plotted in cases when no noticeable differences existed between two different meshes. The results for the Reynolds stresses are shown in Fig. 3.14. Again, agreement with the DNS data is very good. The differences between two meshes start to appear in the results of the components for the turbulent heat flux vector, see Fig. 3.15. While differences between the results for u+_t+ _{are not significant, the peak of w}+_t+ _{is captured better by the} finer mesh. The influence of different boundary conditions for the energy equation can be seen in Fig. 3.16, which shows the distribution of the temperature fluctuation t+

LES, RANS and combined simulation of impinging flows and heat transfer

LES, RANS and Combined Simulation

of

LES, RANS and Combined Simulation of

Impinging Flows and Heat Transfer

Contents

Nomenclature

Roman

Greek

Subscripts

Abbreviations

Chapter 1

Introduction

1.1

Motivation

1.2

Objectives

1.3

Outline

Chapter 2

Numerical method

2.1

Introduction

2.2

Governing equations for heat and fluid flow

2.3

Discretization

2.3.1

Discretization in space

2.3.2

Integration in time

2.4

Data structure

2.5

Auxiliary numerical details

2.5.1

Implementation of bounded convective schemes for

un-structured grids

φ

φ

φ

φ

φ

φ

φ

φ

φ

U

C

j

D

φ

ξ

Chapter 3

Large-eddy Simulation

3.1

Introduction

3.2

Subgrid-scale modelling

3.2.1

Smagorinsky subgrid-scale model

3.2.2

Dynamic subgrid-scale model

3.3

Unstructured-grid Large Eddy Simulation

channel walls

2h

Flow direction

x, U

z, W

y, V

3.4

LES of a plane channel flow

3.4.1

Subgrid-scale model and computation details

3.4.2

Results

3.4.3

Channel flow structures

3.5