Evaluation of Turbulence Models for Prediction of Flow Separation at a Smooth Surface

(1)

Project 9

Evaluation of Turbulence Models for Prediction of Flow Separation at a Smooth Surface

Eric Furbo, Janne Harju & Henric Nilsson

Report in Scientic Computing Advanced Course June 2009

PROJECT REPORT

(2)

(3)

Abstract

This project was provided by the Swedish Defence Research Agency, FOI. The task is to test several RANS (Reynolds-averaged Navier-Stokes) models on a generic case, constructed for flow separation at a smooth surface. This is of very high importance for applications in the auto industry, ship construction and for aeronautics. The models tested are implemented in the open source CFD (Computational Fluid Dynamics)-program, OpenFOAM.

OpenFOAM uses the finite volume method and the simple algorithm as solution procedure.

The main flow features evaluated in this project is the shape, position and size of the flow separation. Most of the models that we have tested have problems describing the complex dynamics of flow separation in this particular case.

(4)

List of Figures

1 Computational domain . . . 4

2 Inlet velocity profile . . . 5

3 The computational grid . . . 6

4 Case directory in OpenFOAM. . . 12

5 One type of inlet boundary conditions for k and ² . . . . 14

6 Results from the k − ² model. . . . 17

7 Results from the k − ². . . . 17

8 Results from the Nonlinear k − ² Shih model. . . . 18

9 Results from the Nonlinear k − ² Shih model. . . . 18

10 Results from the RNG k − ² model. . . . 19

11 Results from the RNG k − ² model. . . . 19

12 Results from the k − ω SST model. . . . 20

13 Results from the k − ω SST model. . . . 20

14 Results from the realizable k − ² model. . . . 21

15 Results from the realizable k − ² model. . . . 21

16 Results from the Spalart Allmaras model. . . 22

17 Results from the Spalart Allmaras model. . . 22

18 Comparison with LDV experiment. . . 23

21 Comparison with LDV experiment and LES data. . . 26

22 Comparison with LDV experiment and LES data. . . 27

List of Tables

1 Inlet conditions for k and ². . . . 14

2 Table of working turbulence models . . . 15

3 Table of turbulence models that encountered problems. . . 16

4 Table of separation and re-attachment points. . . 23

(6)

1 Introduction

The computational prediction of flow separation of a turbulent boundary layer from a smoothed surface is a process of primary concern. The fluid flow becomes detached from the surface and instead takes the form of eddies and vortices. Despite the simple shape (see figure 1) of the hill of which the flow separation occurs, it is a challenging computational problem and experiments indicate that the flow in the recirculation zone is complex and strongly time-dependent.

In the simulations described in this report we apply turbulence models based on RANS (Reynolds-average Navier-Stokes) equations. The domain that the flow have been simulated in is constructed to be very smooth without edges and has no real life application. However, the physical phenomena that arise is a subject of interest for many engineering components and systems. Streamlined car bodies, low-pressure turbine blades and highly loaded aircraft wings are some examples of where flow separation can have significant influence of the ability of the device in question to perform effectively.

The main objective is to evaluate the tested RANS-models prediction of the flow over the hill. We mainly evaluate the results with respect to the location and shape of the separation zone and the velocity distribution in the near wake of the wall. Our results will be compared to results from experiments made by Byun and Simpson [9] and LES simulations made by Persson et al. [7]

1.1 Problem Formulation

In order to simulate the turbulent flow and separation we have used the open source CFD- program OpenFOAM [3]. We have only used RANS turbulence models and OpenFOAM has a variety of such models implemented. Since all of the 14 RANS turbulence models have been tested, no extensive work have been done to models that in some way encountered problems to converge.

The geometry is a three-dimensional axisymmetric hill placed on the floor of a wind tunnel.

The computational domain is given in figure 1.

Figure 1: Computational domain seen from the inlet boundary.

(7)

The shape of the hill is defined by y(r)

H = − 1

6.04844 h

J0(Λ)I0

³ Λr

a

´

− I0(Λ)J0

³ Λr

a

´i

(1) where Λ = 3.1926, a = 2H. H = 0.078 m, is the height of the hill,a is the radius of the circular base of the hill and r² = x²+ z². J0 and I0 are the Bessel function of the first kind and the modified Bessel function of the first kind, respectively. A Cartesian coordinate system is used for the simulations, and in the domain coordinates (x, y, z) in this report also denoted by (x1, x2, x3) the velocity field is U = (ux, uy, uz) also denoted by U = (u1, u2, u3). The x-axis is parallel to the centerline of the domain and pointing downstream. The y-axis is normal to the floor, pointing upward. The origin is fixed at the center of the bump, with y = 0 corresponding to the wind tunnel floor.

The top and bottom patch seen in figure 1 are set to be wall boundaries. This implies the use of no slip conditions, i.e. U = 0 at these boundaries. The left, right and outlet patches are used to limit the computational domain. Therefore we will use homogeneous Neumann boundary conditions for these boundaries, hence ^∂U_∂n = 0.

The inlet velocity has a maximum flow velocity of 27.5 m/s. The inlet velocity is not uniformly distributed, it is instead a profile with lower velocities near the top and bottom wall and higher in the middle of the inlet patch corresponding to turbulent boundary layers of approximate thickness of H/2. The inlet velocity has only one non-zero component, that is in the x-direction.

The complete inlet profile can be seen in figure 2.

Figure 2: Velocity profile for the inlet patch of the computational domain The Reynolds number is a dimensionless quantity, defined as

Re = U L ν ,

and it represents the relation between inertial and viscous forces. Here U is a characteristic velocity, L is a characteristic length and ν is the kinematic viscosity. The Reynolds number for

(8)

our case, based on the maximum inlet velocity of 27.5 m/s, the height of the hill, H = 0.078 m and the kinematic viscosity of air ν = 1.65 · 10⁻⁵ m²/s is Re = 1.3 · 10⁵.

