Direct simulation and regularization modeling of turbulent thermal convection

(1)

(2)

(3)

Direct simulation and regularization modeling

of turbulent thermal convection

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 12 maart 2007 om 10.00 uur.

door

Maarten VAN REEUWIJK

(4)

Dit proefschrift is goedgekeurd door de promotor: Prof. dr. K. Hanjali´c

Toegevoegd promotor: dr. H.J.J. Jonker

Samenstelling promotiecommissie:

Rector Magnificus voorzitter

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

dr. H.J.J. Jonker Technische Universiteit Delft, toegevoegd promotor Prof. dr. ir. B.J. Geurts Universiteit Twente

Prof. dr. D.D. Holm Imperial College London, United Kingdom Prof. dr. R.M. Kerr University of Warwick, United Kingdom Prof. dr. D.J.E.M. Roekaerts Technische Universiteit Delft

Prof. dr. K.R. Sreenivasan International Centre for Theoretical Physics, Italy

This work was financially sponsored by theStichting voor Fundamenteel Onderzoek der Materie (Foundation for Fundamental Research of Matter), FOM.

Copyright c 2007 by Maarten van Reeuwijk

All rights reserved. No part of the material protected by this copyright notice may be reproduced or utilized in any form or by any means, electronic of mechanical, including photocopying, recording of by any information storage and retrieval system without written permission from the author.

ISBN 978-90-9021652-2 cover design by Anymotion

(5)

v

SUMMARY

Direct simulation and regularization modeling

of turbulent thermal convection

This dissertation focuses on the simulation of turbulent thermal convection. The first part comprises a study of an ensemble of direct simulations of Rayleigh-B´enard convection in a domain with aspect-ratio four. The second part deals with a study of the Leray-α model for wall-bounded flows and the last part focuses on the effects of truncation errors on the energy cascade. The main objectives and results are summarized below.

Wind and boundary layers of Rayleigh-B´enard convection

The wind of Rayleigh-B´enard convection and its effect on the boundary layers is studied by direct numerical simulation in a Γ = 4 aspect-ratio domain with peri-odic sidewalls at Ra = {105_{, 10}6_{, 10}7_{, 10}8_{} and at Pr = 1. Symmetry-accounting} ensemble-averaging has been applied to decompose independent realizations into wind and fluctuations. It is shown that the characteristic peak in the squared mean horizontal velocities, is nearly entirely due to the wind. Both the wind and the fluc-tuations scale approximately as Ra0.5_{. Deep inside the thermal boundary layer,} horizontal heat-fluxes exceed the average vertical heat-flux by a factor 3 due to the interaction of the wind and the mean temperature field. The horizontal heat-transport causes spatial temperature differences, which creates buoyant regions that cause pressure gradients. The pressure gradients drive the wind, by which a feedback-cycle is complete.

From a study of the momentum and temperature budgets it follows that the boundary layers are clearly turbulent, although standard turbulence indicators like the shape factor and friction factor point towards laminarity. The reason for this ambiguity is that the flow is buoyancy dominated, and the induced shear is secondary. The joint effect of shear and buoyancy is characterized by the Obukhov length, and is found to scale at the same rate as the thermal boundary layer thick-ness. Overall, it is found that buoyancy forces significantly alter the structure of the boundary layer, and can therefore not be neglected in characterization and parametrization of the boundary layer.

Leray-α regularization for wall-bounded flows

(6)

vi Summary

107and Pr = 1, and the side-heated vertical channel at Ra = 5 × 106_and_{Pr = 0.7.} The simulations are compared to results from a well resolved and coarse DNS. It is found that for all three flow cases, the variability of the velocity field in-creases as the filter width parameter a is increased. Here, a is connected to the

filter width as α_i = a∆x_i, with ∆x_i the local grid size. Furthermore, the bound-ary layers tend to thicken relative to the coarse DNS results as a function of a. In

the cases where coarse DNS overpredicts wall gradients (as for Rayleigh-B´enard convection and the side-heated vertical channel), the thickening is benificial. How-ever, for the plane channel flow, coarse DNS underpredicts the wall-shear velocity, and increasinga only degrades the results. The inclusion of body forces into the

Leray-α model has to be done with care, as the coupling between the production of TKE by buoyancy and the heatflux can become indirect. Overall, we conclude that the performance of the Leray-α model in the present formulation is not con-vincingly superior to coarse DNS for the flow types under consideration. For these flows, additional (dissipative) modeling seems unavoidable to capture accurately all effects from turbulence.

The effect of truncation errors on the energy spectrum

(7)

CHAPTER

1

Introduction

Life on earth would look rather different without fluid turbulence. Typical mixing timescales would be in the order of decades instead of minutes. The atmosphere and oceans could not nearly maintain a temperature as uniform as now, resulting in much larger temperature differences. Provided that a hydrological cycle per-sists, rivers would flow at rates of the speed of sound. In all likelihood we would not exist or at best look very different, if there were no turbulence. Turbulence is caused by fluid elements transporting their properties such as linear momentum and internal energy along. Although this is a simple local process, it is nonlinear, which makes an accurate macroscopic description notoriously difficult.

The focus of this thesis is on turbulence generated by buoyancy and shear in simple domains bounded by two flat plates. Adopting a coordinate system aligned with the mean flow, where thez-direction represents the wall-normal direction, the

ratio between buoyancy- and shear-production of turbulent kinetic energy (TKE) is described by the flux-Richardson number

Ri_f = βg w0Θ0

w0u0_∂_zu (1.1)

where β is the expansion coefficient,g the gravitational acceleration, w0Θ0_the ver-tical turbulent heat-flux and w0_u0_∂_z_{u the shear-production term. The · operator} denotes a suitable (time, space or ensemble) average. This thesis solely deals with unstably stratified flows wherew0Θ0 >0, i.e. the thermal fluctuations act as a pro-duction term of TKE. The asymptotic limit of Ri_f → 0, represents forced convec-tion where all turbulence is produced mechanically by the straining of the fluid. Free (thermal) convection is the asymptotic limit ofRi_f → −∞.

The forced convection case Ri_f _{→ 0 represents turbulent plane channel flow,}

and is one of few real success-stories of turbulence, with its logarithmic law of the wall (e.g. Pope, 2000) ∗_{. At the other asymptotic limit of} _Ri

f → −∞ is Rayleigh-B´enard convection, which is generated when a destabilizing temperature differ-ence is enforced between the bottom and the top plate. Although Rayleigh-B´enard convection has been of great importance to chaos theory (Lorenz,1963, and many others afterwards) and the study of non-equilibrium systems (Glansdorff & Pri-gogine,1971), the asymptotic behavior at large Rayleigh numbersRa is still a point

of active discussion (e.g.Chavanneet al.,2001;Niemela & Sreenivasan,2003). One

∗_{Although even this is a topic of discussion.}_Barenblatt₍₁₉₉₃_{) argues that the overlap layer may}

better described by a power law.

(12)

2 Chapter 1. Introduction

of the main controversies is whether the heat-flux at very large Ra scales

