channel flow in the PETSc toolkit

(1)

faculty of mathematics and natural sciences

Simulation of wall- bounded turbulent

channel flow in the PETSc toolkit

Master's thesis in Applied Mathematics

August 2011

Student: Leo van Kampenhout

Advisors: R.W.C.P. Verstappen, J. Top

Institute of Mathematics and Computing Science

(2)

(3)

Chapter 1 Introduction

Turbulence modeling is part of the broad scientific field of computational fluid dynamics, or CFD for short. The driving forces behind the research in this field can be roughly recognized into two categories. First, there is the wish for a deeper understanding of the underlying physics behind fluid dynamics. For well over a century now, turbulence has been one of the most intriguing and least understood phenomena in physics.

Secondly, because the engineering and medical applications of CFD are numerous, every advance in CFD leading to more accuracy can mean a potential business advantage. There exist many both commercial and academic CFD codes for all kinds of applications. What most of them have in common is that one or more modeling approximations or heuristics have been applied to enable them to deal with real-life problems. DNS¹ can be recognized as the only CFD instrument that does not make assumptions whatsoever about a particular flow. It is therefore considered equally good as (or even better than) real-life experiments.

Before the dawn of computers however, doing experiments was the only option available. Physicists set up experiments in order to gain a deeper understanding of fluids and their behavior. For example, Ludwig Prandtl discovered the existence of boundary layers while doing his experiments in 1904 [7]. Nowadays, experiments are still a very reliable and valuable tool for engineers in the design process. Probably every car sold in the last 30 years has been tested in a wind-tunnel to monitor its aerodynamic properties. However, experiments this scale require tremendous resources in time, money and technology, which are not always available. Also, for certain types of flow, experiments can only be done in scale (e.g. ships) or not at all (blood flow through the heart). These are examples where CFD is indispensable.

1Direct Numerical Simulation

(8)

Lx

h = 2δ

Lz

U0

Figure 1.1: Schematic view of the channel flow.

1.1 Problem setting

In this thesis we consider fully developed turbulent channel flow. The geometry for this common academic problem is quite simple, see Figure 1.1. Two perpendicular solid walls below and above the fluid are considered to have no-slip boundary conditions and have got distance h between them. The streamwise (x) and spanwise (z) directions are considered periodical and a linear pressure difference in imposed between the surfaces x = 0 and x = Lx

to act as an external forcing term. Since the flow is considered to be fully developed, the velocity statistics do not vary in the periodical directions:

they are statistically homogeneous.

This kind of channel flow can therefore be thought of as statistically one-dimensional and is in a way ’easier’ to calculate than a fully 3D inhomogeneous one. For example, the common approach here is to perform a (relatively cheap) Fast Fourier transform in the periodical stream and spanwise directions, thereby greatly simplifying the resulting Poisson equations.

Instead of having the fullO(N) complexity, the remaining part of the Pois- son matrix has only aboutO(N^1/3) complexity. This approach is called the spectral method and is used by [5], among others.

By contrast, the most difficult problems would arise in flows which do not have any homogeneous directions. These flows are common in many engineering applications, especially those with moving parts involved, such as the flow in the cylinder of a reciprocating engine ( [6], p338).

Despite the fact the channel flow is statistically one-dimensional, in this thesis we do not exploit this behavior and treat it as a fully inhomogeneous flow. The reader should mind this when comparing (timing) results with those using spectral methods.

(9)

1.2 Modern computer architecture

The aim of this section is to provide the reader with some basic background information about modern computer architecture. There are several distinct architectural types, and each brings with it some specific challenges. Very important to us, some are better scalable than others.

1.2.1 Multicores

Virtually all modern consumer PCs sold nowadays are equipped with so- called multicore CPUs, such as the Intel Core Duo. Next to usual improve- ments in the architectural structure of the core, increase in clock speed and better manufacturing processes, the multicore paradigm puts multiple copies of the same core together and let them act as one. The computational work directed to the CPU is distributed among the cores, so we might have a word processor running on the first core, a web browser on another whilst playing a Flash application on the third. We see that this approach is very suitable for a typical desktop environment, where we typically have separate and relatively small applications.

The downside of the multicore approach is that programs that are very computational demanding are still confined to only one core, effectively ren- dering the remaining cores of the CPU worthless. This is a very fundamental problem and can only be tackled when the program is, in some sense, “split- table” or parallizable. For instance, a computer game application is parallizable since it can be split into several components, such as ‘Video system’,

‘Audio system’, ‘Artificial Intelligence’ and ‘Physics Engine’. Therefore, computer games usually take well advantage of the ever increasing speed of multicore processors.

Applications that have not yet been made parallel are called serial and those that have parallel.

1.2.2 Shared vs. distributed memory

One of the most important distinctions between modern computer architectures is the way memory is used. The multicore architecture described above is an example of the shared memory approach, in which all the available cores use the same system memory, see Figure 1.2. The advantage of sharing the memory between processors is that it very much simplifies communication:

data computed by one processor is readily available to another.

The Symmetric Multiprocessing (SMP) architecture is very similar to the multicore approach in the way they deal with memory. However, instead of a single CPU with multiple cores, we have multiple CPUs on a single motherboard, called a node. Since the CPUs are located on the same physical circuit board, therefore having all the same distance to the memory,

(10)

core System memory core

core

Figure 1.2: Schematic view of the multicore paradigm, an example of the shared memory architecture. Similar is the SMP architecture, where the cores are no longer on the same chip, but rather form separate CPUs (with private caches).

node 1 memory

network bus CPU

node 2 memory CPU

node 3 memory CPU

node 4 memory CPU

Figure 1.3: Schematic view of distributed memory. Common variations involve individual nodes being multicore or SMP. The network can also have various layouts, often depending on the application.

which makes the shared memory approach preferable here.

None of this is the case with distributed memory, displayed in Figure 1.3.

In this approach each processor has its own private memory space, to which no other processor has direct access. The advantage of this approach is that the memory bus is not shared, so that each processor has full speed on its private memory. Also, processors need not be on the same node anymore, so it enables the construction of cluster-computers, in which multiple nodes are interconnect through a high-speed network, such as Ethernet or Infiniband.