As for computational mesh we where given a mesh with about one million computational cells from FOI [7]. In the simulations made by FOI, a mesh with about 4 million cells was primary used. The mesh used in our simulations have the same distribution of cells as that one. The mesh is a structured mesh that have been refined near the wall boundaries to better resolve the boundary layers and the flow near this regions. The number of cells in x, y and z-direction is 69 × 119 × 119 respectively. The y+ value is related to the distance between a boundary and the nearest computational point. The y+ value for our mesh is in the range between 1.44 and 2.25, depending on which patch is considered. A wireframe plot of the mesh can be seen in figure 3.

On this computational domain the incompressible Reynolds-averaged Navier-Stokes (RANS) equations have been solved together with the different turbulence model. The Navier-Stokes equations will be briefly discussed in section 2.1, the RANS equations and turbulence models will be described in section 2.3.3. All the computations have been carried out in OpenFOAM.

OpenFOAM uses the Finite Volume Method and we have chosen the SIMPLE algorithm as solution procedure. This will be discussed in section 2.2.1 and 2.2.2 respectively.

(a) (b)

(c)

Figure 3: (a) The computational grid showing the cells on the surface of the hill and in the plane z = 0. (b) Zoomed in at the hill, notice the refined mesh near the hill boundary. (c) Mesh seen from the side.

(9)

1.2 Previous Research

This particular case has been extensively investigated by a number of computational research groups such as Garcia-Villalba et al. [4], Krajnovi´c [2] and Person et al. [7]. Since most of the previous research and the references therein for this case have relied on other turbulence models than RANS, such as DES (Detached Eddy Simulation), LES (Large Eddy Simulation) and hybrid LES-RANS, our objective is to compare how the different implemented RANS models in OpenFOAM performs compared to the other more computational costly models. In comparison we also have LDV (Laser doppler velocity) measurements for this case conducted by Byun and Simpson [9].

2 Theory

2.1 The Navier-Stokes equations

The fundamental basis for fluid dynamics are the Navier-Stokes equations. The fluid is assumed to be described as a continuum and the governing equations of the physics of fluid dynamics are that mass, momentum and energy are conserved. Derivation begins with application of the conservation of momentum, using control volumes and Reynolds transport theorem. The incompressible form of N-S equations and the conservation of mass is described as

ρ(dU

dt + U · ∇U) = −∇p + µ∇²U (2)

∇ · U = 0 (3)

where U is the flow velocity, p is the pressure, µ viscosity and ρ is the density. See Pope [8] for a complete derivation of the Navier-Stokes equations.

(10)

2.2 Numerical solution method

2.2.1 The Finite Volume Method (FVM)

The solution domain is subdivided into a finite number of small control volumes (polyhedra) and the conservation equations are applied to each control volume. FVM uses the integral form of a general convection-diffusion equation for a quantity as its starting point which is described as

d dt

Z

V

φdV + Z

S

φu · ndS = Z

S

Γ(∇φ) · ndS + Z

V

qdV (4)

where V is a control volume, S its bounding surface, n unit normal, Γ diffusion coefficient and q some external term. The i:th component of the momentum equation is of the same form as equation 4 with φ = ui.

The centroid of each control volume is assigned to be the computational node at which the variable values are to be calculated. Interpolation is used to express variable values at the control volume surface in terms of the centroid nodal values. Surface and volume integrals are approxi- mated using suitable discretization methods. As a result, one obtains an algebraic equation for each control volume, in which a number of neighbour nodal values appear. One advantage with FVM is that it can handle complex geometries, however one disadvantage compared to other computational methods is that methods of higher order than 2nd, are more difficult to develop in 3D. A more detailed description of the Finite Volume Method can be found in chapter 4 of [1].

2.2.2 The solution algorithm SIMPLE

When the RANS approach for turbulence is used, a stationary problem arises. RANS and many other methods for steady problems in computational fluid dynamics can be regarded as unsteady problems until a steady state is reached. If an implicit method is used in time, the discretized momentum equations at the new time step are non-linear. Due to this and that the underlying differential equations are coupled, the equations system resulting from discretization cannot be solved directly. Iterative solution methods are the only choice.

The momentum equations are usually solved sequentially for each component. The pressure used in each iteration is obtained from the previous time step and therefore the computed velocities normally do not satisfy the discrete continuity equation. In order for the velocities to fulfill this equation one have to modify the pressure field. This can be done by solving a discrete Poisson equation for the pressure.

After solving this new equation for the pressure the final velocity field at the new iteration is calculated. This new velocity field satisfies the continuity equation, but the velocity and pressure fields do not satisfy the momentum equations. Therefore, the procedure described above is iterated until a velocity field is obtained that satisfy both the momentum and continuity equations.

Methods of this kind which first construct velocity fields that do not satisfy the continuity equation and then correct them are known as projection methods. The SIMPLE algorithm is such a method, and it is the solving procedure used in OpenFOAM for our computations.

2.3 Turbulent Flows

Laminar flow changes into turbulent flow when the Reynolds Number becomes larger than a critical value. This value is dependent of the flow, thus the geometry of the domain. As mentioned in section 1.1 the Reynolds number for our case is Re = 1.3 · 10⁵. Since this indicate that turbulence will occur we have to account for the turbulence in some way.

(11)

Turbulent flows are highly unsteady. The direction of the velocity field fluctuates rapidly in all three spatial dimensions. These flows (are also vortical) and fluctuate on a broad range of time and length scales. Therefore it can be very hard to correctly simulate these types of flows.

2.3.1 Direct numerical simulation (DNS)

The most natural and straight forward approach to turbulence simulations is to solve the Navier- Stokes equations without any approximation of the turbulence other than numerical discretiza- tions. These type of simulations are called Direct Numerical Simulation (DNS). In such simula- tions all the motions contained in the flows are resolved. In order to capture the structures of the turbulence one have to use extremely fine grids and time steps for these type of simulations.