propor-tional to Ra1/3 _or_Ra1/2_{, which is intimately connected with the way the velocity} field transports the heat. If the heat-flux becomes independent of the distance H

between the plates, a dimensional argument shows that the Nusselt number Nu

will scale as Nu ∝ Ra1/3 _(Priestley, _1954; _Adrian _{et al.}_, _{1986). If the boundary} layers communicate by the presence of a large-scale wind structure (e.g. Krishna-murti & Howard, 1981), H will remain a control parameter (by which the system

becomes independent of molecular properties) and a scaling ofNu ∝ Ra1/2_{can be} expected (Kraichnan,1962;Grossmann & Lohse,2000;Lohse & Toschi, 2003). The controversy reflects an insufficient knowledge of the exact heat-transfer mecha-nism of Rayleigh-B´enard convection. The wind and boundary layers are central in this process, and providing insight into these two topics is one of the main aims of chapter4and5.

After a description of the code for direct simulation (chapter 2) and its vali-dation (chapter 3), a detailed study is made for moderate aspect-ratio Rayleigh-B´enard convection (L/H = 4) at Pr = 1 and with Ra = 105_{− 10}8_{. In} particu-lar, a symmetry-accounting ensemble-averaging method is introduced in chapter 4 which can distinguish wind from fluctuations in a domain with periodic side-walls. Then, a detailed study is made of the wind structure and the boundary layers in chapter5.

The second part of the thesis focuses on regularization modeling of turbu-lent flows for largy-eddy simulation (LES) (Geurts & Holm, 2003, 2006b). The classical approach to large-eddy simulation is to model the unresolved subfilter scales by replacing them with a diffusive term. This approach is succesful away from boundary layers but proves to be quite cumbersome when applied to wall-bounded flows due to the strong anisotropy of the flow near the wall and the smallness of the large eddies. The regularization approach modifies the advective terms, which reduces the nonlinearity and therefore enhances the computability of the problem. In this way, the regularization acts as an subfilter model. The model we will study is the simplest of the regularization models, the Leray-α model, which replaces the advective transport velocity with a smoothed velocity.

(13)

3 The thesis is a collection of papers, either published or in preparation. There-fore, all chapters are self-contained, which makes a certain degree of redundancy inevitable. As the different papers address different audiences, the meaning of symbols may differ per chapter. The most important difference is that ˜u represents

a symmetry-accounted average in chapter 4 and 5 but ˜u represents the filtered

(14)

(15)

CHAPTER

2

Theory and numerical method

This chapter provides a brief overview of the governing equations of turbulent convection in the Boussinesq approximation. Several aspects of the numerical code are discussed and the discretization in space and time is presented.

2.1 Transport equations

The compressible Navier-Stokes equations are given by

∂tρ + ∂_ju_jρ = 0 (2.1)

∂tρui+ ∂jujρui+ ∂ip = ∂jµ∂jui+ fi+ ρgi (2.2) Here, ∂_j = _∂x∂

j and implicit summation over repeated indices is used unless stated

otherwise. These equations ensure the conservation of mass and linear momen-tum, respectively. The variables ρ,u_i and p are the density, velocity, and pressure,

respectively. The parameter µ represents the molecular viscosity, g_i the gravita-tional acceleration and f_ia body force.

The first law of thermodynamics states that changes in the internal energy ˆU

result either from the addition or removal of heat or by work done by or on the system. It is convenient to use the equation of enthalpy h = ˆU + p/ρ and its

transport equation is given by (Birdet al.,2002)

ρ∂_th + ρ∂_ju_jh = ∂_jλ∂_jT + ε_u+ (∂_tp + ∂_ju_jp). (2.3) Radiation effects and internal sources have been neglected here, and Fourier’s law has been used for the conductive energy transport. The temperature and thermal conductivity are denoted byT and λ, respectively. The source term ε_u represents the conversion of kinetic energy into heat given by εu = µ(∂_ju_i)2. For our purposes the dependence of h on the pressure can be neglected (which is exact for ideal

gasses) and the specific heat capacityc_p can be assumed constant, so thath = c_pT

up to an integration constant.

Using the continuity equation, decomposing ρ in a reference density and fluc-tuations as ρ = ρ0+ ρ0and dividing by the reference density ρ0, we obtain

(16)

6 Chapter 2. Theory and numerical method

Here the kinematic viscosity ν and the thermal diffusivity κ are given by ν = µ/ρ0 and κ = λ/(ρ0cp), respectively. Note that this formulation is still exact: no further assumptions have been made.

The Boussinesq approximation states that the density fluctuations ρ0 _{are small} compared to the reference density ρ0 so that ρ0/ρ0 1, and that all density vari-ations can be neglected except for the buoyant forcing in the momentum equation (Nieuwstadt,1998;Birdet al.,2002). Furthermore, dependencies of the fluid prop-erties β, µ andPr = ν/κ on temperature and pressure are neglected (Turner,1973). To estimate the typical magnitude of the various terms, we introduceU, L, t∗ _and ∆T as a typical velocity, lengthscale, timescale and temperature difference, respec-tively. For the continuity equation it follows that ∂_ju_j ≈ U/H is much larger than the term for density fluctuations, so that the latter term can be safely ne-glected. In the momentum equation, the Boussinesq approximation requires that (∂_t +u_j∂_j)u_i _g_i, i.e. the assumption is thatU2/H g. In the energy equation, the dissipation rate can be neglected as long asU2 _c

p∆T. For the same reason, the work done by pressure can be neglected. This results in

∂_ju_j =0, (2.4) ∂_tu_i+u_j∂_ju_i+ ∂_i p ρ0 − gjxj = ν∂2_ju_i+ fi ρ0 + ρ0 ρ0gi, (2.5) ∂_tT + ∂_ju_jT = κ∂2_jT. (2.6)

Here the pressure is independent of the density and temperature, reflecting an infinite speed of sound. The continuity equation has become diagnostic and the fluid is incompressible.

As we are interested in thermal convection, the density can be assumed to be purely dependent on temperature. A Taylor expansion of the density around the reference temperatureT0results in

ρ(T) = ρ0− ρ0β(T − T0)) +O(T2), (2.7)

where the expansion coefficient β is defined as

β =₋ 1 ρ₀ ∂ρ ∂T _T 0 . (2.8)

From (2.7) it follows that ρ0

ρ0 = −βΘ, where Θ = T − T0is the relative temperature.

Introducing a modified pressure ˜p = p/ρ0− gjxj and body force ˜fi = fi/ρ0, we arrive at the incompressible Navier-Stokes equations in the Boussinesq approxi-mation:

∂_ju_j =0, (2.9)

∂tui+uj∂jui+ ∂ip = ν∂˜ 2jui+ ˜fi− βgiΘ, (2.10)

(17)

2.2. Domain and boundary conditions 7

x

L

x

L

z

y

_L

y

T

0

T

1

g

Figure 2.1: Definition sketch of the domain under consideration

The tilde on the pressure and the body force will be dropped from this point on-ward.

2.2 Domain and boundary conditions

Underlying the work of this thesis is a numerical code for direct simulation. The coordinate system of the code is shown in Fig. 2.1. The domain is periodic in the