Processors in a distributed memory environment only have access to their local memory and thus need to communicate with others in order to get global information. For example: during a typical CFD flow simulation, each processor owns only a certain partition of the domain, say [xs, xe] in the x-direction. In order to solve the differential equations that govern the dynamics of the flow, boundary values at xs− 1 and xe + 1 are required.

These have to be communicated from neighboring nodes over the network.

(11)

type number of nodes memory per node (GByte)

interactive (login) 1 128

batch 105 128 (87) or 256 (18 nodes)

test (system maintenance) 2 128

Table 1.1: List of nodes on the Huygens cluster.

1.2.3 Programming

Programming for a parallel architecture requires the use of one or more APIs. The leading APIs are:

OpenMP – the de-facto standard for shared memory programming and available on most architectures and operating systems.

MPI – The Message Passing Interface is the standard API for programming on distributed memory machines. This API takes care of the network-side of the communication (protocols, buffering, etc).

For this thesis however, PETSc is used, an high level numerical library specialized in advanced linear system solving. PETSc² is open-source, built upon the MPI standard and intended for use in scientific computations (so it is fast). By default it comes with a large set of Krylov solvers and a few basic preconditioners. Interfaces are present to install more advanced preconditioners, such as Hypre.

1.2.4 Huygens

All simulations featured in this thesis have been performed on the Dutch national supercomputer “Huygens”, maintained by SARA, the national center for High Performance Computing and e-Science.

The Huygens (IBM pSeries 575) is a clustered system of 108 SMP nodes.

Each node features 16 dual core processors (IBM Power6, 4.7 GHz, 64 bit) and either 128 GByte or 256 GByte. The total disk space amounts 700 TByte, divided between scratch space and home file systems. The total peak performance is 60 Teraflop/sec. Table 1.1 shows the distribution of the nodes.³

We are very grateful to NWO–NCF and SARA for the use of this fan- tastic machine.

2Portable, Extensible Toolkit for Scientific Computation. http://www.mcs.anl.gov/

petsc/petsc-as/

3All data retrieved from the SARA website, August 2011.

(12)

(13)

Chapter 2 Theory of Turbulence

This chapter reviews the most basic aspects of turbulence in an incompressible fluid. Starting off with the Navier-Stokes equations, it subsequently treats the Kolmogorov hypotheses and shows how energy is distributed over the Fourier modes.

2.1 Equations of motion

The flow of any Newtonian fluid is fully governed by the conservation laws of mass, energy and momentum, which are referred to as the equations of fluid motion. Since turbulence is merely a specific type of flow, it is also subject to these equations. All properties of turbulence can be led back to these equations and therefore they are the starting point of this chapter.

Conservation of mass The equation for conservation of mass, written in differential form, reads

∂p

∂t +∇ · (p U) = 0, (2.1)

where the pressure at point (x, t) is given by p(x, t) and the velocity vector by U (x, t). Throughout this thesis we will assume that the fluid is incompressible, in which case the energy equation is replaced by the condition of incompressibility. This condition enforces that the density of a particle is constant along it’s pathline, i.e.

Dp Dt ≡ ∂p

∂t + U · ∇p = 0, (2.2)

where _Dt^D denotes the material derivative. Combined with conservation of mass, the velocity field can be shown to be divergence free

∇ · U = 0. (2.3)

(14)

Conservation of momentum Using this result, the equation for conservation of momentum can be shown to be

DU Dt =−1

ρ∇p + ν∆U, (2.4)

where p(x, t) denotes the (modified) pressure and ν the kinematic viscosity coefficient. Equations (2.3) and (2.4) are commonly referred to as the Navier-Stokes (NS) equations.

2.2 Structure of turbulence

We recall that the (dimensionless) Reynolds number is given by Re = UL

ν , (2.5)

with ν the kinematic viscosity,U the largest velocity and L the largest length scale in the system. The Reynolds number quantifies the ratio between viscous forces and inertial forces. It is here important to realize that the Reynolds number is a property of the flow and not of the particular fluid involved. A high Reynolds number indicates that the inertial forces are strong, and turbulence starts to occur.

2.2.1 The energy cascade

In this chapter we will only consider fully turbulent flow, which has a very high Reynolds number. The structure of a fully turbulent flow is complex and consists of large-scale motion, which we call large eddies, and smaller scale structures. The large eddies in the flow are in size comparable to L, the smallest eddies can be as small as molecular scale.

Energy usually enters the system on the largest scale, which is for example the case for the driven cavity, wake and jet flow. The same amount of energy is removed from the system by means of dissipation, given that the total amount remains constant. Dissipation is known to be proportional to (size)⁻², which leads to a increase in energy in the largest eddies, since they are not very good in dissipating energy. This process renders the large eddies unstable and results in their breakup. The energy is thereby transported down to smaller scales. Eddies on the smaller scales will undergo the same destabilization and the process will repeat itself to yet smaller scales, until the dissipation effect is strong enough to prevent breakups. This concept, of the over and over again breakup of eddies, was first proposed by Richardson in 1922 and is called the energy cascade.

The Reynolds number of the smaller eddies decreases, until the eddy motion becomes stable and the kinetic energy is dissipated through viscous forces. This results in (some) thermal energy, i.e. the fluid heats up. The rate of dissipation is denoted as ε.

(15)

2.2.2 The Kolmogorov hypotheses

To investigate what precisely happens at the smallest scales of motion, we present the theory of Kolmogorov. In a general flow, the largest eddies are anisotropic. The first hypothesis of Kolmogorov states that this direction- dependent information is lost in the energy cascade.

Hypothesis 1 (Kolmogorov hypothesis of local isotropy). At sufficiently high Reynolds number, the small-scale turbulent motions are statistically isotropic. [6]

Taken this hypothesis to be valid, we can make a cut somewhere in the range of length scales, splitting the anisotropic eddies from the isotropic eddies. The former region which containing the largest eddies is called the energy-containing range, which suggest that most energy is stored here. This we will investigate later on. The region containing the small eddies is called the universal equilibrium range.