Due to this DNS is very costly from a computational point of view and this method is only useful for flows with low Reynolds numbers. The applications of DNS is therefore limited to turbulence research and results from DNS simulations can be important to verify results from other turbulence models.

2.3.2 Large Eddy Simulation (LES)

A simulation that is less computationally costly is Large Eddy Simulation (LES). This type of simulation uses the fact that the large scale motions are generally more energetic than the small scale ones. These large eddies are more effective transporters of the conserved properties and therefore LES threats the large eddies more exactly than the small scale ones. LES is preferred over DNS when the Reynolds number is too high or the computational domain is too complicated.

2.3.3 Reynolds-averaged Navier-Stokes models (RANS)

Both DNS and LES-modelling give a detailed description of turbulent flows. If one is not inter- ested in that much details of the turbulence, a Reynolds-Averaged Navier-Stokes model (RANS- model) can be a good choice. This type of model is much less costly compared to both DNS and LES. In the RANS approach to turbulence, all of the unsteadiness in the flow is averaged out and regarded as part of the turbulence. When averaging, the non-linearity of the Navier-Stokes equations gives rise to terms that must be modeled. These terms can be modeled differently leading to various types of RANS models, each behaving differently depending on the Reynolds number of the flow, the computational domain etc.

As mentioned above, the unsteadiness of the flow is averaged out. The flow variable, in this example it is one component of the velocity, is divided into two terms, one time-averaged value and the fluctuation about that value:

ui(xi, t) = ui(xi) + u⁰_i(xi, t) where

ui(xi) = lim

T →∞

1 T

Z _T

0

ui(xi, t)dt

xi are the spatial coordinates, also known as x, y and z. T is the averaging interval and must be large compared to the typical time scale of the fluctuations.

If the flow is unsteady, time averaging cannot be used and it has to be replaced with ensemble averaging. The concept of this is to imagine a set of flows in which all the variables that can be

(12)

controlled (energy, boundary conditions etc.) are identical but the initial conditions are gener- ated randomly. This will give flows that differ considerably from the others. An average over a large set of such flows is an ensemble average. In mathematical form written as

ui(xi, t) = 1 N

XN n=1

uni(xi, t) (5)

where N is the number of members of the ensemble. The term Reynolds averaging refers to any of the processes above and applying it to all the variables in the Navier-Stokes equations yields the Reynolds-averaged Navier-Stokes (RANS) equations which are written as

∂(ρui)

∂t + ∂

∂xj

³

ρuiuj+ ρu⁰_iu⁰_j´

= − ∂p

∂xi + µ ∂

∂xj

µ∂ui

∂xj +∂uj

∂xi

¶

(6)

∂ui

∂xi = 0 (7)

In the equations above, ρ is the density of the fluid and µ is the dynamic viscosity of the fluid.

The equation for the mean of a scalar quantity, temperature for example, can be written as

∂(ρφ)

∂t + ∂

∂xj

³

ρujφ + ρu⁰_jφ⁰

´

= ∂

∂xj

µ Γ∂φ

∂xj

¶

. (8)

As mentioned before the RANS equations give rise to a stationary problem. The time derivatives in equation 6 and 8 are used to advance the equations in an ”imaginary” time to reach steady state.

The RANS equations contain terms such as ρu⁰_iu⁰_jcalled the Reynolds stresses and u⁰_jφ⁰which is the turbulent scalar flux. The presence of these new quantities in the conservation equations means that they contain more variables than there are equations, The system of equations is not closed. Therefore one have to introduce approximations for the Reynolds stresses and the turbulent scalar flux. The approximations introduced are called turbulence models. To model the Reynolds stress the eddy-viscosity model is used, which is written as

ρu⁰_iu⁰_j = −µ_t µ∂ui

∂xj

+∂uj

∂xi

¶ +2

3ρδ_ijk (9)

where the eddy-diffusion model is used for the turbulent scalar flux

ρu⁰_jφ⁰= Γt∂φ

∂x_j. (10)

The turbulent kinetic energy k is described by k = 1

2

¡u⁰_xu⁰_x+ u⁰_yu⁰_y+ u⁰_zu⁰_z¢

(11) and µtis the eddy viscosity and will be defined later.

The exact equation for the turbulent kinetic energy k, is given by

∂(ρk)

∂t +∂(ρujk)

∂xj = ∂

∂xj

µ µ∂k

∂xj

¶

− ∂

∂xj

³ ρ

2u⁰_ju⁰_iu⁰_i+ p⁰u⁰_j

´

− ρu⁰_iu⁰_j∂ui

∂xj − µ∂u⁰_i

∂xk

∂u⁰_i

∂xk. (12)

(13)

The terms on the left-hand side of this equation and the first term on the right-hand side need no modelling. The second term on the right-hand side represents the turbulent diffusion of kinetic energy and is modelled according to

ρ

2u⁰_ju⁰_iu⁰_i≈ µt

σk

∂k

∂xj

(13) where σk is the turbulent Prandtl number whose value is approximately 1. The third term in equation 12 represents the rate of production of turbulent kinetic energy by the mean flow. Using the eddy-viscosity model for the Reynolds stress, this term is modelled to

−ρu⁰_iu⁰_j∂ui

∂xj ≈ µt

µ∂ui

∂xj +∂uj

∂xi

¶∂ui

∂xj. (14)

The equation for the dissipation of the kinetic energy, ², in it most commonly used form, is written as

∂(ρ²)

∂t +∂(ρuj²)

∂xj

= C²1Pk²

k − ρC²2²² k + ∂

∂xj

µµt

σ²

∂²

∂xj

¶

. (15)

In this model, the eddy viscosity is expressed as

µt= ρCµk²

² . (16)