x and y directions and enclosed in the z-direction by two plates of length Lx and width Ly. The distance between the plates is Lz. The gravity vector is restrained to g_i _{= −g(sin φ, 0, cos φ), where g is the gravitational acceleration and φ is the}

angle with the z-axis. The body force f_i in the momentum equation is defined as

f_i =Pδ_i1, and P can be used to maintain a constant mass-flux if desired.

The boundary conditions on the sidewalls are periodic, i.e. X(x + nLx,y +

mLy,z) = X(x, y, z) for each integer n, m and field variable X. On the bottom and top plates, the velocity boundary condition can be either no-slip or free-slip. For the temperature, both a constant temperature and constant heat-flux can be enforced.

2.3 Numerical method

2.3.1 Spatial discretization

The equations (2.9-2.11) are discretized in a Cartesian coordinate system with a staggered grid-arrangement (e.g. Ferziger & Peri´c, 2002), as sketched in Fig. 2.2. The indices i, j, k are associated with the x, y, z directions respectively and scalars

(18)

8 Chapter 2. Theory and numerical method w_i,j,k v_i,j,k u_i,j,k c_i,j,k x_i−1 x_i+1 ∆x ∆zk c_i,k u_i,k w_i,k z_k+1 z_k z_k−1 x_i (a) (b)

Figure 2.2: Definition sketch of the staggered grid. a) the control volume and variable

arrangement; b) 2D vertical view of the grid and variables, where c can represent the

pressurep, relative temperature Θ or other scalars.

In the two homogeneous directions x and y, the grid is equidistant so that x_i = i∆x and y_j = j∆y. The code supports grid clustering in the wall-normal direction. Special attention has been given to the conservation of variance and ki-netic energy on a non-uniform grid by implementing a symmetry-preserving dis-cretization (Verstappen & Veldman,2003). The idea is that when the symmetries of the continuous operators are preserved upon discretization, important conserva-tion properties are conserved. For incompressible flows, the result of preserving the symmetries is that both momentum and energy is conserved on the discrete nonuniform grid. The symmetries of the advective (skew-symmetric) and diffu-sive operator (symmetric) are reviewed in appendixA.

The most important feature of the method is that advection does not take part in the total variance dissipation or production, and it will be shown below that this is achieved when the discretized advection operator is skew-symmetric. For con-venience we consider the transport of a scalarc by advection, although the same

argument holds for the kinetic energy. The vector c = [c_1,1,1,c_1,1,2...] represents the entire spatially discretized field, which evolves according to

dtc = Ac, (2.12)

whereA is the matrix of the discrete advective operator. Taking the inner product

with c on both sides results in an evolution equation of the total variance given by

d_t(1

2c · c) = cTAc (2.13)

(19)

2.3. Numerical method 9 the discrete advection operator A is skew-symmetric, it purely redistributes the

variance. If A has a symmetric contribution variance can either be dissipated as

is typical for upwind schemes or variance can be produced which may result in blow-up.

The code uses a second order symmetry-preserving discretization in space. Be-low, the discretized equations are given, where they-dependence has been omitted

for brevity. Continuity equation u_i+1,k− ui,k ∆x + w_i,k+1− wi,k ∆z_k = 0 (2.14) u-momentum equation du_i,k dt +

u_i+1/2,ku_i+1/2,k− ui−1/2,kui−1/2,k

∆x +

w_i−1/2,k+1u_i,k+1/2− wi−1/2,kui,k−1/2 ∆z_k

− ν  

ui+1,k−ui,k

∆x −

ui,k−ui−1,k

∆x

∆x +

ui,k+1−ui,k

∆z_k+1/2 −

ui,k−ui,k−1

∆z_k−1/2 ∆z_k   + pi,k− pi−1,k ∆x = P + βg sin φΘ (2.15) w-momentum equation dw_i,k dt +

u_i+1,k−1/2w_i+1/2,k_{− u}_i,k−1/2w_i−1/2,k

∆x +

w_i,k+1/2w_i,k+1/2_{− w}_i,k−1/2w_i,k−1/2

∆z_k−1/2 − ν   wi+1,k−wi,k ∆x − wi,k−wi−1,k ∆x ∆x + wi,k+1−wi,k ∆z_k − wi,k−wi,k−1 ∆z_k−1 ∆z_k−1/2   + pi,k− pi,k−1 ∆z_k−1/2 = βg cos φΘ (2.16) Temperature equation dΘ_i,k dt +

u_i+1,kΘ_i+1/2,k _{− u}_i,kΘ_i−1/2,k

∆x +

(20)

10 Chapter 2. Theory and numerical method

Here, ∆z_k−1/2 = z_k_{− z}_k−1 and ∆z_k = z_k+1/2 _{− z}_k−1/2 (see also Fig. 2.2). For completeness, it is noted that this discretization can also be obtained by a formal coordinate transformation to computational space, integrating over (now equidis-tant) control volumes and using central differences to estimate the fluxes.

Lastly, we would like to stress that energy conservation is a necessary but not sufficient condition for accurate solutions. Energy conservation ensures stable so-lutions but there are other numerical requirements, such as the cell-Reynolds num-berU∆x/ν, which needs to be sufficiently small to prevent wiggles from occuring.

Furthermore, turbulence is a delicate interaction between a large range of energy-exchanging scales. Due to truncation errors, the energy transfer is reduced at the high wavenumbers, which can also result in inaccurate solutions. The effect of truncation errors on the energy cascade are the topic of chapter8.

2.3.2 Temporal discretization

Defining a vector u_i,j,k = (u_i,j,k,v_i,j,k,w_i,j,k), the discrete momentum equations and continuity equation are given by

du_i,j,k

dt = Fi,j,k− ∇upi,j,k, (2.18)

∇cui,j,k=0. (2.19)

Here the vector Fi,j,k contains advection, diffusion and the body forces. The pres-sure gradient has not been included in F_i,j,k as it is required to make the velocity field divergence free at the next time step. Due to the staggering, two discrete di-vergence operators ∇uand ∇cexist; the former is defined at velocity locations and the latter at the control volume centers.

Time-stepping is performed by means of an explicit second order Adams-Bash-forth scheme, by which the system of equations becomes:

u_i,j,kn+1− u_i,j,kn ∆t = 3 2 Fn_i,j,k− ∇upn_i,j,k −1₂Fn−1_i,j,k − ∇upn−1_i,j,k (2.20)

2.3.3 Solving the pressure The pressure pn

i,j,k is used to make the velocity field divergence free on the time-leveln + 1. Taking the divergence ∇c of (2.20) gives that

(21)

2.3. Numerical method 11 Requiring that ∇cun+1_i,j,k = 0 results in a Poisson equation for the yet unknown pressure given by ∇c∇upi,j,kn = ∇cFni,j,k−1₃ ∇cF_i,j,kn−1− ∇c∇upn−1_i,j,k + 2 3 ∇cun_i,j,k ∆t . (2.22)

Note that the divergence of the velocity field at time-level n has not been