Kolmogorov continues and claims that all geometric information about the large eddies (i.e. their shape) is lost as well in the energy cascade, which implies that small-scale turbulence in any flow is in a sense universal. Of course this micro-turbulence still depends on the rate of energy transfer into the universal equilibrium range and the viscosity. In [6] (equation 6.13) it is shown that there exists a global balance between incoming energy at the largest scales, and the dissipation at the small scales on the other hand. In other words, the rate of energy transfer equals the dissipation rate ε.

Hypothesis 2(Kolmogorov’s first similarity hypothesis). In every turbulent at sufficiently high Reynolds number, the statistics of the small-scale motions have a universal form that is uniquely determined by ν and ε.

From these two parameters we can deduce, using dimensional analysis, unique length-, time and velocity scales called the Kolmogorov scales. For instance, the Kolmogorov length scale η is given by

η = (ν/ε)^1/4. (2.6)

The Kolmogorov length scale is the smallest scale in the energy cascade.

In this region energy is no longer transported downwards, instead all energy is being dissipated. The Reynolds number based on these scales is unity. A consequence of these hypotheses is that on the small scales, all high-Reynolds-number turbulent velocity fields are statistically similar; that is, they are statistically identical when they are scaled by the Kolmogorov scales.

Until now, we have talked about the energy cascade by which the energy passes from the largest flow scale (L) down to the Kolmogorov length (η). These length scales usually differ a number of orders of magnitude.

(16)

Energy-containing

Inertial subrange Dissipation range

Universal equilibrium range

η `DI `EI `0

range

L Figure 2.1: The ranges of scales that Kolmogorov distinguishes. Note that the

x-axis is shown in log-scaling.

Therefore, there exists a range of scales ` which are small compared to L but large compared to η. These scales are within the universal equilibrium range so the first similarity hypothesis applies, but dissipation hardly occurs.

It can be thought of that these scales merely transport the energy down the cascade. The inertial forces dominate and viscosity does not play a role, therefore this range is called the Inertial subrange. The third Kolmogorov hypothesis claims the existence of this range.

Hypothesis 3 (Kolmogorov’s second similarity hypothesis). In every turbulent flow at sufficiently high Reynolds number, the statistics of the motions in the range L ` η have a universal form that is uniquely determined by ε and `, independent of ν. [6]

The partitioning of scales into the ranges we have seen until now is pictured in Figure 2.1. The length scale where the energy-containing (E) range ends and the inertial (I) range starts is denoted by `EI. Likewise, the length scale between the inertial range and the dissipation (D) range is denoted `_DI.

2.3 Energy distribution

An interesting question is how the energy is distributed over the various length scales. This enables one to calculate the relative importance of the small length scales and estimate the error when these scales are modeled or omitted altogether.

2.3.1 The energy spectrum function

Before addressing this question, the following theory and notation needs to be introduced. Let U (x, t) denote the time-dependent velocity field. Then

(17)

the Reynolds decomposition into the mean and fluctuation is given by u(x, t)≡ U(x, t) − hU(x, t)i, (2.7) where h·i denotes the statistical mean. The one-point, one-time covariance of the velocity (the Reynolds stresses) are given by

hui(x, t)u_j(x, t)i. (2.8) The field U (x, t) is statistically homogeneous if all statistics are invariant under a shift of position. The two-point correlation is the two-point, one time auto-variance statistic

Rij(r, x, t)≡ hui(x, t)uj(x + r, t)i. (2.9) For homogeneous turbulence the two-point correlation is independent of x and the information it contains can be transformed into the wavenumber space. The velocity spectrum tensor is exactly that transformation:

Φ_ij(κ, t) = 1 (2π)³

∞

Z Z Z

−∞

e^−iκ·rR_ij(r, t) dr, (2.10)

and has an inverse transform which reads R_ij(r, t) =

∞

Z Z Z

−∞

e^iκ·rΦ_ij(κ, t)dκ. (2.11) Setting r = 0 yields the following equality

R_ij(0, t) =

∞

Z Z Z

−∞

Φ_ij(κ, t)dκ =huiu_ji. (2.12) This reveals that Φ_ij(κ, t) is the contribution to the covariancehuiu_ji from the velocity modes with wavenumber κ. The velocity spectrum tensor energy spectrum function Φij(κ, t) is a second-order tensor function of a vector and therefore contains a big deal of information. A more practical quantity is the energy-spectrum function E(κ) which may be viewed as Φij(κ, t) stripped of all directional information. This is done as follows:

E(κ, t)≡

∞

Z Z Z

−∞

1

2Φii(κ, t)δ(|κ| − κ)dκ (2.13) Integration over all scalar wavenumbers κ yields

Z ∞ 0

E(κ, t)dκ = 1

2R_ii(0, t) = 1

2huiu_ii, (2.14) where we identify the right hand size as the turbulent kinetic energy. We observe that E(κ, t)dκ represents the contribution to the turbulent kinetic energy from all modes with vector length κ. If the turbulence is isotropic, the velocity spectrum tensor is completely determined by E(κ, t).

(18)

2.3.2 The −⁵₃ spectrum

The question at hand concerns the distribution of energy over the various length scales. Kolmogorov uses the second similarity hypothesis to predict the form of the energy spectrum function E(κ) in the inertial range [6].

Dimensional analysis shows that the famous Kolmogorov−⁵₃ spectrum

E(κ) = Cε^2/3κ^−5/3, (2.15)

is the only way of ε and κ could form the universal spectrum. C is a dimensionless constant called the Kolmogorov constant. In this notation, dependence on time is omitted for convenience.

An interesting question is how the energy is distributed over the various length scales. Recalling that motions of length scale ` correspond to wavenumber κ = 2π/`, this question is equivalently answered when the distribution over the wavenumbers is known. This makes life easier, since the kinetic energy in the range (κa, κb) is the integral