Equation 15 for epsilon together with equation 12 is called the k−² model and is a commonly used RANS-model for turbulent flows. This model contains five parameters and the most commonly used values are

Cµ= 0.09 C²1= 1.44 C²2= 1.92 σk= 1.0 σ²= 1.3 (17) The most important difference in using a turbulence model like the k − ² model compared to just using the ordinary laminar equations is that two new equations have to be solved. One also have to replace the ordinary Navier-Stokes equations with the RANS-equations derived earlier.

One other thing that should be mentioned is that because the time scales associated with turbulence are much shorter than those connected to the mean flow. Therefore the equations connected to any turbulence model are much stiffer than the laminar equations. This has to be taken in to consideration in the numerical solution procedure.

Near the physical walls of a computational domain the profiles of the turbulent kinetic energy and its dissipation show a more significant peak than the mean velocity profile. One way to resolve the turbulent quantities near the wall is to use a finer grid. This will of course generate more computational points and therefore be more computational costly. One can also use so called wall functions for the velocity profile near walls. When these wall functions are used for the velocity, the normal derivative of k at the walls is assumed to be zero, i.e.

∂k

∂n = 0. (18)

At inflow boundaries it is quite hard to know what values to set for k and ². One common way for k at inlets are to set some small enough value times the square of the mean inlet velocity of the flow. To set the proper condition for ² one can use an approximation, written as

² ≈ k^3/2

L (19)

(14)

where L is the characteristic length off the flow. There are several other two-equation models that try to describe the turbulence. Many of them are based on the k − ² model described above.

Such models are the k − ω model, the RNG k − ² model and the Nonlinear k − ² Shih model.

They are in general more complicated and contain more parameters than the k − ² model.

There are also one-equation models like the Spalart Allmaras model. This model solves a transport equation for a viscosity-like variable ˜ν.

The last family of models that can be used to model the turbulence and are used together with the RANS equations are the Reynolds stress models. The assumption about the eddy- viscosity according to equation (9) is often a valid approximation for two dimensional problems.

For three-dimensional flows this approximation is not always accurate enough, since the eddy- viscosity may become a tensor quantity. One way to model this is to have a dynamic equation for the Reynolds stress tensor itself. These equations can be derived from the Navier-Stokes equations.

3 Simulation Parameters

3.1 OpenFOAM implementation

Open Field Operation and Manipulation (OpenFOAM) [6] Computational fluidmechanics toolbox is an open source program that can simulate a variety of physics, mainly fluid flows. The program is based on finite volume methods and is an object oriented C++ program consisting of modules that can be used to create a desired setup for a solver. The program also comes with tools for both pre- and post-processing of the solution.

3.1.1 Input data for a computational case

To use OpenFOAM a case must be set up, the case is basically a file structure that contains all the data that the program needs to process the problem. The picture below shows the file structure of an OpenFOAM case. To begin with, the case must contain mesh-files that the program can process to obtain a domain for the solution. Further more, the case must also contain files that specifies boundary conditions, solvers, time to start, stop and step, models to use and some more case specific data. When building the case the user has great freedom of choice regarding what solvers, preconditioners etc to use.

Figure 4: Case directory in OpenFOAM.

The three directories that the case consists of are: system, constant and time directories.

The system directory is for setting solution procedure parameters and must contain at least

(15)

three files. The controlDict, fvScheme and fvSolution. The constant directory contains the mesh information and specified physical properties of the domain. The time directories contain the initial values and further the solution over the domain.(for all the time steps that the user has specified to be written to file.) [5]

3.1.2 Our case

The parameters, preconditioners and solvers used in the simulations are specified in some files that are shortly be presented file by file below. Detailed examples of the structure of these files can be found in appendix A.

controlDict: In controlDict the solver is specified and in all of the simulations presented, sim- pleFoam is used as the solver. SimpleFoam is a steady state solver for incompressible, turbulent flow of non-Newtonian fluids, and it uses the solution algorithm SIMPLE described earlier in section 2.2.2.

fvSolution: To begin with, in the fvSolution file all solvers for the different parts of the equation have been specified. For solving the pressure, the preconditioned conjugate gradient method is used with diagonal incomplete Cholesky preconditioner. For the velocity and turbulence terms the preconditioned biconjugate gradient method is used, where diagonal incomplete LU- factorization is the preconditioner. Relaxation factors determine the convergence rate and are set for all quantities. The solution from the previous time iteration is multiplied with the relaxation factor before being employed in the current iteration. Lower relaxation factors are more likely to grant convergence but lead to longer computational time. Higher factor could speed up the process, but might also lead to divergence.

fvSchemes: In the fvSchemes files the schemes for discretizing different types of terms are specified, and for the _dt^d terms steady state solver is used. For gradient terms Gauss linear is being used, divergence terms are being discretized using Gauss linear or Gauss upwind and laplacian terms are being discretized with Gauss linear corrected and interpolation linear scheme.

RASproperties: RASproperties file is for specification of the turbulence model and the parameters included in the specified model. Also wall function coefficients are specified in this file.

3.2 Boundary Conditions for the turbulent quantities

All of the turbulence models that we have used except two (Spalart Allmaras and k −ω SST), use the turbulent kinetic energy, k, and its dissipation, ², to describe the turbulence. The boundary conditions for k and ² are set to a constant value on the inlet patch. For k the value is 0.756 m²/s²where we have used 27.5 m/s as reference velocity and computed the boundary condition according to section 2.3.3. The value for ² is set to 1.08 m²/s³ at the inlet, which is computed from the k inlet value according to equation 19. We used one tenth of the width of the domain as the characteristic length to compute the epsilon inlet value. This setup is the type 1 setup for the inlet conditions. For all the other boundary patches, Neumann boundary conditions were used.

To investigate the sensitivity of the flow to inlet turbulence, we also tried this model with four times higher and four times lower values of k, and the value of ² according to equation 19