elimi-nated from the right-hand side. It has been kept there to be able to handle non-solenoidal initial conditions and to prevent the emergence of non-divergence-free velocity fields by the accumulation of round-off errors.

The Poisson equation (2.22) is solved by a direct method which performs Fast Fourier Transforms (FFTs) in homogeneous directions by which the equations de-couple in these directions and a differential equation in the z-direction remains.

This system can be solved mode by mode by a double-sweep algorithm. Details are presented in appendixF.

2.3.4 Enforcing a constant massflux

The body-force P in (2.15) can be used to maintain a constant average velocity

U0 in the x-direction, which is calculated as follows. The horizontal momentum equation of (2.20), after averaging over the entire volume h·iV, is given by

hun+1iV− huniV

∆t = hruiV+ 3

2Pn. (2.23)

Here hruiV contains all body and boundary terms, and the (unknown) pressure term vanishes by the Gauss divergence theorem (using the periodicity of the side-walls). Requiring that hun+1_i_V ₌ _U₀_{gives that}_Pn_{can be calculated by}

Pn = 2 3 U0− huniV ∆t − hruiV . (2.24)

(22)

(23)

CHAPTER

3

Validation of the code

In this chapter, the code for direct simulation is validated against a DNS database of turbulent plane channel flow and a side-heated vertical channel. Results are shown for the average velocity and temperature, Reynolds stresses, turbulent heat-flux and budgets of turbulent kinetic energy and temperature variance. Both the simulations of the plane channel and the side-heated vertical channel show an excellent match with the DNS database, although there are some small deviations, particularly for the side-heated vertical channel, which seem to be caused by differences in the numerics.

3.1 Plane channel flow

3.1.1 Case description

Turbulent plane channel flow is one of the standard flow cases to consider when validating a numerical model. Here we will use results from direct simulation per-formed by Moseret al. (1999) at a shear-Reynolds number ofReτ = uτδ/ν = 590. Here,uτ is the shear velocity, defined asuτ = pτw/ρ and the parameters δ, ν, τw and ρ represent the channel half-width, kinematic viscosity, wall-shear stress and density respectively. The shear-velocityuτand the viscosity ν are the characteristic scaling properties of this flow type, and from these two quantities, a lengthscale

ν/uτ and timescale ν/u2τ can be defined. Variables non-dimensionalized by these scales are refered to as plus-units: y+ ₌_yu

τ/ν andu+ =u/uτ.

For channel flow the coordinate system is usully defined as sketched in Fig.3.1, wherex, y, z are the streamwise, wallnormal and spanwise directions, respectively.

The medium is air with kinematic viscosity ν = 10−5 _m2_{/s and the domain size}

Lx× Ly× Lz is 2πδ × 2δ × πδ with the channel halfwidth δ = 1 m. The average

L

z

L

y

x

L

x

z

Figure 3.1: Sketch of the flow geometry and coordinate system for plane channel flow.

(24)

14 Chapter 3. Validation of the code 0 5 10 15 20 25 0.1 1 10 100 u + y+ present DNS (a) 0 1 2 3 4 5 6 7 8 0 100 200 300 400 500 600 σ 2 y+ present DNS u0_u0 v0v0 w0w0 (b)

Figure 3.2: Channel flow simulation atReτ = 610 compared against reference DNS of

Moseret al.(lines). a) average velocity profile for channel flow; b) variance profiles.

velocity is kept constant at U = 0.115 m/s, which results in a shear Reynolds

number of Reτ = 610. We use the same 384x384x256 resolution as Moser et al. (1999). However, our code has second order numerics whereas the code ofMoser et al. has spectral numerics. Thus it can be expected that differences will occur in the spectra and higher order moments. Grid clustering has been used in the wall-normal direction to ensure that 8 cells are present in the viscous sublayer belowy+ ₌_{5. All statistics presented are averaged over approximately 30 typical} turnover timest∗ _{= δ/}_{U and over homogeneous directions when appropriate.} 3.1.2 Average velocity, Reynolds-stresses and energy spectra

Shown in Fig.3.2a is the average streamwise velocityu in plus-units. The present

results are denoted by the circles and the reference DNS by the drawn line. Clearly visible is the viscous sublayer below y+ ₌ _{5, a buffer layer and a region with} logarithmic scaling fory+ _> _{30. The average velocity matches very well with the} reference DNS. The same holds for the velocity variances, which are shown in Fig. 3.2b, nondimensionalized byu2

τ. The present results are represented by the circles, squares and diamonds, and those of the reference DNS by the lines. The peak in the stream-wise variance u0u0 _{is captured well, and the same holds for the quick} drop-off after.

For plane channel flows, the turbulent shear-stress τt is a linear function of y (e.g.Pope,2000;Nieuwstadt,1998):

τt(y) ≡ ν∂yu − u0v0 =u2_τ y − δ

δ

(25)

3.1. Plane channel flow 15 0 0.2 0.4 0.6 0.8 1 0 100 200 300 400 500 shear -str ess y+ present DNS ν∂zu −w0u0

Figure 3.3: The turbulent shear-stress for channel flow. circles, squares: present DNS;

lines: reference DNS.

Close to the wall, τtis dominated by viscous forces and further away from the wall it is due to Reynolds stresses, as shown in Fig.3.3. The results (circles and squares) are in good agreement with the reference DNS (lines). Near the center of the chan-nelv0u0_{is slightly underestimated due to small transients on large timescales.}

Spatial energy density spectra at y+ ₌ _{39 (log-layer) and} _y+ ₌ ₅₈₇ (center-line) are shown in Fig.3.4. These energy spectra have been obtained by applying a one-dimensional Fast Fourier Transform (FFT) in one homogeneous direction and averaging over the other direction, as well as over time. As is usual for numerical codes, ∑_iEuu(ki) = u0u0, so that Euu represents the energy density up to a factor ∆k (the wavenumber increment). At the low and intermediate wavenumbers, the streamwise spectra (Fig.3.4a and c), are well predicted both near the wall and in the midplane. At the higher wavenumbers, the energy is first slightly overpre-dicted after which a steep drop-off can be discerned. This behavior seems typical for finite difference approximations, and led to the theoretical study of the effect of truncation errors on the energy spectrum in chapter8. The spanwise spectra are shown in Fig.3.4b and d. Aty+ ₌_{39, the low and intermediate wavenumbers are} well predicted, but the variance at high wavenumbers is overestimated, hinting at a slight underresolution. In the midplane y+ ₌ _{587, the energy densities are well} predicted for the entire range of scales.

3.1.3 Budgets for second order moments

(26)

bud-16 Chapter 3. Validation of the code 100 10−1 10−2 10−3 10−4 10−5 10−6 1 10 100 kx E _{present DNS}_E_uu Evv Eww MKM-590Euu Evv Eww (a) 100 10−1 10−2 10−3 10−4 10−5 10−6 1 10 100 kz E _{present DNS}_E_uu Evv Eww MKM-590Euu Evv Eww (b) 100 10−1 10−2 10−3 10−4 10−5 10−6 1 10 100 kx E _{present DNS}_E_uu Evv Eww MKM-590Euu Evv Eww (c) 100 10−1 10−2 10−3 10−4 10−5 10−6 1 10 100 kz E _{present DNS}_E_uu Evv Eww MKM-590Euu Evv Eww (d)