k_(κ_a_,κ_b₎= Z κb

κa

E(κ)dκ. (2.16)

The total kinetic energy is denoted as k =

Z ∞ 0

E(κ)dκ, (2.17)

cf. Equation (2.14). It can be shown [6] that the amount of energy decreases as the wavenumber increases. This is consistent with the assumption that the bulk of the energy is contained in the large scale eddies. Furthermore, it can be shown that contributions to the dissipation rate ε is also dependent on the wavenumber and decreases as ε∼ κ^4/3 as κ decreases towards zero.

(19)

Chapter 3 Simulation of turbulence

Direct numerical simulation (DNS) is both the most straightforward way of simulating turbulent flow and the best one, since it does not involve any model for the turbulence. All physical scales are resolved exactly. However, this requires an extremely fine grid that can capture even the smallest scales of motion. Veldman demonstrates that the computational complexity for a fully 3D simulation scales with Re^11/4. A Reynolds number 10 times larger would require almost a factor 1000 more computing power [9]. DNS is therefore inapplicable to very high Reynolds-number flows. Other methods have been invented that make up for this practical constraints, such as LES.

3.1 LES

Large-eddy simulation (LES) has been invented to avoid the enormous computation costs that come with DNS. It is therefore more practical for most engineering problems and is widely used. In LES only the largest scales of motion, which contain most of the energy, are resolved. For the small scales a model is used, which should reduce the computational cost while main- taining the universal properties of the small-scale turbulence. There are the four conceptual steps in LES.

(1) The velocity U (x, t) is filtered such that the large eddies motion is still captured in the filtered velocity U (x, t) and the residual velocity u⁰(x, t) contains the smaller scales. The filter operation can be chosen freely, but it should at least meet some properties of locality, as we will see later.

(2) The Navier-Stokes equations are used to set up the evolution equations. In this equations the residual motions are captured in the so- called residual-stress tensor.

(3) A model is used in order to resolve the residual-stress tensor. This provides closure for the evolution equations.

(20)

(4) The flow equations are solved numerically for U (x, t), which gives us one realization of the flow. The numerical method is chosen with care to avoid unexpected side-effects. For instance, when both the model and the numerical method introduce artificial dissipation, the total amount of dissipation becomes unpredictable and difficult to control.

There are several new parameters introduced by this approach, such as the boundary-value between ’large’ and ’small’ scales, the filter operation used to remove the small scales from the flow, the particular model used for the small scales and it’s accompanying tuning parameters.

3.1.1 Filtering

With LES we do no longer require an extremely fine grid, since the filter operation eliminates the small-scale motions. A filter which does that is called a low-pass filter. The minimum grid spacing that will now suffice for the filtered velocity field is proportional to the specified filter width

∆. In order to resolve the most interesting motions, which are the energy- containing motions, we would want filter width to be smaller than `_EI.

The general filtering operation is defined by U(x, t) =

Z

G(r, x)U (x− r, t)dr (3.1) Z

G(r, x)dr = 1, (3.2)

where integration is over the entire flow domain. For the sake of simplicity we will restrict ourselves for now to the one dimensional case. The one- dimensional filters presented all have a three-dimensional counterpart, for which similar theory can be derived. Furthermore, we assume that the filter function is homogeneous, i.e. independent of x. For a random scalar field U (x) the filtering operation will therefore reduce to

U (x, t) =

∞

Z

−∞

G(r)U (x− r, t)dr. (3.3)

We can now define the transfer function bG(κ) as 2π times the Fourier transform of the filter function, which enables us to write the Fourier transform of the filtering operation as

U (x, t) =ˆ F{U(x)} = bG(κ) ˆU (κ), (3.4) where we use the convolution theorem for Fourier transforms.

In Figure 3.1 the effect of the Gaussian filter on the one-dimensional spectrum E11(κ) becomes apparent. The high wavenumbers, which correspond to the small-scale motion, are no longer contributing to the spectrum,

(21)

Figure 3.1: The one-dimensional spectrum E11(κ) (solid line)

i.e. they are filtered out. The filtered spectrum is obtained as follows, E₁₁(κ) =| bG(κ)|²E₁₁(κ). (3.5) Thus, filtering with filter function G results in attenuating the energy of mode κ by a factor bG(κ)². These attenuation factors are shown in Figure 3.2 for three commonly used filters. The sharp spectral filter appears to be most effective at attenuating high wavenumber modes, but is on the other hand non-local in physical space. The Gaussian filter is compact in both physical and wavenumber space.

3.1.2 Filtered equations of motion

Applying the filtering operation defined by Equation (3.1) to the Navier- Stokes equations leads to a new system of equations expressed in U (x, t) and p(x, t). For a spatially uniform filter, differentiation and the filtering operation commute, through which we can write the filtered continuity equation as

∂U_i

∂x_i

= ∂U_i

∂x_i = 0. (3.6)

cf. Equation(2.3). Note that the Einstein or suffix notation is used here, i.e.

∇ · U = ∂Uⁱ/∂xi. The filtered momentum equation is similarly shown to be

∂U_j

∂t +∂U_iU_j

∂xi

= ν ∂²U_j

∂xi∂xi −1 ρ

∂p

∂xj

, (3.7)

(22)

Figure 3.2: Attenuation factors bG(κ)²for these filters: the box filter (dashed line), the Gaussian filter (solid line), sharp spectral filter (dot-dashed line).

cf. Equation (2.4). This filtered equation introduces a new troublemaker, namely the filtered product UiUj, which is different from the product of the filtered velocities U_iU_j. This difference is called the residual-stress tensor

τ_ij^R≡ UiU_j− UiU_j (3.8) and can be split into an isotropic part ²₃krδij and a anisotropic part τ_ij^r. One can easily get rid of the former part by including it in the pressure, then called the modified filtered pressure, ¯p≡ p + ²₃kr. The filtered momentum equation then becomes

D Uj

Dt = ν ∂²Uj

∂x_i∂x_i −∂τ_ij^r

∂x_i − 1 ρ

∂ ¯p

∂x_j. (3.9)

In order to solve these equations for U (x, t) and ¯p(x, t) we need an estimate of the anisotropic residual stress tensor τ_ij^r(x, t). In LES, this estimate is provided through the use of a residual stress model, or subgrid- scale (SGS) model. This model will in general depend on the type of filter used, and the filter width.

3.1.3 The Smagorinsky model

The simplest model for τ_ij^r is the Smagorinsky model, proposed by meteo- rologist Joseph Smagorinsky in 1963.

(23)

It is based on the linear eddy-viscosity model, which expresses the hypothesis that the residual stress is linearly proportional to the filtered rate of strain

τ_ij^r =−2νrS_ij, (3.10) where S_ij denotes the filtered rate of strain tensor

Sij = 1 2

∂Ui

∂xj

+∂UJ

∂xi

. (3.11)

We will be using the Smagorinsky model in Chapter 5, when we will compare LES results to direct numerical simulations.

3.2 Other methods

Other methods have been proposed that (partially) solve the Navier-Stokes equations, are less computational expensive than DNS and still give very good results in some application areas. They fall outside of the scope of this thesis, but just to name a few:

Linear eddy viscosity models, such as the Cebeci-Smith model;

Two equation models, such as k-epsilon models;

Nonlinear eddy viscosity models;

Reynolds stress model. These four are all based on the Reynolds averaged Navier-Stokes equations (RANS). Lastly there is,

Hybrid methods between RANS and LES. This is alternatively called Detached eddy simulation (DES).

An excellent starting point to learn more about the differences between these models is the “CFD Online” wiki¹.

1http://www.cfd-online.com/Wiki/Turbulence_modeling

(24)

(25)

Chapter 4 Discretization

Until now we have treated the Navier-Stokes equations 2.3–2.4 as expressions involving continuous fields. While this the classical and most natural way to look at it, continuous fields are impractical for our purpose of computer simulations. Since the fundamental building blocks of computer hardware are discrete on–off switches called bits and we have only finitely many of them, it is theoretically impossible to represent a continuous field in the memory of a computer. (Basically for the same reason why it is impossible to represent an interval of real numbers with a finite amount of integers.)

This chapter is concerned with the approximation of continuous fields as discrete fields, intended for use on computer hardware. These discrete fields do not contain equally much information as the continuous fields, rather they only take values at certain points in space. The process of turning a continuous field into a discrete fields is called sampling, which is in itself not a terribly difficult task. It needs a lot of thought however, of how to represent operators (such as gradients) on these discrete fields. We need these operators in order to time-step the momentum and continuity equation.

4.1 Spatial discretization

4.1.1 Positioning of variables

In this work we will restrict ourselves to the finite-volume discretization on a structured or regular grid. This inevitably imposes certain restrictions on the geometry of the domain, e.g. walls can only have right angles and curved walls are ruled out. Methods exist however how to transfer the ideas in this section to unstructured grids, e.g. the concept of the shift operator idea which enables conversion of collocated-mesh operators to a unstructured, staggered grid [4].

The ‘location’ of variables is relative to the gridlines. In order to avoid the red-black pressure decoupling problem as described in [9], we follow the

(26)

vi,j

pi,j

ui−1,j

vi,j−1

ui,j

Figure 4.1: Location of the velocities u, v and the pressure p on the staggered grid. The control volume Ω is depicted as the solid line.

Marker-And-Cell (MAC) or staggered grid scheme. This scheme has some additional advantages when using the finite-volume approach. Velocities are by default known on the faces of a cell and need not be interpolated anymore.

This of course depends on the control volume at hand. In Figure 4.1, the location of the variables in two dimensions is shown. One minor disadvantage of the staggered grid is that boundary conditions do not always correspond directly with the location of the variables. However, the introduction of mirror points is a very elegant way to deal with this inconvenience [9].

4.1.2 Symmetry-preserving operators

The idea propagated by Verstappen and Veldman [8], among others, is that the discrete operators should mimic the continuous operators as closely as possible. To put it more precisely, the discrete operator should have the same symmetry properties as the underlying differential operator. It turns out that this kind of symmetry-preserving discretization has some very favourable properties:

it is stable on any grid;

it preserves total mass and momentum;

it preserves kinetic energy, when dissipation is turned off.

The resulting semi-discrete¹ system of equations can be written as follows:

Ωdu

dt + Cu + V u + Gp = 0, (4.1)

Du= 0, (4.2)

cf. Equations 2.3–2.4. In this discrete notation, Ω is a diagonal matrix representing the sizes of the control volumes and G is the discrete gradient

1Time is not discretized yet.

(27)

operator. The convective fluxes are represented by the matrix C, which is skew-symmetric by construction (i.e. C^∗ = −C). The diffusive (viscous) fluxes are represented by the matrix V , which is symmetric positive definite.

For the construction of these operators we refer to [8]. Please note that since the variables are defined on different spatial locations, their control volumes are also different.

The above discretization is stable when the kinetic energy norm, defined as ||u||² = u^∗Ωu, does not increase over time. Substituting Equation 4.1 into the evolution equation for the energy norm, we have

d

dt(u^∗Ωu)=−(Cu + V u + Gp)^∗u− u^∗(Cu + V u + Gp)

=−u^∗(C + C^∗)u− p^∗G^∗u− (G^∗u)^∗p− u^∗(V + V^∗)u

=−2u^∗V u≤ 0. (4.3)

The skew-symmetry of C makes sure the convective terms vanish. Fur- thermore, we use that the adjoint of the gradient equals the divergence, G^∗ =−D. Combined with 4.2, the terms with p vanish.

Now, what this exercise tells us is, energy is preserved if and only if the flow is inviscid i.e. the diffusion is turned off. For viscous flow, the energy is strictly decreasing, which means it is stable. In their article, Verstappen and Veldman show that it is possible to construct these operators on both uniform and non-uniform grids.

4.1.3 Accuracy and higher order methods

The standard approach in determining the accuracy of a discretization is by looking at the local truncation error. Second order behaviour in the local truncation errorO(h²) will lead to second order accurate global truncation error. However, this is merely a lower bound and second order accuracy can be achieved with first-order local truncation error schemes, too. Ver- stappen and Veldman show that this is true for the symmetry-preserving discretization, which has first order local truncation error [8].

An higher order discretization is obtained by applying Richardson extrapolation on two differently sized control volumes, see Figure 4.2. The largest one has 3^d times the volume of the smaller, with d the spatial dimension. For the convection C and divergence D, the same (second-order accurate) discretization is applied on the control volumes. Let C₁ and C₃ denote the convection fluxes on control volumes Ω1 and Ω3, respectively.

Then the following linear interpolation can be made:

C = 3^2+dC₁− C3. (4.4)

For the diffusion, we will also use the most straightforward approach

V = 3^2+dV₁− V3, (4.5)

(28)

Ω1

Ω3

Figure 4.2: Superposition of two control volumes for the construction of a fourth- order discretization. In this case the velocities are defined on the control faces, so no interpolation is needed. Shown is the control volume corresponding to the pressure.

Note that this differs slightly from the operator proposed by Verstappen and Veldman [8] and no longer mimics the exact symmetry properties of the continuous operator.

4.2 Time integration and stability

Preferably, the discretization of the time derivative is done is such a way that the energy is preserved in time. In this way the timestep ∆t would be of no influence on the stability, hence the discretization would be uncon- ditionally stable. The only way this can be achieved however, is with an implicit scheme. Implicit time schemes require a huge amount of work to be done per time step and to ensure accuracy, many timesteps have to be done. Therefore, implicit schemes are nearly infeasible and we need to re- sort to less favourable explicit schemes instead. In particular, we will use the parametrized second order one-leg method, which results in the momentum equation

(β +1

2)uⁿ⁺¹− 2βuⁿ+ (β− 1

2)uⁿ⁻¹+ ∆t Ω⁻¹Gpⁿ⁺¹ = ∆t F , (4.6) where F represents the convective and viscous forces evaluated at the previous timelevels:

F =−Ω⁻¹C((1 + β)uⁿ− βuⁿ⁻¹) + V ((1 + β)uⁿ− βuⁿ⁻¹). (4.7) Furthermore, the continuity equation is taken at the new timelevel:

Duⁿ⁺¹ = 0. (4.8)

It can be shown that the quantity u^n+1∗Ωuⁿ is conserved for this type of schemes. This quantity does not represent the energy, since it is not

(29)

a norm, but it comes close to it. The parameter β can be chosen freely and determines the stability region. For more stability analysis and the determination of the optimal β, we refer to [8].

Let ∆y_mindenote the minimum mesh size (which is somewhat suggestive, but in our simulations it indeed concerns the mesh size in the y-direction).

The timestep is then subject to the following stability conditions:

U ∆t < ∆ymin, (4.9)

which stems from the convective part of the momentum equation, and

2∆t < Re∆y_min² , (4.10)

which stems from the diffusive part. Although the latter inequality scales with ∆y_min² , Veldman [9] shows that (for properly chosen meshes) the convective limit is the most stringent. Therefore, the optimal timestep ∆t should be picked such that it satisfies both stability limits and scales with ∆y_min.

4.2.1 An equation for the pressure

The discretized Navier-Stokes Equations 4.6–4.8 clearly show which operators work on which timelevel. However, this formulation is not very practical yet, since both equations can never be solved simultaneously. Therefore, we rewrite them a bit to make it look like an algorithm. The first step is to update the velocity with convective and diffusive terms,

(β +1

2)˜u= 2βuⁿ− (β −1

2)uⁿ⁻¹+ ∆t F . (4.11) This yields a non-physical velocity field ˜u from which we can compute the new timelevel

uⁿ⁺¹= ˜u− ∆t Ω⁻¹Gpⁿ⁺¹. (4.12) In order to add this so-called pressure-correction term, the pressure pⁿ⁺¹ is required. This pressure can be solved from the Poisson equation which results from substituting 4.6 into 4.8. In our notation, it reads

DGpⁿ⁺¹ = ΩD ˜u. (4.13)

The Poisson equation is the most difficult equation of the three, and requires an iterative method to solve. For large grid sizes, by far the most computational effort is spent on this equation and it is therefore the limiting factor of most simulations (as we will see later on in Section 5.4).

(30)

Figure 4.3: Average number of iterations of the FGMRES Krylov solver, with and without pressure extrapolation. Taken from fully developed DNS Reτ= 1000.

4.2.2 Pressure extrapolation

When using an iterative method for the pressure-equation 4.13, it is very helpful to have a good initial approximation to the solution. An obvious choice is to take the pressure at the previous timelevel pⁿ as the initial solution, which is the preferred choice when there is no other historic information. However, during a typical simulation, one has access to multiple

‘old’ timelevels and can make the following extrapolation:

p₀= 3pⁿ− 3pⁿ⁻¹+ pⁿ⁻². (4.14) The idea behind this extrapolation is that, in general, Krylov solvers have difficulties eliminating the low frequency error, so to speak. The pressure field that is to be solved ’follows’ large scale structures (eddies) in the flow, that are moving with a certain velocity. This velocity is ‘measured’ by the extrapolation and eliminated from the initial solution. The Krylov solver should now need less iterations for solving the Poisson equation, which is indeed the case (Figure 4.3).

4.2.3 A remark

As mentioned before in Section 1.1, we deliberately do not use any spectral field representation of the Poisson equation. Since a channel flow in a rectan- gular box is statistically homogeneous in the mean-flow (x) and spanwise (z) directions, Fourier representation in those directions could speed up things significantly. However, our code is in theory applicable to more complex flows and geometries, even to those which are time dependent. These might not be statically homogeneous at all, in which case the spectral approach is useless.

(31)

Chapter 5 Results

This chapter contains the main results obtained during this research. First the general preparation of the channel and measurement methods is described, before showing the DNS in Section 5.2. Then, in Section 5.3 the performance of various PETSc solvers is discussed and in 5.4 the scalability of the code is investigated.

During the preparation of this thesis, a parallel running project was performed by Siebrich Kaastra dealing with Large Eddy Simulations. Her results have been obtained using the same codebase and the same hardware.

For all LES-related results and their discussion, the reader is referred to her thesis [3].

5.1 Preparation

A schematic view of geometry of the channel has been presented in Section 1.1. In our research we normalize against the height of the channel, i.e.

δ = 0.5 and h = 2δ = 1. The other lengths are calculated analogously as in Moser, Kim and Mansour [5], which yields Lx = π and Lz= π/2.

A pressure difference is imposed in the streamwise direction of the flow which should keep the bulk velocity constant U = 1. This is achieved by specifying the massflow parameter, which is the integral of the bulk velocity over the cross section of the channel.

All simulations are performed with 4th order discretization scheme (un- less otherwise specified) and take β = 0.05 in the parametrized timestepping method, see Section 4.2.

(32)

5.1.1 The nonuniform grid

The grid is chosen uniform in the streamwise (x) and spanwise (z) directions and stretched in the wall-normal (y) direction, according to

yj = sinh(γj/Nj)

2 sinh(γ/2). (5.1)

The stretching parameter γ is taken equal to 6.0.

Next, the number of gridpoints is calculated more or less by analogy. In order to solve all relevant length scales (up to the Kolmogorov scale) we need a fine enough mesh. It is shown that for one particular Reynolds number, the chosen mesh is sufficiently good. Subsequently, for higher Reynolds numbers we choose the grid s.t. the resolution in the y⁺ coordinate is kept constant.

The y⁺ coordinates are defined as

y_j⁺= yj · u^∗· Re, (5.2)

where y_j denotes the j-th gridline in wall-normal direction, u∗ the friction velocity or wall shear velocity and Re the Reynolds number based on the channel width and bulk velocity, Re = 2U δ/ν = 1/ν. We have found that an y₁⁺ value of about 1.2 yields good enough resolution to resolve all length scales, see for example Figure 5.3.

5.1.2 Friction velocity and Re_τ

An estimate of the friction velocity u∗ is still needed in order to compute y⁺ coordinates. Closely related is a quantity known as Re_τ. This (scaled) Reynolds number is common in engineering applications and is another way to classify wall-bounded flows based on the friction velocity and channel halfwidth alone. Defined as

Re_τ = u∗· δ

ν , (5.3)

it can be expressed in the global Reynolds number using Re = 1/ν:

Reτ = 1

2 · u∗· Re. (5.4)

Close to the wall, the shear stresses are quantified as τ_w = ρ· ν · (∂u

∂y)_y=0, (5.5)

which has dimension [kg/(s² · m)]. Now, using only variables nearby the wall, dimension analysis shows that the following should hold for the friction velocity,

u∗ = rτ_w

ρ . (5.6)

(33)

ReL Reτ Box size

Nom. Nom. Lx× h × L^z Ms Me

5,600 180 π× 1 × π/2 90 180 14,000 395 π× 1 × π/2 56 112 21,500 590 π× 1 × π/2 50 100 40,000 1000 π× 1 × π/2 23.8 50 56,000 1400 π× 1 × π/2 18 33.3 Table 5.1: Physical parameters for each DNS simulation.

The task of estimating u∗ is now equal to the task of estimating the shear stress τ_w (taking ρ = 1 for incompressible flow). Here is where experimental

‘hard’ data comes in. We use the skin friction coefficient proposed by Dean in [1],

C_f ≈ 2 · τw = 0.073Re⁻¹⁴. (5.7) This suffices to derive the optimal grid for any given Reynolds number.

Starting with a target value Reτ, the corresponding Re is computed using Dean’s skin friction coefficient and Equation 5.4. Subsequently, y⁺ coordinates are computed, thus the resolution is found and can be compared to the reference.

5.1.3 Spatial and time averaging

During the simulations, mean velocity fields and their statistics are collected and regularly written to disc for post processing purposes. These statistics have all been averaged over the channel streamwise and spanwise directions (i.e. the periodical directions) in order to reduce memory. As mentioned before in Section 1.1, the channel flow at hand is statically homogeneous in these directions, so we do not lose any information.

Furthermore, the statistics are averaged over time. After a start-up period of M_stime-units the flow is supposedly fully developed and the measurement is started. Statistics are averaged over the interval [Ms, Me] in time. The actual values of Ms and Me differ per simulation, depending on the Reynolds number, see Table 5.1. The heuristic used here is that a high Reynolds flow has better mixing properties and therefore needs less passes¹ before the flow is fully developed (Ms). Moreover, the measurement time window can be made shorter for this reason.

The aim was to integrate 200,000 timesteps in every simulation, of which 100,000 went into the start-up period and 100,000 into the measurement.

Multiplying these values by ∆t yields Ms and Me.

1One pass is defined as the time it takes for one ‘particle’ to traverse the channel once.

For a channel length of Lx= π and bulk speed unity, this takes about 3 time units.

(34)

ReL Reτ Reτ Grid size

Nom. Nom. Act. nx× n^y× n^z N ∆ymin ∆t y⁺₁

5,600 180 172 128× 128 × 128 2.1· 10⁶ 2.3· 10⁻³ 1.0· 10⁻³ 0.83 14,000 395 393 200× 200 × 200 8.0· 10⁶ 1.5· 10⁻³ 7.0· 10⁻⁴ 1.20 21,500 590 571 400× 300 × 400 4.8· 10⁷ 1.2· 10⁻³ 5.0· 10⁻⁴ 1,16 40,000 1000 1008 700× 500 × 700 2.5· 10⁸ 6.0· 10⁻⁴ 2.5· 10⁻⁴ 1.21 56,000 1400 1399 900× 700 × 900 5.7· 10⁸ 4.3· 10⁻⁴ 1.8· 10⁻⁴ 1.20 Table 5.2: Numerical parameters for each DNS simulation. The total number of

meshpoints is N = nx× n^y× n^z.

5.2 DNS results

Five direct simulations have been performed with the 4th order symmetry preserving discretization. Table 5.2 reveals which Reynolds numbers have been used and lists some of their properties. Calculated velocity profiles are to be compared with previous findings on the subject. In particular, we will use the following data-sets for comparison:

DNS results for Reτ = 180, 395, 590 channel flow by Moser, Kim and Mansour (referred to as MKM hereafter) from their 1999 article [5];

DNS results for Reτ = 934 channel flow by ´Alamo and Jim´enez from their 2003 article [2];

Experimental data for Reτ = 960 by Niederschulte, Adrian and Han- ratty;

Experimental data for Reτ = 1017, 1655 by Wei and Willmarth;

DNS results for Reτ = 1451 channel flow by Hu, Morfey and Sandham (referred to as HMS) from their 2006 article [10].

We note that all three DNS data-sets have been obtained by using the spectral method.

5.2.1 Discussion of velocity profiles

First we shall compare how well our results fit with existing data, which sources have just been mentioned. Based on the degree how well the data match we can say something about the reliability of the simulations. Mean velocity and root-mean square (RMS) profiles for Reτ = 180, the smallest Reynolds number of the lot, are shown in Figure 5.1. We observe that the profiles have good agreement with the corresponding external DNS results by MKM.

Velocity profiles for the second Reynolds number of Reτ = 395 are displayed in Figure 5.2. Again the mean velocity shows a good agreement with

(35)

10⁰ 10¹ 10² 0

5 10 15 20 25

U+

y⁺

(a) Mean velocity

0 20 40 60 80 100 120 140 160 180

0 0.5 1 1.5 2 2.5 3

u‘+,v‘+,w‘+

y⁺

(b) Rms

Figure 5.1: Reτ = 180, mean ux velocity and rms profiles. Dashed lines show corresponding DNS data by MKM. The upper, middle and lower profiles in the RMS chart represent u⁰⁺, the streamwise rms fluctuations;

w⁰⁺, the spanwise fluctuations; and v⁰⁺the wall-normal fluctuations, respectively.

MKM. However, our simulation seems to have difficulties in obtaining the correct maximum value for the streamwise rms fluctuations u⁰⁺. This is the only Reynolds number where this discrepancy has been observed, so we consider it to be an anomaly.

The largest nominal Reynolds number still available by MKM is Reτ = 590, which can be seen in Figure 5.3. The actual value turned out to be slightly lower in our simulation, a measured Re_τ = 571 (about 3% lower).

Again both the mean velocity and rms fluctuation profiles show good agreement with existing data.

There are multiple ways to investigate the quality of the statistics and therefore the simulation. For instance, one way is to look at the channel halves separately to see how much they differ. Until now we have averaged over both halves so any discrepancy will be hidden in the mean flow. In the long run however, we want the statistics of both halves to converge to the same mean profile. In Figure 5.4(a) we see the rms profiles for both our simulation and MKM in detail (mean velocity is not shown: they are virtually identical already).

Now then, the spread of two graphs is a measure of how well the simulation has converged in time. We see that the difference for urms(top graphs) is larger in our simulation but only slightly and in a certain range around y⁺ = 300. The wall-normal rms component has better convergence in our simulation (bottom graphs) but the spreads are both relatively small. For the spanwise wrms (center graphs) we observe again a larger spread in our own simulation.

In short we can conclude that, using this quality measure. the two

(36)

10⁰ 10¹ 10² 0

5 10 15 20 25

U+

y⁺

(a) Mean velocity

0 50 100 150 200 250 300 350 400

0 0.5 1 1.5 2 2.5 3

u‘+,v‘+,w‘+

y⁺

(b) Rms

Figure 5.2: Reτ = 395, mean ux velocity and rms profiles. Dashed lines show corresponding DNS data by MKM.

10⁰ 10¹ 10²

0 5 10 15 20 25

U+

y⁺

(a) Mean velocity

0 100 200 300 400 500 600

0 0.5 1 1.5 2 2.5 3

u‘+,v‘+,w‘+

y⁺

(b) Rms

Figure 5.3: Reτ = 590, mean ux velocity and rms profiles. Dashed lines show corresponding DNS data by MKM.

(37)

0 100 200 300 400 500 600 0

0.5 1 1.5 2 2.5 3

urms,vrms,wrms

y⁺

(a) Rms Reτ = 590

0 100 200 300 400 500 600 700 800 900 1000 0

0.5 1 1.5 2 2.5 3

urms,vrms

y⁺

(b) Rms Reτ = 1000

Figure 5.4: Rms velocity profiles where upper and lower channel halves are shown separately. Dashed lines show corresponding DNS data by MKM (left) and ´Alamo & Jim´enez (right).

simulations have comparable quality where perhaps the MKM simulation is the better of the two, by a narrow margin. Of course, this measure is not taking into account general convergence in time, for that affect both halves simultaneously.

For the next simulation with nominal value Reτ = 1000, we are comparing data-sets from a number of sources. Up to y⁺ = 100, the mean velocity profile follows the DNS by Jim´enez and the experimental data, as is shown in Figure 5.5(a). Between y⁺ = 100 and y⁺ = 400 the mean is slightly higher than the others, whereafter it is settling for a steady decrease in slope. At the centerline (y⁺ = 1000) this results in a substantial lower mean velocity than the external data. Basically what we observe here is the result of the flow not being fully turbulent. Since the initial condition at the centerline is u_c= 1 and during the start-up period we expect this value to go up, it is still not where it should be. The effect at the centerline is the most pronounced since it is furthest away from the geometry.

The rms profile might suffer from the same shortcoming, see Figure 5.5(b). Experimentation with the time-interval over which is averaged reveals that the rms profile is expected to go further down near the centerline, as the flow becomes more settled (not shown). This would result in an almost perfect fit with the experimental data from Wei & Willmarth at Reτ = 1017.

Furthermore we observe that the other experimental data (by Nieder- schulte) has problems getting the near wall region y⁺ < 20 right. This de- ficiency could (in theory) give rise to the incorrect calculation of u∗, which would mean the entire scaling of the graph is wrong. On this basis, we are hesitant about drawing conclusions which of the two simulations is the better one.

channel flow in the PETSc toolkit

Simulation of wall- bounded turbulent