(16)

Types of inlet conditions for k and ² Turbulent variable type 1 type 2 type 3 type 4

k [m²/s²] 0.756 3.024 0.189 non uniform profile, see figure 5 a

² [m²/s³] 1.08 8.64 0.135 non uniform profile, see figure 5 b Table 1: Table of the four types of inlet condition that have been used.

resulting from this new k value. These are the type 2 and type 3 setups and were only used for the k − ² model.

A fourth setup for the inlet values of k and ² that was not constant and hopefully more physical when modelling the turbulent boundary layer upstream of the hill, was also tried. To achieve this we ran a k − ² simulation in the channel flow environment without the hill. After convergence to steady state we used a plane parallel to the inlet, but inside the domain, with cell values of k and epsilon. The distribution of k and ² for this plane can be seen in figure 5.

These inlet conditions were used for the k − ² model, the Non linear k − ² Shih model and the realizable k − ² model. This setup is known as the type 4 inlet condition. For an overview of the types of inlet conditions, see table 1.

As for the k − ω SST model it uses identical boundary conditions for k as the k − ² model.

The inlet boundary condition for ω is simply ω =^k_².

The Spalart Allmaras model, which is a one equation model as mentioned in section 2.3.3, has only one turbulent quantity ˜ν. The boundary condition for ˜ν at the inlet boundary was set to 3 · 10⁻⁵.

0 0.05 0.1 0.15 0.2 0.25

0 1 2 3 4 5 6 7

y k [m²/s²]

(a)

0 0.05 0.1 0.15 0.2 0.25

0 2000 4000 6000 8000 10000 12000 14000 16000 18000

y epsilon [m²/s³]

(b)

Figure 5: Inlet boundary conditions according to type 4. (a) Distribution of k at inlet. (b) Distribution of ² at inlet.

(17)

4 Results

4.1 Short description of the physical experiment

Results from the LDV measurements by Byun and Simpson [9] states that there is no separation in front of the bump but that the flow decelerates there and then accelerate until the top of the bump. The mean flow on the lee side is closly symmetric around the centerline and complex vortical separation occurs downstream from the top and merge into large-scale turbulent eddies with two large streamwise vorticies. The flow along the streamwise centerline at x/H = 3.63 is a downwashing reattachment flow. The LDV experiment shows, with resulting velocity vectors in the plane of z/H = 0, that the mean location of separation is at x/H = 0.96. A more detailed description of the results from the LDV experiment and comparison with our results is discussed in 4.4.

4.2 Results from our simulations

All simulations of the different turbulence models have been conducted on the UPPMAX (Upp- sala Multidisciplinary Center for Advanced Computational Science) cluster called Os. A detailed list of the technical specifications of the Os cluster can be seen in appendix B. The wall clock time for performing a simulation on 4 cores for 10 time-steps was 3 min and 44 seconds.

The models that converge to a stationary solution are shown in table 2. Notice that even the simplest k − ² model needed almost 3000 time-steps ( about 18 hours 40 minutes of wall-clock time using 4 cores ) to converge. Keep in mind that the size of the time step is of no importance since a steady state solver was used.

Turbulence Model Description Comments

kEpsilon Standard k − ² model with wall functions

Convergence reached after approximately 3000 time-steps.

RNGkEpsilon RNG k − ² model with wall func- tions

Convergence reached after approximately 3000 time-steps.

NonlinearKEShih Non-linear Shih k − ² model with wall functions

Convergence reached after approximately 11 000 time-steps.

SpalartAllmaras Spalart-Allmaras 1-eqn mixing- length model for external flows

kOmegaSST k − ω SST two-equation eddy- viscosity model with wall functions

realizablekEps Realizable k − ² model with wall functions

convergence reached after approximately 17 000 time-steps.

Table 2: Table of working turbulence models

(18)

Turbulence Model Description Comments LienCubicKE Lien cubic k − ² model with wall

functions

Unable to reach convergence

QZeta q − ζ model without wall func-

tions

Solution did converge though the result was inconsistent with the physical solution.

LaunderSharmaKE Launder-Sharma low-Re k − ² model without wall functions

Unable to reach convergence after 30 000 time-steps (no difference in residuals for last 10 000 time-steps)

LamBremhorstKE Lam-Bremhorst low-Re k − ² model without wall functions

Solution starts to diverge after approximately 1000 time-steps LienLeschzinerLowRE Lien-Leschziner low-Re k − ²

model without wall functions

Solution start to diverge after approximately 800 time-steps

LRR Launder-Reece-Rodi Reynolds

stress model with wall functions

Solution starts to diverge after approximately 40 time-steps LaunderGibsonRSTM Launder-Gibson Reynolds stress

model with wall-reflection terms and wall functions

Solution starts to diverge after approximately 20 time-steps

Table 3: Table of turbulence models for which the simulation encountered problems, preventing it to produce a converged, steady solution.

The models that in any way failed to give accurate results are shown in table 3. Recall that the main objective was to try as many RANS turbulence models as possible for this particular case. The time frame and the priorities of the project excluded the possibility to more closely investigate models for which problems were encountered. It is quite possible that with an extra effort, the turbulence models in table 3 may work for this case as well.

The k − ² types of models that failed to converge were tested with more suitable under relaxation factors and different types of inlet boundary conditions without success. The Reynolds stress models that diverged early in the simulations were tested numerous ways. We tried running a k − ² method first, start the simlutation after 1000 time-steps with the Reynolds stress models without success. We also tried suitable relaxation factors which also did not make any difference either.

4.3 Global behaviour of the models tested

All of the models that converged capture the magnitude of the velocity quite reasonably compared to the physical solution. Separation of the flow is found in all models except the k − ² model.