Figure 3.4: Spectra of streamwise and spanwise velocity components for channel flow. a)

y+ ₌ _{39 (log-layer), streamwise spectrum; b)} _y+ ₌ _{39, spanwise spectrum; c)}_y+ ₌ ₅₈₇

(27)

3.1. Plane channel flow 17 gets ofu0u0_,v0v0_andw0w0_{, of which the equations given by, respectively,}

0 = Pαα +Παα +Tαα −εαα

0 = −2v0u0_∂_yu +2p0_∂_xu0 _+∂_y_ν∂_yu02_{− v}0u02 _−2ν(∂_ju0₎2_{, (3.2)} 0 = +2p0_∂_y_v0 _+∂_y_ν∂_y_v02_{− v}0_v02_{− p}0_v0 _−2ν(∂_j_v0₎2_{, (3.3)} 0 = +2p0_∂_zw0 _+∂_y_ν∂_yw02_{− v}0w02 _−2ν(∂_jw0₎2. (3.4)

No implicit summation is performed over the Greek subscripts α, which denote the budget it corresponds to: T11 is the transport term ofu0u0, T22 is the transport term of v0_v0_{and so on.}

The budgets ofu0u0_,v0v0_, w0w0 _{are shown as a function of}y+ _{in Fig.}_3.5, nondi-mensionalized by u3

τ/ν. The production P11 of turbulent kinetic energy is in the streamwise component of the Reynolds-stress u0_u0_{, which is shown in Fig.} _3.5a. The production matches very well with the results ofMoseret al.(the dashed line), although the peak of the production seems marginally closer to the wall. The other components of the budget match well, although the dissipation rate ε11 and the transport at the wall T11 seem slightly overpredicted. The redistribution term Π11 drains energy from the u0u0 _{budget, which is given to the other two components} as Π11+ Π22 + Π33 = 0 by the incompressibility of the fluid. The wall-normal fluctuationsv0v0 _{are shown in Fig.}_{3.5b. Compared to the reference DNS (lines), it} can be seen that the values are slightly underestimated. However, we note that the vertical scale here is ten times smaller than in Fig.3.5a. The pressure redistribution term Π22 is the main source of production in the bulk, although it drains energy near the wall due to wall-blocking. Note that the dissipation ε22is zero at the wall due to the no-slip boundary conditions. The budget of the spanwise fluctuation

w0w0 _{is shown in Fig.} _{3.5c. Again, the profiles match very well with the results} fromMoseret al.(lines). The source of production is by the pressure redistribution Π₃₃, and the characteristic dip in the dissipation ε₃₃ is well predicted.

The budget of turbulent kinetic energy is obtained by summing (3.2-3.4) and dividing by 2, resulting in 0 = −v0u0_∂_yu | {z } P + ∂y(ν∂ye − v0e0− v0p0₎ | {z } T − ν(∂ju0_i)(∂ju0_i) | {z } ε (3.5) Here, e = u0

(28)

18 Chapter 3. Validation of the code -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0 50 100 150 budget y+ present DNS T11 P11 ε11 Π11 (a) -0.04 -0.02 0 0.02 0.04 0 50 100 150 budget y+ present DNS T22 ε22 Π22 (b) -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0 50 100 150 budget y+ present DNS T33 ε33 Π33 (c) -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0 50 100 150 budget y+ present DNS T P ε (d)

Figure 3.5: Budgets of turbulent budgets for channel flow compared to the DNS ofMoser et al. (lines). a) budget ofu0u0_{; b) budget of}v0v0_{; c) budget of}w0w0_{; d) budget of turbulent}

(29)

3.2. Side-heated vertical channel 19

z

y

Lz

Ly

L

x

T0

g

T1

Figure 3.6: Definition sketch of side-heated vertical channel.

3.2 Side-heated vertical channel

3.2.1 Case description

The side-heated vertical channel is a flow case where both buoyancy and shear production are important. The flow has several unusual features, such as neg-ative shear-production in the near-wall region and counter-gradient heat fluxes. The side-heated vertical channel has been studied intensively both experimentally (see e.g.Betts & Bokhari, 2000, and references therein)) and by direct simulation (Boudjemadi et al., 1997; Versteegh, 1998; Versteegh & Nieuwstadt, 1998, 1999). The flow geometry for the side-heated vertical channel is sketched in Fig.3.6. The coordinate system usually has thex-direction in the wall-normal direction and the z-direction pointing upward. For the study below, we will use the DNS database

fromVersteegh & Nieuwstadt(1999)∗_.

The code of Versteegh & Nieuwstadt is quite similar to our code for direct simulation. It has a staggered grid arrangement and uses second order central differences for the spatial discretization and second order Adams-Bashforth time-integration. The advection scheme is variance preserving (Piacsek & Williams,

∗_{The database is available at}

(30)

20 Chapter 3. Validation of the code

1970), although this may not be the case on non-uniform grids. There are two main differences between our results and those ofVersteegh & Nieuwstadt. The side-heated vertical channel has very long transients of about 100 typical timescales

t∗ ₌ _H/p_{βg∆ΘH. As these transients were quite inhibitive at that time,} Ver-steegh & Nieuwstadt introduced a pressure gradient which enforces a constant mass-flux which suppresses these transients, thereby reducing the required time to gather statistics. This additional body force has not been used for the present results. Furthermore,Versteegh & Nieuwstadt (1999) used Richardson extrapola-tion to construct soluextrapola-tions which are ’fully converged’ - results from simulaextrapola-tions with a coarse and fine grid were used to assemble a fully converged solution. The results from our DNS are for the fine resolution only.

We use an identical domain asVersteegh & Nieuwstadtof sizeLx× Ly× Lz =

H × 6H × 12H. The extent of this domain (aspect-ratio 12 in streamwise and 6 in

spanwise direction) is unusually large, but this size is required to prevent long-range spatial correlations from influencing the statistics (Versteegh & Nieuwstadt, 1998). The medium is air with a viscosity ν = 1.0 × 10−5 _m2_{/s, expansion} coef-ficient β = 3.3 × 10−3 _K−1_{, a Prandtl number} _{Pr = ν/κ = 0.709, and a channel} width H = 0.2 m. The temperature difference ∆Θ between the plates is 2.7 K

with T0 = ∆Θ/2 and T1 = −∆Θ/2. All statistics presented are averaged over approximately 120 typical turnover times and over homogeneous directions when appropriate. The resolution of the grid is 128 x 216 x 432 cells, with mild clustering to have sufficient resolution near the wall. This resolution is slightly higher than that ofVersteegh & Nieuwstadt, who use 96 x 216 x 432 gridcells. The extra reso-lution is not strictly necessary, but was used to check the grid independence of the solution.

3.2.2 Average quantities, variances and heat-fluxes

The side-heated vertical channel is driven by the buoyancy resulting from the heat-ing and coolheat-ing of the fluid on the sidewalls. This results in an upward flow near the hot plate and a downward flow near the cold plate, and the resulting average vertical velocity w is shown in Fig. 3.7a, nondimensionalized by the free-fall ve-locityU = pβg∆ΘH. The results from the reference DNS of Versteegh & Nieuw-stadt are represented by the lines, and those of the present DNS by the circles. The average velocity is about 5% higher than the reference DNS at the maximum velocity. Various components of the Reynolds stress are shown in Fig.3.7, nondi-mensionalized byU2 _{= β}_{g∆ΘH. The profiles of u}0u0 _(circles), v0v0 _{(squares) and}

w0u0 (plusses) match well with the results ofVersteegh & Nieuwstadt (lines). The profile ofw0_w0_{(diamonds) is slightly higher than the reference DNS.}

(31)

3.2. Side-heated vertical channel 21 0 0.1 0.2 0.3 0.4 0.5 0.6 0 0.1 0.2 0.3 0.4 0.5 w x/H present DNS 0 0.05 0.1 0.15 0.2 0 0.1 0.2 0.3 0.4 0.5 σ 2 x/H present DNS u0_u0 v0_v0 w0_w0 w0_u0