Generally the separation region is too large compared to the physical solution, the start of the separation is however reasonably consistent with the physical solution. In all of the cases the motion of the turbulent eddies is incorrectly resolved, such as the flow direction and the magnitude of this flow. Near the centerline two large eddies appear which should not be there.

A more detailed description of the solutions computed with the models in table 2 is listed below and comparisons to the LDV measurements and LES data are found in section 4.4.

(19)

4.3.1 The k − ² model

The first model that we tried for our domain was the basic k − ² model. It is often used as first model to test with. The result from this model was not accurate, the simulated flow was not separated which is shown in figure 6 and 7 , and no eddies or vortices were captured.

The result from the simulation that used the inlet setup 2 and setup 3 did not differ substantially from the result with inlet conditions according to the first setup described in 3.2. When we used setup 4 for the inlet conditions we saw only very minor differences in the solution compared to the two other inlet setups for this model. Since the difference was none or very small with these conditions, no results in form of graphs will be presented.

(a) (b)

Figure 6: (a) The magnitude of the velocity at z = 0. (b) The velocity field in the wake on the leeward side of the hill. Observe that there is no recirculation or separation.

(a) (b)

Figure 7: (a) The turbulent kinetic energy k, at z = 0. Compared to other models the k − ² model produces a much more spread out profile of this quantity. (b) Streamlines of the velocity in x-direction through the line from [0.125 0.02 0.15] to [0.125 0.02 -0.15].

(20)

4.3.2 The Nonlinear k − ² Shih model

The results of this model differ from the results from other models. It was tested both with inlet conditions accordning to setup 1, the same as the rest of the models, and also with setup 4.

As seen in figure 9, both cases gave unrealistic results, with different asymmetric solutions and occurrence of only one large vortex either to the right or to the left of the center line downstream of the hill. The results in both cases can be considered bad compared to the physical solution and the LES-simulation. It is however an interesting result that this turbulence model ”supports”

an asymmetric solution (and it’s mirror on the centre plane, of course).

(a) (b)

Figure 8: (a) The magnitude of the velocity at z = 0. (b) The velocity field in the wake on the leeward side of the hill.

(a) (b)

Figure 9: The streamlines from this model are completely different from other models that were tested. (a) Streamlines of the velocity in x-direction through the line from [0.125 0.02 0.15] to [0.125 0.02 -0.15]. (b) Streamlines of the velocity in x-direction when using inlet setup 4. through the line from [0.125 0.02 0.15] to [0.125 0.02 -0.15].

(21)

4.3.3 The RNG k − ² model

This model was used with inlet conditions according to setup 1 and setup 4 but converged only for setup 1! As seen in figure 10 (a), the flow layer near the bottom of the domain decelerates before the hill and then accelerates on the top to finally decelerate again on the leeward side of the hill. The flow separates around x/H ≈ 0.96 and is re-attached around x/H ≈ 2.12. The separation and re-attachment can be seen in figure 10 (b). Two small and two bigger vortices develop downhill of the hill in the positive x-direction. This is shown in figure 11 (a). The turbulent kinetic energy is shown in figure 11 (b). As expected more turbulent kinetic energy appears in the separation zone where the eddies are as greatest.

(a) (b)

Figure 10: (a) The magnitude of the velocity at z = 0. (b) The velocity field in the wake on the leeward side of the hill.

(a) (b)

Figure 11: (a) Streamlines of the velocity in x-direction through the line from [0.125 0.02 0.15]

to [0.125 0.02 -0.15]. (b) The turbulent kinetic energy k, at z=0.

(22)

4.3.4 The k − ω SST model

This model was used with inlet conditions according to setup 1 and setup 4 but converged only for setup 1! The model managed to clearly capture the flow separation. The flow separation occurs around x/H ≈ 0.71 and is re-attached around x/H ≈ 2.24. This can be seen in figure 12 (a).Two vortices on the leeward side of the hill occur and they are almost symmetrically placed around the plane z = 0. This is shown in figure 13 (a).

(a) (b)

Figure 12: (a) The magnitude of the velocity at z=0. (b) The velocity field in the wake on the leeward side of the hill.

(a) (b)

(23)

4.3.5 The realizable k − ² model

This model was used with inlet conditions according to setup 1 and setup 4 and converged for both setups. The results from the two inlet setups were almost identical and only the results from setup 1 will be presented in the graphs. Similar to the k − ω SST model the realizable k − ² model also managed to capture the flow separation. The flow separation occurs at x/H ≈ 0.64 and is re-attached around x/H ≈ 2.24. This can be seen in figure 14 (a).Two vortices on the leeward side of the hill turn up and they are almost symmetrically placed around the plane z = 0.

This is shown in figure 15 (a).

(a) (b)

Figure 14: (a) The magnitude of the velocity at z=0. (b) The velocity field in the wake on the leeward side of the hill.

(a) (b)

(24)

4.3.6 The Spalart Allmaras model

This model was only used with inlet setup 1. The result of this model is a bit different compared to the other turbulence models that managed to capture the flow separation. In the plane z = 0, no separation is appearing and the wake is divided in two symmetrically placed separation

”bubbles”. This is shown in figure 16. In difference to all other models this model generates four vortices, two at each side of the centerplane, at the leeward side of the hill.

(a) (b)

Figure 16: (a) The velocity field in the wake on the lee side of the hill. Notice that no separation of the flow is found at this cut. (b) z=0.02

(a) (b)

to [0.125 0.02 -0.15]. Notice that this model captures 4 vortices at the lee side of the hill. (b) The magnitude of the velocity at z=0.

(25)

4.4 Comparison to experimental and LES data

Results of the LDV investigations of the flow around the hill where reported by Byun and Simpson 2005 [9]. Their result for the velocity field lee ward of the hill in the z/H = 0 plane is shown in figure 18. The flow separates at x/H = 0.96. Note that the k − ² model did not have any flow separation and the Non linear k − ² Shih model showed obvious strange behaviour and therefore these models have not been included in the comparison section. The Spalart Allmaras model shows no separation at all for the plane z/H = 0, however separation is found in the plane z/H = ±0.02 at x/H = 0.90.