Figure 3.7: Average velocity (nondimensionalized by U = pβg∆ΘH and Reynolds

stresses (nondimensionalized by U2_{. lines: reference DNS; circles, squares, diamonds,}

plusses: present DNS. 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5 Θ , Θ 0Θ 0 x/H present DNS Θ 20 Θ0Θ0 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0 0.1 0.2 0.3 0.4 0.5 σ 2 x/H present DNS u0Θ0 w0Θ0

Figure 3.8: Average temperature (nondimensionalized by ∆Θ), temperature variance

(nondimensionalized by ∆Θ2_{) and turbulent heat fluxes (nondimensionalized by} _U∆Θ)

(32)

22 Chapter 3. Validation of the code

the ambient fluid due to the buoyancy (w0 >0), resulting in a positive stressw0u0. The average temperature Θ and temperature variance Θ0Θ0 match well with the results from the reference DNS (Fig.3.8a). The profiles are normalized by ∆Θ, and the values of Θ0Θ0 _{have been multiplied by 20 to bring them to the same} or-der as Θ. The turbulent heatfluxes are presented in Fig. 3.8b and are normalized by pβg∆Θ3_{H. The horizontal heat-flux u}0Θ0 _{is constant in the bulk of the flow,} whereas the vertical heat-flux w0Θ0 _{is nearly three times larger and has a mild} peak in the near-wall region. Compared to the reference DNS, u0_Θ0 _{is slightly} higher andw0Θ0 _{has the same typical value but is less peaked.}

3.2.3 Budgets for second order moments

The equations for the budgets of u0u0, v0v0 and w0w0 for the side-heated vertical channel are given by, respectively,

0 = Pαα +Παα +Tαα −εαα 0 = +2p0_∂_x_u0 _+∂_x_ν∂_x_u02_{− u}0_u02_{− p}0_u0 _−2ν(∂_j_u0₎2_{, (3.6)} 0 = +2p0_∂_yv0 _+∂_x_ν∂_xv02_{− u}0v02 _−2ν(∂_jv0₎2, (3.7) 0 = −2u_+2β0_gww0∂₀xw Θ0 +2p0∂zw0 +∂x ν∂_xw02_{− u}0w02 _−2ν(∂_jw0₎2. (3.8) As before, no implicit summation is performed over the Greek subscripts α, which denote the budget it corresponds to: T11 is the transport term of u0u0, T22 is the transport term ofv0v0_{and so on. The budget of turbulent kinetic energy is obtained} by summing (3.6-3.8) and dividing by 2, resulting in

0 = −u0w0_∂_xw + βgw0Θ0 | {z } P + ∂x(ν∂xe − u0e0− u0p0₎ | {z } T − ν(∂ju0_i)(∂ju0_i) | {z } ε (3.9) Here,e = u0

iu0iis the turbulent kinetic energy ande0 =u0iu0i.

(33)

3.2. Side-heated vertical channel 23 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0 0.1 0.2 0.3 0.4 0.5 budget x/H present DNS T11 ε11 Π11 (a) -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0 0.1 0.2 0.3 0.4 0.5 budget x/H present DNS T22 ε22 Π22 (b) -0.2 -0.1 0 0.1 0.2 0.3 0 0.1 0.2 0.3 0.4 0.5 budget x/H present DNS T33 P33 ε33 Π33 (c) -0.2 -0.1 0 0.1 0.2 0 0.1 0.2 0.3 0.4 0.5 budget x/H present DNS T P ε (d)

Figure 3.9: Budgets of turbulent budgets for channel flow compared to the DNS of

Ver-steegh & Nieuwstadt(lines). a) budget of u0u0_{; b) budget of} v0v0_{; c) budget of} w0w0_{; d)}

(34)

24 Chapter 3. Validation of the code -0.01 0 0.01 0.02 0.03 0.04 0.05 0 0.1 0.2 0.3 0.4 0.5 budget x/H present DNS u0_w∂xw βgw0Θ0

Figure 3.10: The two production terms of TKE for the side-heated vertical channel

com-pared to the reference DNS (lines).

The separate contributions of shear and buoyancy to the TKE production are shown in Fig. 3.10. The buoyancy production βgw0Θ0 _{is positive definite with a} small peak in the near-wall region. The shear-production −u0w0_∂_xw is a complex function ofx. In this flow case, u0_w0_{is positive definite. Therefore, the sign of the} shear-production is purely dependent on the average velocity gradient ∂xw. In the bulk, ∂xw is constant and negative to that −u0w0∂xw is positive; For 0.01 < x/H < 0.05, a region is present with countergradient fluxes, as u0_w0 _{and ∂}_x_{w have the} same sign. As a result, the shear-production term is actually negative, reflecting an energy transfer from the fluctuations to the mean flow. Very close to the wall (x/H < 0.01), there is a tiny region where the shear-production is positive again.

The small differences between the present results and those of Versteegh & Nieuwstadt may be caused by differences in the implementation of the advec-tive operator.Versteegh & Nieuwstadtuse an energy-conserving implementation, based on the skew-symmetric form of the advection operator (Versteegh,1998; Pi-acsek & Williams,1970). This formulation may not be fully energy-conserving on non-uniform grids. Furthermore, Versteegh (1998) reports small residuals up to 5% in the budget terms, whereas for the current results, the residuals are negligi-ble.

Lastly, the budget of turbulent temperature variance is given by 0 = −u0Θ0_∂_xΘ | {z } P + ∂x(ν∂xΘ0Θ0− u0Θ0Θ0₎ | {z } T − ν(∂jΘ0)(∂jΘ0) | {z } ε (3.10)

(35)

3.2. Side-heated vertical channel 25 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0 0.1 0.2 0.3 0.4 0.5 budget x/H present DNS T P ε

Figure 3.11: turbulent temperature variance budget for side-heated vertical channel

(36)

(37)

CHAPTER

4

Identification of the wind in Rayleigh-B´enard

convection

Using a symmetry-accounting ensemble-averaging method, we have identi-fied the wind in unbounded Rayleigh-B´enard convection. This makes it pos-sible to distinguish the wind from fluctuations and to identify dynamic fea-tures of each. We present some results from processing five independent three-dimensional direct numerical simulations of a Γ = 4 aspect-ratio domain with periodic side-boundaries at Ra = 107_{and Pr = 1. It is found that the wind}

boundary layer scales linearly very close to the wall and has a logarithmic re-gion further away. Despite the still noticeable molecular effects, the identifica-tion of log-scaling and significant velocity and temperature fluctuaidentifica-tions well within the thermal boundary layer clearly indicate that the boundary layer cannot be classified as laminar.

4.1 Introduction

Recent theories for the prediction of the scaling behavior of the non-dimensional heat fluxNu and the Reynolds number Re in Rayleigh-B´enard convection presume

the existence of ’wind’ (Shraiman & Siggia, 1990; Grossmann & Lohse, 2000), a large scale circulation that autonomously arises in the system. Indeed, this circula-tion has been observed by various groups and received quite some attencircula-tion over the last years (Krishnamurti & Howard,1981;Castainget al.,1989;Siggia,1994;Qiu & Tong, 2001;Niemelaet al.,2001;Sreenivasanet al.,2002;Parodiet al., 2004). The theories use wind primarily as a conceptual notion; little is known about the ac-tual magnitude of the wind as compared to velocity fluctuations and if it is present for all Rayleigh numbers (Sreenivasanet al.,2002). But most important, a straight-forward definition of ’wind’ seems to be missing. The natural candidates are: the ensemble average, a time average over a sufficiently long period, a spatial aver-age over a sufficiently large area, and a combination of the three. However, the problem is that all these averages yield a vanishing windu = 0. For the ensemble

average and the spatial average this may be obvious from symmetry considera-tions (see below), but for the time average it is more subtle. In this context the experiments (Niemela et al., 2001; Sreenivasan et al., 2002) are most instructive. Here it is shown that the ’wind’ erratically reverses its direction on time scales far exceeding the convective turnover time.Sreenivasanet al.(2002) propose a simple

Published inPhys. Fluids 17, 051704

(38)

28 Chapter 4. Identification of the wind in Rayleigh-B´enard convection

dynamical picture to explain these wind reversals, with two meta-stable states i.e. clockwise and counter-clockwise mean wind motion; only extreme events, such as (a combination of) energetic plumes can flip the system from one state to an-other. The relevant point here is that the reversals actually conserve theergodicity

of the system: the (long) time average becomes equal to the ensemble average — which is zero. The purpose of this Letter is therefore to come up with a suitable and useful definition of wind. Based on this definition we then determine the wind based on five direct numerical simulations (DNS) of moderate aspect-ratio and with periodic side-boundaries. An analysis of the wind field and fluctuations shows evidence of a boundary layer that is not laminar.

Rayleigh-B´enard convection is generated when a layer of fluid with thickness

H is subjected to a positive temperature difference ∆ between the bottom and top

wall. Within the Boussinesq approximation, the only control parameters are the Rayleigh numberRa = βg∆H3_(νκ)−1 _{and the Prandtl number}_{Pr = ν/κ. Here ν is} the kinematic viscosity, κ the thermal diffusivity and β the expansion coefficient. In domains with finite size, the aspect-ratio Γ = L/H, with L the horizontal extent

of the layer, is an additional control-parameter. The direction of gravity pointing in the negativez-direction, the equations for momentum, continuity and heat transfer

are given by

∂_tu_i+u_j∂_ju_i+ ρ−1∂_ip − ν∂2_ju_i = βgΘδ_i3, (4.1)

∂_ju_j =0, (4.2)

∂tΘ +uj∂jΘ − κ∂2jΘ = 0 (4.3)

with ρ the density,u_ithe velocity, Θ the temperature andp the pressure.

Since definitions for the processes occurring in Rayleigh-B´enard convection are not unambiguously defined, a small glossary is given here. We prefer to use the termconvective structure, which generalizes the terms wind and large scale

circula-tion, in that it involves both the velocity and the temperature field. This convective structure normally features convection rolls, which are the steady roll-like struc-tures. Plumes are the unsteady structures erupting from the boundary layers and propagating to the other side. Spatial averages will be denoted by hiAand hiH for plane- and height-averaging, respectively.

4.2 Accounting for symmetries

(39)

4.2. Accounting for symmetries 29 1 + + + + = u_i= 0 g

Figure 4.1: convection cell with zero mean flow. The zero ensemble mean solution is the

result of a superposition of conjugate symmetric modes.

symmetries will show up, although – due to the non-linear interactions – for large Rayleigh numbers only in an average sense (Frisch,1995).

Instead ofu_i =0 being the trivial solution, one may have an image in which the zero ensemble mean consists of groups of superimposed equi-probable conjugate symmetrical modes, as shown in Fig.4.1. Here a two-dimensional cell is sketched, together with some possible convective structures (modes). Focusing on the situ-ation with one roll, the clockwise and counter-clockwise modes are qualitatively and after a reflective operation also quantitatively identical. A more thorough and precise treatment of these arguments can be found in Chossat & Golubitsky (1988);Golubitsky & Stewart(2002).

Thus although mathematically correct, the ensemble mean u_i = 0 does not necessarily represent a physical mode of the system. Therefore we argue that for a useful interpretation of results it is important to account for these symmetries and present a generalized ensemble-averaging method that is able to do so. In classical ensemble-averaging the average is defined asX ≡ 1

N ∑Nα=1X(α)withN the number of realizations andX = X(x). To extend this to a symmetry-accounting

ensemble-averaging method, we have to apply some operations before ensemble-averaging:

e

X ≡ _N1

∑

N α=1

S(α)_X(α)_. _(4.4)

HereS(α) _{is an operator}_{S for realization α. The operator S is composed of one} or more elementary symmetry operatorsS(α) ₌_S₁_{◦ S}₂_{◦ . . ., such as horizontal or} vertical reflections. The elementary operators follow directly from the symmetries of the domain and boundary conditions. Classical ensemble averaging now re-duces to the case when S is the identity. Like in classical ensemble averaging, we

can decompose the fields in a symmetry-accounted mean and fluctuating part as

(40)

30 Chapter 4. Identification of the wind in Rayleigh-B´enard convection

4.3 Chasing the wind

Applying these concepts to our problem with periodic side-walls, the most impor-tant symmetry that must be accounted for is the translation invariance inx, y. This

operation is denoted by Sr with r ≡ (rx,ry) as the relative displacement. This would make the displacement r the only unknown per realization, but unfortu-nately the convective structure is not known eithera priori, which we address by

using an iterative pattern-recognition technique (Eiff & Keffer, 1997). With this technique the convective structure and the displacements are determined simul-taneously, gradually improving the estimation for the convective structure in suc-cessive iterations. The only assumption needed for this method is that – among all realizations – one and only one persistent structure (mode) is present inside the domain.

To start the iterative pattern-recognition process a reference fieldX0(x) is needed, for which a randomly picked realization is used – the convective structure is present in every realization so the starting point should not make a difference. Using a correlation functionC(X, Y), every realization can be compared to X0(x), and the location of maximum correlation is picked as the displacement vector:

d(α) _{← max}

r C(SrX

(α)_,_X

0). (4.5)

There is some freedom in choosing how to calculate the overall 2D correlation field, as it is constructed fromX ∈ {ui, Θ,p}. In this case we opted for the instan-taneous height-averaged temperature hΘiH which is closely related to the convec-tive structure as hΘiH > 0 wherew > 0 and vice versa. After calculating d(α) for all realizations and using (4.4), a new and improved estimation can be determined by e Xn+1 = _N1 N

∑

α=1 S_d(α)X(α). (4.6)

Repeatedly applying (4.5) with X0 replaced with Xn and (4.6) until eXn+1(x) = e

Xn(x) = eX(x) results in the convective structure, or symmetry accounted average (4.4), as well as the relative displacements d(α)_.

4.4 Results

(41)

4.4. Results 31

(a) (b)

(c) (d)

Figure 4.2: Development and growth of large scale structures in time as indicated by the

instantaneous height-averaged temperature hΘiH. Figure (a) ˆt = 12.5; (b) ˆt = 16.75; (c)

ˆt = 22.25; (b) ˆt = 37.5.

Ra = 1 × 107_, _{Pr = 1 and with an aspect-ratio Γ = 4 domain, the grid consisting} of 2563 _{cells. For this} _{Ra, the convective turnover-time based on the maximum} variance of horizontal velocity ist∗ ₌_{40 s; the non-dimensional time is defined as} ˆt ≡ t/t∗_{. In total five independent simulations have been performed.}

Figure 4.2 shows four successive time shots of one of the simulations for the height-average temperature hΘiH, which clearly show the persistence of the con-vective structure and its growth in time. At ˆt = 12.5 (Fig. 4.2a.), the flow is orga-nized into two up- and downward regions. Then follows an intermediate situation (Fig. 4.2b), resulting in a configuration with one up- and downward region (Fig. 4.2c and4.2d). The growth of convective structures has been observed before (Par-odiet al., 2004; de Roode et al., 2004; Hartlep et al., 2003), and goes on long after the process is statistically stationary. Although the details of the field are very unsteady as plumes rise and fall, the large-scale pattern is remarkably steady as Figs.4.2c and4.2d clearly show, indicating the presence of a persistent convective structure.

(42)

struc-32 Chapter 4. Identification of the wind in Rayleigh-B´enard convection (a) (b) λu λΘ h]w0_w0_i_A h ew ewiA h gu0_u0_i_A heu euiA k z/ H k/(βg∆H) × 102 5 4 3 2 1 0 1 2 0 −1₂ (c)

Figure 4.3: Symmetry-accounted wind field; a) streamlines at the edge of bottom thermal

boundary layer with contours of height-averaged temperature h eΘiH; b) the wind field

(43)

4.4. Results 33 1 0 _u_ˆ 1 0 λΘ 10 lo gˆz 10_{log ˆu} 0 -1 0 -1 λΘ 10 lo gˆz 10_{log ˆu} 0 -1 0 -1 z+ u+ 20 10 0 5 0 -5 ˆz ˆ u 1 0.5 0 1 0.5 0 (a) (b) (c)

Figure 4.4: Horizontal velocity profile in boundary layer: a) normalized as ˆu = eu/euλand

ˆ

z = z/λu with uλ(x) = u(x, λu(x)). Inset: non-dimensionalized with wall-shear stress.

b) log-log plot of ˆu, showing a linear near-wall profile. c) semi-log plot of ˆu, showing a

logarithmic region.

ture). The realizations are taken from five independent simulations at intervals of ∆ˆt = 0.5, from the moment the flow has developed to its largest scale. In the

lower boundary layer, the streamlines (Fig. 4.3a) clearly show an attracting re-gion where the flow is upward; one repelling rere-gion where the flow is downward and two saddle points. The height-averaged temperature h eΘiH (Fig. 4.3a) is con-sistent with this picture, as the relatively hot fluid is carried upwards and vice versa. Figure 4.3b shows a side-view of the average field after averaging over the y-direction. The contour lines are of relative temperature, which is the devi-ation from the plane-averaged temperature h eΘiA(z). Clearly visible in the figure is the projection of the two rolls onto the side view. This result may clarify why two-dimensional simulations are able to predict the Ra-Re-Nu behavior

reason-ably well (Kenjereˇs & Hanjali´c,2000).

(44)

dom-34 Chapter 4. Identification of the wind in Rayleigh-B´enard convection

inate in the vertical transport. In the horizontal, the average and fluctuating parts are of similar magnitude, with a striking difference that the squared mean hori-zontal velocity heueuiA peaks at the hydrodynamic boundary layer, but the mean squared velocity hgu0_u0_i_A _{is quite uniformly distributed outside the boundary} lay-ers, revealing a well-mixed bulk. Here we distinguish a thermal boundary layer

λΘ, defined as the location of maximum temperature variance hΘ]0Θ0iAand a hy-drodynamic boundary layer λu, defined by the location of maximum squared mean horizontal velocity heueuiA. It can be appreciated that even deep inside the thermal boundary layer λΘ(inset Fig. 4.3c), the mean squared horizontal velocity h gu0_u0_i_A_{is significant and of the same order as the squared mean velocity heueui}_A_.

Focusing on the hydrodynamic boundary layer, the inset of Figure4.4a shows the horizontal velocity profile non-dimensionalized with the wall shear-stress, whe-reu+ _{≡ u/u}

∗, z+ ≡ zu∗/ν, with the friction velocityu∗ ≡ pτs/ρ and local wall shear-stress τs. The velocity profiles are obtained by taking equidistant vertical cuts from the average field shown in fig. 4.3b. Clearly, scaling with the wall-shear stress does not yield a universal profile, indicating that the effects of buoyancy are not negligible - as would be expected from a buoyancy driven flow. However, all profiles do collapse upon normalization ofz by the local hydrodynamic boundary

layer thickness λu(x) and u by the wind maximum uλ(x) = u(x, λu(x)) (Fig.4.4a). The universal profile shows a linear near-wall profile (Fig. 4.4b), followed by a short region with logarithmic scaling (Fig. 4.4c).

Direct simulation and regularization modeling of turbulent thermal convection

Direct simulation and regularization modeling

of turbulent thermal convection

Maarten VAN REEUWIJK

SUMMARY

Direct simulation and regularization modeling

of turbulent thermal convection

CONTENTS

CHAPTER

1

Introduction

CHAPTER

2

Theory and numerical method

2.1 Transport equations

x

L

L

z

y

L

T

T

g

2.2 Domain and boundary conditions

2.3 Numerical method

CHAPTER

3

Validation of the code

3.1 Plane channel flow

L

L

y

x

L

z

z

y

Lz

Ly

L

x

T0

g

T1

3.2 Side-heated vertical channel

CHAPTER

4

Identification of the wind in Rayleigh-B´enard

convection

4.1 Introduction

4.2 Accounting for symmetries

∑

4.3 Chasing the wind

∑

4.4 Results

_L