Model Separation point Re-attachment point

LDV, figure 18 (a) x/H = 0.96 x/H ≈ 2.0

RNG k − ², figure 10 (b) x/H = 0.96 x/H = 2.12 Realizable k − ², figure 14 (b) x/H = 0.64 x/H = 2.24 k-Omega SST, figure 12 (b) x/H = 0.71 x/H = 2.24

Table 4: Table of separation and re-attachment points for different models and LDV experiment.

Note that no obvious re-attachment point is found in the LDV-experiment.

As we can see from the comparison in table 4 our results in the same plane show that the RNG k − ² model have similar separation region. The realizable k − ² and k − ω model separate the flow to early at x/H = 0.64 and x/H = 0.71 respectively.

(a) (b)

Figure 18: (a) Velocity field from LDV measurements in the plane z = 0. (b) Contours of the magnitude of the velocity U , normalized with reference velocity Uref = 27.5 m/s and the vector field of Uy and Uz in the plane x/H = 3.63 for the LDV measurements. Note that the view is from behind the hill. Plots made by Byun and Simpson [9]

If we look at the results from the LDV-experiment in the y, z-plane with x/H = 3.63, there are two vortices centred at z/H ≈ ±1.35. This can be seen in figure 18 b. All of our most successful models fail to predict this result, however some similar results can be found from a couple of our RANS models.

The Spalart Allmaras model produces two large vortices centred at z/H ≈ ±0.45, and two very small ones at z/H ≈ ±1.92. This is shown in figure 20 a.

(26)

(a)

(b)

Figure 19: (a) Surface plot of the magnitude of velocity U and the transverse vector field(in the y, z-plane), for the RNG k − ² model. (b) k − ω SST model. Note that the view is from behind the hill.

The k − ω SST model manage to produce two vortices at z/H ≈ ±1.41 which are quite similar to the LDV experiment. However, this model behaves strangely around z/H ≈ ±0 with two small vortices that collide. This is shown in figure 19 b.

The realizable k − ² model produces a result that is quite similar to the k − ω SST model.

There are two vortices at z/H ≈ ±1.54 and they are a bit smaller than the ones produced by the k − ω SST model at around the same place. This model, similarly to the k − ω SST model, fail by producing two vortices that collide around z/H ≈ 0. These vortices can be seen in figure 20 b.

The RNG k − ² model produces two vortices at z/H ≈ ±1.40, which is consistent to the LDV experiment. In difference with the LDV experiment some inaccurate behaviour occurs around z/H = 0.38. These vortices can be seen in figure 19 a.

(27)

Figure 20: Surface plot of the magnitude of velocity U and the vector field of Uy and Uz. Note that the view is from behind the hill. (a)The Spalart Allmaras model (b) The realizable k − ² model.

The line plots in figure 21 and figure 22 show the ux and uz respectively at x/H = 3.63 for different values of z/H. Note that these plots show the velocity components behind the separation region. Therefore ux is positive for all y/H-values in each plot. However, there are still vortices in this region, but they appear in the y, z-plane according to figure 19 and figure 20.

The results for the ux velocity from our RANS turbulence models are quite similar with LDV-experiment except some larger deviations for the Spalart Allmaras and the Non linear k − ² Shih models. The latter of these two models produces for almost all z/H points inconsistent results compared to the reference data.

If we look at the u_z velocity almost all of the RANS turbulence models produce different results compared to the LDV-experiment. As one can see in figure 22 for z/H = 0.8143 and z/H = 0.3257, the uz velocity points in opposite direction for small values of y/H. The flow direction of the eddies are opposite compared to the LDV and LES data. If we look in the centerplane, hence z/H = 0, we can see in figure 22 that both the realizable k − ² model and the k − ω SST model produce incorrect results, which have been mentioned before.

(28)

0 0.2 0.4 0.6 0.8 1 1.2 1.4 0

0.2 0.4 0.6 0.8 1 1.2 1.4

U_x/U₀

y/H

RANS−models compared with LDV measurements and LES data, z/H = 0.8143 LDV data

LES data kEpsilon SpalartAllmaras NonlinearKEShih RNGkEpsilon kOmegaSST realizablekEps

0 0.2 0.4 0.6 0.8 1 1.2 1.4

U_x/U₀

y/H

RANS−models compared with LDV measurements and LES data, z/H = 0.3257 LDV data

0 0.2 0.4 0.6 0.8 1 1.2 1.4

U_x/U₀

y/H

RANS−models compared with LDV measurements and LES data, z/H = 0 LDV data

0 0.2 0.4 0.6 0.8 1 1.2 1.4

U_x/U₀

y/H

RANS−models compared with LDV measurements and LES data, z/H = −0.3257 LDV data

0 0.2 0.4 0.6 0.8 1 1.2 1.4

U_x/U₀

y/H

0 0.2 0.4 0.6 0.8 1 1.2 1.4

U_x/U₀

y/H

Figure 21: Comparison of RANS models with LDV experiment and LES data. Ux for RANS and LDV is measured at x/H = 3.63 and at x/H = 3.69 for LES. Data is for six different values of z/H.

(29)

−0.10 −0.05 0 0.05 0.1 0.15 0.2

0.4 0.6 0.8 1 1.2 1.4

U_z/U₀

y/H

RANS−models compared with LDV measurements and LES data, z/H = 0.8143

LDV data LES data kEpsilon SpalartAllmaras NonlinearKEShih RNGkEpsilon kOmegaSST realizablekEps

−0.20 −0.1 0 0.1 0.2 0.3

0.2 0.4 0.6 0.8 1 1.2 1.4

U_z/U₀

y/H

RANS−models compared with LDV measurements and LES data, z/H = 0.3257

−0.250 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.2

0.4 0.6 0.8 1 1.2 1.4

U_z/U₀

y/H

RANS−models compared with LDV measurements and LES data, z/H = 0

−0.40 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.2

0.4 0.6 0.8 1 1.2 1.4

U_z/U₀

y/H

RANS−models compared with LDV measurements and LES data, z/H = −0.3257

−0.40 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.2

0.4 0.6 0.8 1 1.2 1.4

U_z/U₀

y/H

−0.150 −0.1 −0.05 0 0.05

0.2 0.4 0.6 0.8 1 1.2 1.4

U_z/U₀

y/H

Figure 22: Comparison of RANS models with LDV experiment and LES data. Ux for RANS and LDV is measured at x/H = 3.63 and at x/H = 3.69 for LES. Data is for six different values of z/H.

(30)

5 Conclusions

The results from the RANS models implemented in OpenFOAM of the flow over a three- dimensional axisymmetric hill have been reported, analysed and compared to experiment- and LES-data. Many other research groups like Wang et al. [10] and Person et al. [7] have tested RANS turbulence models on the same geometry. Their results are confirmed by our study, the RANS turbulence models tested for this geometry do not give accurate results.

Even though none of the tested models give good enough result some models do capture the separation region fairly well. Our conclusion is that the best performing model in these tests is the RNG k − ² model. The wake in the centerplane, z/H = 0 is most similar in size and position compared to the results from LDV experiment. In addition to the RNG k − ² model two more models give results that are acceptable. These are the realizable k − ² and the k − ω SST models.

All of the tested RANS models, even the three that produces acceptable results, produces eddies in the plane x/H = 3.63 near the centerline that do not exists in the experimental data.

Improvements may be possible if the meshes are refined, especially in the turbulent boundary layer near the hill. Another possible way to improve the results is to try other types of inlet boundary conditions for the turbulent quantities. Finally one could also further develop better adapted wall functions to handle the turbulent boundary layer close to the walls.

References

[1] Milovan Peri´c Joel H. Ferziger. Computational Methods for Fluid Dynamics. Springer, 3rd edition, 2002.

[2] Siniˇsa Krajnovi´c. Large eddy simulation of the flow over a three-dimensional hill. Flow Turbulence Combust, 81:189–204, December 2008.

[3] OpenCFD LTD. http://www.opencfd.co.uk/openfoam/.

[4] W. Rodi M. Garcia-Villalba, N. Li and M. A. Leschziner. Large-eddy simulation of separated flow over a three-dimensional axisymmetric hill. J. Fluid Mech., 000:1–42, 2009.

[5] OpenCFD Limited. OpenFOAM - Programmer’s Guide, 1.5 edition, July 2008.

[6] OpenCFD Limited. OpenFOAM - User Guide, 1.5 edition, July 2008.

[7] Benson R. E. Persson T., Liefvendahl M. and Fureby C. Numerical investigation of the flow over an axisymmetric hill using les, des and rans. Journal of Turbulence, 7(4):1–17, 2006.

[8] Stephen B Pope. Turbulent Flows. Cambridge University Press, 2000.

[9] G.Byun Roger L. Simpson. Structures of three-dimensional separated flow on an axisym- metric bump. 43rd AIAA Aerospace Sciences Meeting and Exhibit, 23, 2005. Paper No.

2005-0113.

[10] Jang Y. J. Wang C. and Leschziner M. A. Modelling two- and three-dimensional separation from curved surfaces with anisotropy-resolving turbulence closure. International Journal of Heat and Fluid Flow, 25(3):499–512, 2004.

(31)

A Examples of OpenFOAM case-files

File: controlDict

/*---*\

| ========= | |

| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |

| \\ / O peration | Version: 1.5 |

| \\ / A nd | Web: http://www.openfoam.org |

| \\/ M anipulation | |

\*---*/

FoamFile {

version 2.0;

format ascii;

class dictionary;

object controlDict;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application simpleFoam;

startFrom startTime;

startTime 0;

stopAt endTime;

endTime 1.0;

deltaT 0.0001;

writeControl timeStep;

writeInterval 1000;

purgeWrite 0;

writeFormat ascii;

writePrecision 6;

writeCompression uncompressed;

timeFormat general;

timePrecision 6;

runTimeModifiable yes;

// ************************************************************************* //

(32)

File: fvSchemes

/*---*- C++ -*---*\

| ========= | |

| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |

| \\ / O peration | Version: 1.5 |

| \\ / A nd | Web: http://www.OpenFOAM.org |

| \\/ M anipulation | |

\*---*/

FoamFile {

version 2.0;

format ascii;

class dictionary;

object fvSchemes;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes {

default steadyState;

}

gradSchemes {

default Gauss linear;

grad(p) Gauss linear;

grad(U) Gauss linear;

} divSchemes {

default none;

div(phi,U) Gauss upwind;

div(phi,k) Gauss upwind;

div(phi,omega) Gauss upwind;

div(phi,R) Gauss upwind;

div(R) Gauss linear;

div(phi,nuTilda) Gauss upwind;

div((nuEff*dev(grad(U).T()))) Gauss linear;

}

laplacianSchemes {

default none;

laplacian(nuEff,U) Gauss linear corrected;

laplacian((1|A(U)),p) Gauss linear corrected;

laplacian(DkEff,k) Gauss linear corrected;

laplacian(DomegaEff,omega) Gauss linear corrected;

laplacian(DREff,R) Gauss linear corrected;

laplacian(DnuTildaEff,nuTilda) Gauss linear corrected;

}

interpolationSchemes {

default linear;

interpolate(U) linear;

}

snGradSchemes {

default corrected;

}

fluxRequired {

default yes;

p;

}

// ************************************************************************* //

Evaluation of Turbulence Models for Prediction of Flow Separation at a Smooth Surface

Full text

Project 9