CFD model of fluid flow in reactor: A simulation of velocity and heat distribution in a channel

(1)

UPTEC ES 06 012

Examensarbete 20 p

September 2006

CFD model of fluid flow in reactor

A simulation of velocity and heat distribution

in a channel

(2)

(3)

Abstract

The basic problem of operating a boiling water nuclear reactor (BWR) is that of maximizing the power output while avoiding fuel rod over-heating (dry-out). For the safe operation of BWRs this entails a detailed understanding of the flow of water and steam through the reactor core. In a BWR the water functions not only as a coolant, but also as a moderator for the neutrons emitted in the fission process. To describe the thermohydraulic properties of the reactor a number of parameters are of common interest. Examples of such parameters are void, pressure, temperature, water and steam velocities, pressure, sheer forces and turbulent kinetic energy. There are a few ways of revealing these values such as experiments built up to behave like a reactors and computer simulations using models based on the laws of fluid dynamics and thermodynamics.

This research concerns a computer based model which uses a continuous fluid dynamic, (CFD), calculation program called OpenFOAM (Open Field Operation and Manipulation). OpenFOAM uses Navier Stokes equations for continuity, momentum- and energy conservation to simulate flows. The method used in this research has been to first build a model which describes the adiabatic flow correctly by using an already existing solver which uses the continuity and momentum conservation laws. In order to achieve a model that can solve temperature distributions in the flow the energy equation is added to the program coding in OpenFOAM.

There are totally 12 turbulence models. Some of the models have not produced results on account of that they either diverged or needed input that was not attainable. The models which were tested and used were four k - ε models, one RSTM model and two low-Re models.

A question that is addressed in this report is which of the many turbulence models that describes the experimental flow most correctly. The low-Re model LienLeschzinerLowRe produces the results with best congruence to the experimental data. The k - ε model RNGkEpsilon and the RSTM model LRR were also fairly close.

It is found that there are a few turbulence models that describe the experimental flows sufficiently well. LaunderSharmaKe was the turbulence model which simulated the temperature distribution best and was almost within the error bar limit of 5 K in all of the plots.

It is interesting that the two low-Re models show the best results if only one characteristic at the time is studied. One of the turbulence models describes only the velocity profile well and the other one oppositely describes the temperature distribution the best. It can thereby be stated that if the user wants a turbulence model that describes both velocity profiles and temperature distribution the RSTM model LRR is the best one. If on the other hand computer capacity is a limiting factor it might be profitable to use the simpler k - ε model RNGkEpsilon.

An other conclusions of this thesis is that the LRR and RNGkepsilon models are suited for the simulations of the geometries described in this work provided that the channel is wide enough for the model to simulate a correct temperature profile. With the use of gmsh a case geometry with wider channel area could easily be created. It would of course be necessary to use experimental data to validate the assumption that more realistic results can be obtained on wider channels.

(4)

(5)

Sammanfattning

Det grundläggande problemet med att driva en kokvatten reaktor (BWR) idag är att maximera effekten samtidigt som överhettande av bränslestavarna undviks (dry-out). För att kunna driva BWR:en på ett säkert sätt krävs det detaljerad förståelse över transporten av vatten och ånga genom reaktorkärnan. I en BWR fungerar inte vattnet bara som ett kylmedium utan även som en moderator för neutronerna som emitteras genom fissionsprocessen. För att beskriva de termohydraliska egenskaperna av reaktorn är vissa parametrar av speciellt intresse. Exempel på sådana parametrar är void, tryck, temperatur, vatten- och ånghastighet, skjuvkrafter och turbulent kinetisk energi. För att ta reda på värdena av dessa parametrar förekommer det vissa metoder som till exempel göra experiment som liknar egenskaperna i en reaktor eller datorbaserade simulationer som använder modeller som grundar sig på fluid dynamik och termodynamik.

Fokus på den här rapporten är inriktat på computational fluid dynamics (CFD) genom att använda program coden OpenFOAM (Open Field Operation and Manipulation) för att simulera transporten av vatten genom olika konfigurationer av bränsleknippen. Experiment med uppvärmd luftström utfört av Nuclear Engineering and Design har studerats och använts för validerings data för OpenFOAM modellen.

OpenFOAM används i en trestegs process. Först bygga en modell som beskriver det adiabatiska flödet. Detta görs genom att använda en redan befintlig solver som använder Navier-Stokes bevarande lagar för massa och moment. Därefter utvecklas modellen till att beskriva flödets parametrar när värme är applicerat till bränslestavarna. För att ta fram en modell som kan lösa temperatur distributioner i flödet måste energi ekvation implementeras i solverns programkod. Denna modell skall sedan användas för att estimera flödes parametrar när flödet passerar olika konfigurationer av bränslestavarna som liknar mer de bränsleknippen som används i dagens kärnkraftverk. På grund av att arbetsgången är uppdelad i tre etapper så har även rapporten delats upp på samma sätt genom att först hantera adiabatiska flöden sen uppvärmda flöden och till sist flöden genom 2x2 samt 3x3 geometrier.

Även om de osammanhängande fluktuationerna av flödes variabler är av deterministisk art så fortsätter försök att simulera turbulenta flöden att generera betydelsefulla problem.

Det är möjligt att beräkna turbulens genom att använda den tidsberoende versionen av Navier-Stokes equation genom en process som är känd som the Direct Numerical Simulations (DNS). Eftersom denna process erfordrar CPU-tid i storleks ordning av Re3 och en rumslig upplösning skala i storleks ordning av Re9/4 så är det nödvändigt att ta hänsyn till de turbulenta effekterna på ett approximativt sätt. Med avseende på detta så har en stor variation av turbulensmodeller tagits fram och fortsätter att utvecklas. En huvudsaklig fråga i den här rapporten är därför att avgöra vilken av de olika turbulens modellerna som beskriver den experimentella datan mest korrekt.

Syftet med detta arbete är att ta fram en CFD baserad modell som kan simulera de parametrar som är av stor betydelse för flödets karaktär i en reaktor.

Modellen i detta arbete är baserat på ett experiment där turbulent luft flöde har blåsts igenom en central kanal av ett el-uppvärmt bränsleknippe med 37 rör. Rören var konfigurerade i en triangulär uppbygnad med två olika pitch-diameter (P/D) kvoter. Rörens längd är 11.5 meter varav den första längden på 4.6 meter är ouppvärmd medan den resterande delen på 6.9 meter är uppvärmd. Uppvärmningen sker genom att ström leds genom en tunn metallfolie utanpå rören. Eftersom metallfolien har en yta som inte varierar mer än ±2 μm kommer inte foliens påverkan på fluidströmningen att tas hänsyn till i modellen.

Den huvudsakliga tekniken bakom CFD programmet OpenFOAM är dess flexibla set av C++ moduler. Dessa används för att bygga så kallade solvers för att simulera specifika problem

(6)

inom ingenjörskonst. OpenFOAM använder finita volym beräkningar för att lösa system av partiella differential ekvationer (PDE) tillhörande en godtycklig struktur av en mesh som består av månghörniga celler.

Det finns totalt 12 turbulens modeller men de turbulens modeller som har producerat resultat utan att divergera är av sorterna k - ε (4 modeller), Reynolds stress turbulence models (RSTM) (1 modell) och low-Re (2 modeller).

Resultaten från dessa är plottade tillsammans med den experimentella datan. Det visar sig att de visar sig att kEpsilon och RNGkEpsilon modellerna ger någorlunda samma resultat som RSTM modellen LRR men att LRR är bäst korrigerad med den experimentella datan. Den modellen som beskriver experimentella datan allra bäst är LienLeschzinerLowRE och tillhör low-Re turbulens modellerna.

Samma modeller har används för att beskriva temperaturfördelningen i kanalerna. Det visar sig svårt att estimera temperaturfördelningen och ingen av modellerna är särskilt bra beskrivande för den experimentella datan. Bedömningen har dock gjorts att LaunderSharmaKE modellen som tillhör low-Re har resultat som är närmst de experimentella värdena.

RNGkEpsilon modellen har använts för att beräkna hastighets profilen samt temperatur distributionen i de mer komplexa geometrierna, 2x2 och 3x3. Geometrierna har fler celler och fallen tar därför längre tid att räkna ut en de som använts för att bygga modellen.

Sammanfattningsvis så har LienLeschzinerLowRE bäst för temperatur fördelningen och LaunderSharmaKE bäst resultat för temperatur distributionen. Om dock en modell skall väljas för att modellera alla egenskaperna hos ett flöde så är LRR modellen bäst men kräver lite mer data resurser än RNGkEpsilon som är en enklare modell.

(7)

1 Introduction

1.1 Background

The basic problem of operating a boiling water nuclear reactor (BWR) is that of maximizing the power output while avoiding fuel rod over-heating (dry-out). For the safe operation of BWRs this entails a detailed understanding of the flow of water and steam through the reactor core. In a BWR the water functions not only as a coolant, but also as a moderator for the neutrons emitted in the fission process. To describe the thermohydraulic properties of the reactor a number of parameters are of common interest. Examples of such parameters are void, pressure, temperature, water and steam velocities, sheer forces and turbulent kinetic energy. There are a few ways of revealing these values such as the use of experiments built up to behave like the flows in a reactor core or computer simulations using models based on the laws of fluid dynamics and thermodynamics.

1.2 Problem definition

The focus of this report is on computational fluid dynamics (CFD) using the code OpenFOAM (Open Field Operation and Manipulation) [1] for the simulation of the flow of water through different fuel geometries. Previous performed experiments with heated air flow performed by the group Nuclear Engeneering and design [2] will be studied which will be used as validation data for the OpenFOAM model.

OpenFOAM is used in a three-step process. First a model is created that describes the adiabatic flow. Then a second model, which is developed from the first one, describes the properties of the flow when a heat flux is applied to the fuel rods. The model of the heated flow is then used to estimate the properties of the fluid passing through other channel-like geometries. The structure of this report is based on these steps, chapters 3, 4 and 5 dealing with unheated flow, heated flow and 2x2 and 3x3 geometries respectively.

Although the chaotic fluctuations of the flow variables are of deterministic nature, the simulation of turbulent flows still continues to present significant problems. It is possible to calculate turbulence using the time-dependent Navier-Stokes equation in a process known as the Direct Numerical Simulation (DNS). Because this simulation requires a CPU-time in the order of Re3 and spatial resolution scales as Re9/4 it is necessary to account for the effects of turbulence in an approximate manner. For this purpose, a large variety of turbulence models have been developed and the research still goes on. Some of these models will be used in an attempt to describe the behaviour of the flow.

1.2 Purpose

The purpose of this research is to find a model, which can properly simulate both adiabatic and a heated flows through a nuclear reactor core.

1.3 Delimitation

A delimitation of this work is the use of symmetry sides. By using symmetry sides as is done when building this model one actually simulates infinitely many fuel rods and channels. Another delimitation has been not to use temperature dependent variables for the energy equation but setting them fixed values to an estimated mean value of the fluid body.

(10)

2 The reference experiment

Abstract of experiment

In order to formulate a valid model empirical data are needed. In this thesis, the model is based on an experiment with turbulent air flow in a central channel of fuel bundle dummy consisting of a 37 heated rods. [2] The rod bundle consists of a triangular array at two different pitch-to-diameter ratios (P/D = 1.12 and P/D = 1.06) later referred to as PD112 and PD106. Measurements were performed with a hot-wire probe with x-wires and an additional temperature wire. Time-mean velocities, time-mean fluid temperatures, wall shear stresses and wall temperatures, turbulent quantities such as the turbulent kinetic energy, all Reynolds stresses and all turbulent heat fluxes were measured at the two different pitch-to-diameter ratios in a central channel of the bundle.

Figure 2-1. The cross-section area of the bundle of

heated rods used in the experiment.

2.1 Experimental apparatus

The experimental setup consists of a rod bundle of 37 parallel rods, each with a diameter of 140 mm, arranged in a triangular array in a hexagonal symmetric channel (see figure 1). The orientation of the channel is horizontal. The total length of the working section is L=11.5 m with an unheated entrance length of Liso = 4.60 m and a heated length of Lheat = 6.90 m. The

maximum pitch-to-diameter ratio of the rods is P/D 1.12. The rods are made of epoxy reinforced with fiberglass, sheathed with a 50 μm foil of monel metal, which serves as a Ohmic heating element. The foil is heated by direct current to temperatures in the range 60-110°C. Because the thickness of the metal foil has very accurate thickness (±2μm), the heat flux is uniform around the perimeter of the rods. The heat conduction towards the interior of the rods is very small due to the small thickness of the metal foil and the low conductivity of the rod material. Thus, the circumferential temperature variations due to different heat transfer will not be eliminated by conduction. The wall heat flux was determined from the measurements of the current and the voltage drop along the rods with an estimated error of ±1.5%.

The channel walls are made of aluminum covered with a thick insulation at the outside to minimize the heat losses. The whole bundle is made up of five sections, each 2.30 m long. The rod gap spacers were made of 4.0 mm thick and 15 mm wide (in axial direction) steel with rounded edges. Due to extremely small manufacturing tolerances, the deviations from the nominal bundle geometry are less than 0.2 mm, including the bending of the rods. At the entrance the fluid is air at atmospheric pressure and room temperature. A centrifugal blower drives air through the rod bundle. Before entering the working section, the air passes through a filter to remove particles greater than 1 μm and an entrance section of 5 m length with a honeycomb grid and a number of fine grid screens.

(11)

The measurements are performed at a position 20 mm upstream of the outlet. The time mean values of the axial velocity and the wall shear stresses are measured by Pitot and Preston tubes (OD: d = 0.6 mm) respectively; the mean temperatures are measured by sheathed thermocouples (OD: d = 0.25 mm). The wall temperatures are measured by an infrared pyrometer (Heimann KT 4), which has a target spot of 5 mm diameter at 1 m distance. Due to calibration difficulties, the uncertainty is ±1 K. The turbulent quantities are measured by hot wire anemometry using a three-wire probe. This probe consists of an x-wire probe with an additional cold wire perpendicular to the x-wire plane for simultaneous measurement of two components of instantaneous velocity and temperature. The x-wires have a length of 1.1 mm, a diameter of 2.5 μm and a spacing of 0.35 mm. The cold wire has a diameter of 1 μm, a length of 0.9 mm and is positioned 0.1 mm upstream of the x-wire prong tips. The measurement volume is approximately 1mm3

Subchannel P/D =1.12 P/D =1.06 Amount of measurement positions 147 119 Ub (m s-1) 20.57 20.63 Tb (C°) 39.7 47.0 TEin (C°) 12.3 5.8 Dh (mm) 53.6 33.5 Reb 64590 38754 qw (kW m-2) 1.39 0.98

Table 2.1-1. Experimental parameters that were used for the model calculations described in

this thesis [2]

2.2 Aspects of OpenFOAM

The OpenFOAM CFD Toolbox can simulate complex fluid flows involving chemical reactions, turbulence and heat transfer among other features.

The core technology of OpenFOAM is a flexible set of C++ modules. These are used to build a wealth of: solvers, to simulate specific problems in engineering mechanics: utilities, to perform pre- and postprocessing: libraries, to create toolboxes that are accessible to the solvers/utilities, such as libraries of physical models.

OpenFOAM is supplied with several pre-configured solvers, utilities and libraries. The program is open, not only in terms of source code, but also in its structure and hierarchical design, so that its solvers, utilities and libraries are fully extensible.

OpenFOAM uses finite volume numerics to solve systems of partial differential equations (PDE) ascribed on any 3D unstructured mesh of polyhedral cells. The fluid flow solvers are developed within a robust, implicit, pressure-velocity, iterative solution framework. [1]

A central theme of the OpenFOAM design is that the solver applications, written using the OpenFOAM classes, have a syntax that closely resembles the partial differential equations being solved. For example equation Eq. 2.2.1

p U U t +∇⋅ −∇⋅ ∇ =−∇ ∂ ∂ρ _φ _μ (2.2.1)

(12)

is represented by the code solve ( fvm::ddt(rho,U) + fvm::div(phi,U) - fvm::laplacian(mu,U) == - fvc::grad(p) ); [code 2.2-1] [3. s 68]

Each term in a PDE is represented individually in OpenFOAM code using the classes of static functions finiteVolumeMethod and finiteVolumeCalculus, referred to as fvm and fvc. The purpose of defining these functions within two classes rather than one is to distinguish:

- functions of fvm that calculate implicit derivatives of and return an fvMatrix<type>

- some functions of fvc that calculate explicit derivatives and other explicit calculations, returning a geometricField<type>

Geometry is defined in OpenFOAM in Cartesian coordinates. With such coordinates vectors are defined as (x1 x2 x3) and second order tensors such as Reynolds stress are defined as (R11

R12 R13 R21 R22 R23 R31 R32 R33). Dimensions are taken into account throughout calculations

and they are defined as [kg m s K A mol cd]. Kinematic viscosity for instance has the dimension m2s-1 and is defined in OpenFOAM as [0 2 -1 0 0 0 0] [3. s 24].

This report id divided into 3 parts. In chapter 3, the possibility to describe the velocity profile of an adiabatic flow will be investigated when exposed to the same conditions as in the experiment. In chapter 4 heat calculations are added to the cases treated here and in chapter 5 new geometries are introduced to the models obtained.

(13)

3 Adiabatic flows simulations

3.1 Adiabatic flows

A specific domain is said to be adiabatic if there exists no heat flux across the boundaries of this domain. In this context the flow is referred to be adiabatic because there is no heat flux between the rods and the channels.

The modelling of adiabatic flows mainly consists of finding a geometry that resembles the experimental setup and find equations used by OpenFOAM and validate them for these cases. In the following sections of heated flow and 2x2 and 3x3 geometry the same theory is used as in the previous sections and therefore only the added theory needed for that part is presented. 3.1.1 Junction of experiment and computer based simulation

The experiment is done with air flow in an attempt simulate incompressible water flow in a reactor. It is known from similar experiments that the compressibility of the fluid is negligible approximately below 0.2⋅mach[9] which is satisfied in this research. In the model the density of air is set to an average fixed value that best represents the full fluid body.

Geometry to be used in the computer model

Because the rod bundle consists of many rods arranged in a symmetrical geometry it is possible to construct the model based on a center rod with its adjacent spacing as channel for the air. The symmetry around the rod also allows the model to be made out of one twelfth of the whole rod perimeter. This is illustrated in Figure 3.1.1-1.

Figure 3.1.1-1. The geometry is extracted from a much greater perspective. In the left figure the center rod of

figure 1 is illustrated. The middle figure illustrates a magnification of one twelfth of the center rod, where b ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ _{⋅ D} D P 2 1

is the base, h

(

b⋅tan

( )

π 6

)

is the height, D is the diameter and ϕ_min =0 and ϕ_max =π₆. The picture on the right illustrates the air channel surrounding a portion of the heated rod with three symmetrical sides.

(14)

3.1.2 Governing equations

Continuum equation

This equation is developed by writing a mass balance over a volume element dxdydz, fixed in space, through which a fluid is flowing.

⎭ ⎬ ⎫ ⎩ ⎨ ⎧ − ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ = ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ out mass of rate in mass of rate mass of increase of rate

Figure 3.1.2-1 One arbitrary cell in the geometry

through which a mass flux runs through the sides normal to the x-direction.

) (

)

_z _dz

]

dydz t dxdydz ₊ + + + − + − − = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ρ _ρ _ρ _ρ _ρ _ρ _ρ (3.1.2-1) By dividing the entire equation by dxdydz and taking the limit as dx, dy and dz approache zero, and using the definitions of the partial derivatives one obtains

⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∂ ∂ = ∂ ∂ z y x U z U y U x t ρ ρ ρ ρ . (3.1.2-2) This is the equation of continuity, which describes the time rate of change of the fluid density at a fixed point in space. This equation can be written more concisely by using vector notation as follows

(

U

)

t ρ ρ ₌₋_∇_⋅ ∂ ∂ . (3.1.2-3)

(

ρUx

)

_x

The momentum equation

To obtain the equation of motion the momentum balance is written over the volume element ΔxΔyΔz as follows:

⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ + ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ − ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ = ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ fluid the on force external out momentum of rate in momentum of rate momentum of increase of rate

Figure 3.1.2-1. One arbitrary cell in the

dydz ₊ + + + − + − −φ φ φ φ φ φ (3.1.2-5)

(16)

for the net rate of addition of x-momentum across all three pairs of faces. Because the gravitational force is neglected in the model the external force, ρgxdxdydz, is not implemented

in the calculations. The sum of the x-component terms must be equated to the rate of increase of x-momentum within the volume element:

(

)

t U dxdydz x

∂

∂ ρ . When this is done, the x-component of the momentum balance is obtained. When this equation is divided by dxdydz and the limit is taken as dx, dy and dz approach zero, one obtains

⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∂ ∂ − = ∂ ∂ zx yx xx x z y x U tρ φ φ φ . (3.1.2-6) Similar equations are obtained for the y- and z-components.

By using vector-tensor notation, the three equations can be written as follows:

[

]

i i U tρ =−∇⋅φ ∂ ∂ (3.1.2-7) Where i = x,y,z

This results in the following formula when the ith component of Eq 3.1.2-7 is multiplied by the unit vector in the ith direction and the three components are added together vectorially:

[

φ

]

ρ =−∇⋅ ∂ ∂ U t (3.1.2-8) This is the differential statement of the law of conservation of momentum. The combined momentum flux tensor Φ is the sum of the convective momentum flux tensor ρvv and the molecular momentum flux tensor π, and that the latter can be written as the sum of p and τ, δ when no external forces are present. When this is inserted the following equation of motion is obtained

[

ρ

]

[

τ

]

ρ =−∇⋅ −∇ − ∇⋅ ∂ ∂ p U U U t (3.1.2-9) where the left-hand side of the equation is the rate of increase of momentum per unit volume. On the right-hand side the first term is the rate of momentum addition by convection per unit volume, the second and third terms are the rate of momentum addition by molecular transport per unit volume.

[4. § 3.2]

3.1.3 Theory of turbulence modelling

Estimating the turbulence

In order to get an idea of how turbulent the flow is with a Reynolds number of approximately 50000 a case is studied where a flow passes a cylinder at different critical values of Re number. At Re numbers in the order of 104 there is a disorderly fluctuating motion (turbulence) in the wake of the cylinder. At Re close to 106, turbulence appears upstream of the separation point, and the wake abruptly narrows down downstream.

(17)

The geometry that will represent the experiment does not carry the geometric symmetry of the cylinder mentioned above, but the latter it mentioned in order to obtain an understanding of the magnitude of turbulence outbreak.

[4. §3.7]

Strategy for solving turbulence problems

OpenFOAM is capable of calculating first order closure turbulence models, second order closure turbulence models and large Eddy simulations (LES).

One should be aware of the fact that there is no single turbulence model, which reliably can predict all kinds of turbulent flows. Each of the models works satisfactorily in the case of attached boundary layers, but may fail completely for separated flows. Thus, it is important to consider whether the model includes all the significant features of the flow being investigated or not. point that should be taken into consideration is the computational effort versus the accuracy required by the particular application. In some cases, a simple turbulence model can predict some global measures with the same accuracy as a more complex model.

Reynolds averaging

This methodology is based on the decomposition of the flow variables into a mean and a fluctuating part. The Navier-Stokes equation for momentum is then solved for the mean values (Eq. 3.1.3-4). For engineering applications, often only the mean flow is of interest. Considering incompressible flows, the velocity components and the pressure are substituted by: i i i U U U = + ′, p= p + p′ (3.1.3-1 and 3.1.3-2) For all averaging approaches, the average of the fluctuating part is zero; however,

0 ≠ ′ ′ j iU

U if both turbulent velocity components are correlated. [5. chapter 7.1.1]

Reynolds-Averaged Navier-Stokes Equations

If the terms in Navier-Stokes equations are averaged the following relations for the mass and momentum conservation are obtained:

0 = ∂ ∂ i i x U (3.1.3-3) ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ₋ ′ ′ ∂ ∂ + ∂ ∂ − = ∂ ∂ + ∂ ∂ j i ij j j j i j i U U x x p x U U t U ρ τ ρ ρ (3.1.3-4) These are known as the Reynolds-Averaged Navier-Stokes equations (RANS). This is the same as Navier-Stokes equations except for the term

(18)

(

i j i j

)

j i R ij =−ρ U′U′ =ρ UU − U U τ (3.1.3-5) which constitutes the so-called Reynolds-stress tensor. It represents the transfer of momentum due to turbulent fluctuations. The laminar viscous stresses are calculated using Reynolds-averaged velocity components, i.e.

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ = i j j i ij x U x U μ τ . (3.1.3-6) In the 3D case the Reynolds-stress tensor consists of nine components

( )

⎥⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ = ′ ′ 2 3 2 3 1 3 3 2 2 2 1 2 3 1 2 1 2 1 U U U U U U U U U U U U U U U U Ui j ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ . (3.1.3-7) However, because U ′ and i U ′ can be interchanged, the Reynolds-stress tensor contains only j

six independent components. The sum of the normal stresses divided by density defines the turbulent kinetic entry, i.e.

( )

[

]

. 2 1 2 1 2 3 2 2 2 1 U U U U U k = _i′ _i′ = ′ + ′ + ′ . (3.1.3-8) As is seen, the fundamental problem of turbulence modelling based on the Reynolds averaged Navier-Stokes equations is to find six additional relations in order to close the equations. [5, chapter 7.1.3]

Problem solving using the k-epsilon model

The first order closures are based mostly on the Eddy viscosity hypothesis of Boussinesq which postulates that the Reynolds stress tensor is proportional to the mean strain rate tensor (Eq. 3.1.3-9), but for certain applications also on non-linear Eddy-viscosity formulations.

( ) ( ) dy U d _x t t yx μ τ =− (3.1.3-9) k - ε is a two equation model and is based on solving the Navier-Stokes equations using the Boussinesq Eddy viscosity assumption (Eq. 31.3-9). Normally the Eddy viscosity, μt, is

computed from the two transported turbulence variables, k and ε. The Boussinesq assumption is both the strength and the weakness of two equation models. This assumption is a huge simplification that allows one to think of the effect of turbulence on the mean flow in the same way as molecular viscosity affects a laminar flow. The weakness of the Boussinesq assumption is that it is not in general valid. There is nothing that says that the Reynolds stress tensor must be proportional to the strain rate tensor. It is true in simple flows like straight boundary layers and wakes, but in complex flows, like flows with strong curvature, or strongly accelerated or decelerated flows the Boussinesq assumption is simply not valid. This

(19)

gives two equation models inherent problems to predict strongly rotating flows and other flows where curvature effects are significant.

[6, two equation turbulence models]

Problem solving using the Reynolds stress model (RSM)

The Reynolds stress model is a higher level, second order closure turbulence model. In RSM, the Eddy viscosity approach has been discarded and the Reynolds stresses are directly computed. The exact Reynolds stress transport equation accounts for the directional effects of the Reynolds stress fields.

The Reynolds stress model involves calculation of the individual Reynolds stresses, u ′i′uj ,

using differential transport equations. The individual Reynolds stresses are then used to obtain closure of the Reynolds-averaged momentum equation.

[6, Reynolds stress model]

Implementation of momentum equation in OpenFOAM

In OpenFOAM the momentum equation is stated as

( )

UU +∇⋅R =−∇p ⋅

∇ .

(3.1.3-10) This equation is divided by ρ. In the right-hand side of the equation p is the kinematic pressure, i.e., pressure divided by density. The time dependent term on the left-hand side of the equation has been removed because of steady state. The viscous stress, R, in the second term is defined as υeff∇Uwhere υeff is the kinematic viscosity, calculated from selected

transport and turbulence models. [7. s 54]

A solution to these equations is made in OpenFOAM by a sequence of calculations where initially an approximation of the velocity field is obtained by solving the momentum equation. The pressure gradient term is calculated using the pressure distribution from the previous iteration or an initial guess. The pressure equation is then formulated and solved in order to obtain the new pressure distribution.

[8]

3.2 Method

3.2.1 Choice of application

OpenFOAM has a number of pre-defined applications that fit certain physical cases. There are application solvers that solve flows with characteristics like:

• compressible flows, • incompressible flows, • heat Transfer,

• multiphase.

Each of these categories is a parent-directory for a number of solvers which are distinguished by characteristics such as time variability, turbulence and transport method. In this work it is necessary to choose a solver that can handle incompressible flow because the model is supposed to illustrate water flow in a reactor core. From the experimental data it is given that the flow has a Reynolds number in the scale of 104 which means that the flow is fully

(20)

turbulent. The application should also be chosen to be steady state because there are no changes in time while the calculation proceeds. The one solver that satisfies the requirements is simpleFoam which also does not need to be altered in order to calculate adiabatic flows with characteristics described above.

3.2.2 Creating a case

Within the tutorial there are many examples that are ready to run without alterations. The easiest way to create a new case is to copy an existing case, which can then be altered.. In this work two cases were created from the pitzDaily pre-defined case, PD112 and PD106 for the cases with pitch to diameter ratios of 1.12 and 1.06 respectively.

Choice of kinematic viscosity

The turbulence of the flow needs to be defined given the Reynolds number in the experiment report [2]. The kinematic viscosity should be entered into the file named transportProperties under the folder named constant. In Table 3.2.2-6 it is shown how the kinematic viscosity is calculated together with its value for the two P/D ratios.

subchannel formula P/D = 1.12 P/D = 1.06 ν (m2_s-1₎ b h b D U Re ⋅ 1.707e-5 1.7833e-5

Table 3.2.2-1. Calculation of the kinematic viscosity for the two cases. Entering the geometry and mesh to the case

illustrated in fig. 3.2.2-1 for the geometries of PD112 and PD106.

Figure 3.2.2-1. The 3D geometry of the PD112 and PD106

cases with vertex numbering. The x-, y- and z- directions are used throughout the handling of these cases.

The geometry of the case is entered in the constant/polyMesh/blockMeshDict file. The description of blockMeshDict file is listed in table3.2.2-1.

(21)

entry Description

convertToMeters all coordinates are multiplied with this given value

vertices These are the vertexes that are the corners of the mesh and are numerically ordered with the first one having index zero. The vertices are given in Cartesian coordinates as (x1,x2,x3)

blocks Each block definition is a compound entry consisting of a list of vertex labels, a vector giving the number of cells that shall be created in the mesh in each direction, the type and list of cell expansion ratio in each direction

edges Each edge joining 2 vertex points is assumed to be straight by default. However any edgy may be specified to be curved by entries in a list named edges. For this case arcs are required on one of the patches which need to be specified with start end and a midpoint.

patches The patches define the boundary of the given geometry with a list of vertexes as input. It is also here the boundaries are given appropriate names for recognition.

mergePatchPairs This function is only used in cases where the geometry consists of more than one block which it does not in these cases being handled.

Table 3.2.2-3. Specification of the contents of blockMeshDict. See appendix for specific entries for the cases PD112

and PD106.

When setting the blocks entry it is important to place the gridline that is the closest to the wall at an appropriate distance. The reason for this is to place the boundary layer within the first layer of cells so that the second layer consists only of flow with turbulent characteristics. In order to calculate this distance the normalized distance,y+ , is used which is defined by the following formula: υ τ u Y y+ = ⋅ . (3.2.2-1)

Where Y is the actual distance from the wall and uτ is the friction velocity defined by

ρ τ_w

.

τw is the wall shear stress which is defined by

2 Re 0791 . 0 2 25 . 0 u w ⋅ ⋅ = ρ

τ . Since the actual distance

is required the formula is rearranged in the following way:

ρ τ υ w y Y = ⋅ + (3.2.2-2) The calculation of the required Y value for the different turbulent models and cases is attached in the appendix.

Setting the grid for k – ε and RMST turbulence models

It is a rule of thumb to have a y+ value between 20 and 50 [9] in order to satisfy the placement of boundary layer for k – ε and RMST turbulence models. Using parameters given in table 2.1-1 it is possible to calculate at what distance to put the first gridline which is shown in table 3.2.2-4.

subchannel PD112 PD106

Y (m) 3.33e-5 3.26e-5

(22)

Since it is required that the boundary layer is in the middle of the first cell layer it is necessary to determine the amount of cell layers, cl, based on the thickness of the thinnest part of each PD case. The amount of cell layers is calculated as shown in table 3.2.2-5 usingCL=sT/

(

2⋅Y

)

.

Subchannel PD112 PD106

Length of shortest side, sT, 8.4 4.2

sT/Y 12.60 6.45

Chosen amount of cell layers, cl, (8 12 120) (8 6 80)

Table 3.2.2-5 Setting of the block entry for k - ε and RSTM turbulence models

The cl value is chosen to be the lower integer value do to caution not to set the boundary layer outside the first cell layer.

Setting the grid for low-Re turbulence models

A rule of thumb for setting the first gridline for a low-Re turbulence model is to choose a value y+ around 1. Here it has been decided to use a y+ value of 0.5 in order to get as high resolution as possible. It has been chosen to use grading for these models to avoid building a mesh with unnecessary many cells. The grading, G, is defined as the length of the last cell width,t2, divided by the first cell width,t1. The number of cell layers, cl, has been chosen as

the same number as the grading for easier calculation of the required Y value. The side length of the geometry is calculated by sT i G t2

n i ⋅ ⋅ =

∑

. But since ₂ 1 t₁ G t = ⋅ it is possible to formulate the side length as =

∑

n ⋅

i

t i

sT ₁. Rearranging this formula gives the desired width of the first cell as shown in equation 3.2.2-1.

∑

= _cl i i sT t₁ (3.2.2-3) As before sT is chosen to be the shortest side length for corresponding PD ratio.

Subchannel PD112 PD106

y+ _{0.5 0.5}

Y (mm) from eq. 3.2.2-2 0.00833 0.00814 Chosen G and cl value 30 30 gridline position 0.018065 0.018065 Chosen amount of cell layers (8 30 120) with grading

30 in y-direction

(8 30 80) with grading 30 in y-direction

Table 3.2.2-6 Setting of the block entry for low-Re turbulence models

This placement of the first gridline is at a greater distance than the calculated position of the boundary layer which is twice the y value. The boundary layer is therefore within the first cell layer for both PD ratios.

(23)

Subchannel PD112 PD106 Mesh for k – ε

models with wall functions and RSTM model with wall functin

Mesh for low-Re turbulence models

3.2.2-7 Mesh used for different turbulence models.

After the blockMeshDict file has been edited and saved the mesh is generated from this dictionary file by entering the following command in the terminal window:

blockMesh <root> <case>

For example, if the geometry of PD112 is to be read and the prompt is in the same directory as the root for the case the correct command is:

blockMesh . PD112 Viewing the mesh

The principle behind blockMesh is to decompose the domain geometry into a set of one or more 3D, hexahedral blocks. After the geometry is read it is understood by OpenFOAM from the point, faces and boundary files which has been created by this action. OpenFOAM is supplied with a post-processing utility named paraFoam that uses ParaView, an open source visualization application. ParaFoam is a script that launches paraView using the reader module supplied with OpenFOAM. It is executed like any of the OpenFOAM utilities with the root directory path and the case directory name as arguments The graphical viewer of geometry and results is named paraFOAM. After the mesh has been read by typing the blockMesh command above it is possible to look at it by typing:

(24)

Define the physical attributes of the boundary

When the blockMeshDict file has been read, a file named boundary can be viewed in the same directory as the blockMeshDict file. In this file the six patch names hole, sym1, sym2 sym3, inlet and outlet, which are chosen in blockMeshDict, are listed. It is necessary to define the type and physicalType of these patches. The type of patch is described purely in terms of geometry or a data ‘communication link’. PhysicalType is a type describing the boundary conditions in the physical world on one or more fields. These settings of the boundary file are illustrated in table 3.2.2-3

boundary type physicalType description

hole wall wallFunction For cases which require wall turbulence modelling, a wall must be specified with a wall patch type, so that the distance from the wall of the cell centers next to the wall are stored

sym1 symmetryPlane symmetryPlane For a symmetry plane sym2 symmetryPlane symmetryPlane -||-

sym3 symmetryPlane symmetryPlane -||-

inlet patch inlet The basic type for a patch condition that contains no geometric or topological information about the mesh e.g an inlet or outlet

outlet patch pressureOutlet -||-

Table 3.2.2-8 Description of how to set the six boundaries in these cases. Setting of the initial and boundary conditions

The initial and boundary conditions are set in the t=0 dictionary of the newly created case.

The parameters that need to be set are turbulent dissipation rate (epsilon), turbulent kinetic energy (k), pressure (p), Reynolds stress (R) andvelocity (U). Some parameters that should be set to zero are instead set to a small value in order to avoid singularities in the calculations. The parameters that shall be set to a certain value are shown in table 3.2.2-9.

parameter source value

dependent boundary P/D=112 Non RSM/RSM P/D=106 Non RSM/RSM epsilon 43 42 1 RSM non h 2 / 3 D 0.3⋅ k [9]

(

)

4 4 4 3 4 4 4 2 1 RSM 2 3 33 22 11 3 . 0 Dh R R R ⋅ + + internalField/ _inlet 48.39/136.87 78.11/219.00 k ₀_.₀₀₂ 2 b U ⋅ [9] internalField/ inlet 0.846 0.851 p [9] internalField/ outlet 0 0 R k R R R 3 2 33 22 11 = = = [9] internalField/ inlet (0.564 0.001 0.001 0.001 0.564 0.001 0.001 0.001 0.564) (0.571 0.001 0.001 0.001 0.571 0.001 0.001 0.001 0.571) U [2] internalField/ inlet (0 0 20.57) (0 0 20.63)

(25)

The setting of the initial conditions and boundaries is a mixture of using functions and numerical values. The total settings of these parameters for the different boundaries are shown in table 3.2.2-10

parameter Hole Sym1 Sym2 Sym3 Inlet outlet

epsilon zeroGradient symmetryPlane symmetryPlane symmetryPlane Table 3.2.2-4 zeroGradient k zeroGradient symmetryPlane symmetryPlane symmetryPlane Table 3.2.2-4 zeroGradient p zeroGradient symmetryPlane symmetryPlane symmetryPlane zeroGradient Table 3.2.2-4 R zeroGradient symmetryPlane symmetryPlane symmetryPlane Table 3.2.2-4 zeroGradient U Table 3.2.2-4 symmetryPlane symmetryPlane symmetryPlane Table 3.2.2-4 zeroGradient

Table 3.2.2-10 Description of the initial and boundary values and functions set to the property parameters in these

cases.

Choice of turbulent model

OpenFOAM has twelve different turbulence models to choose from in order to find a flow with parameters that resembles the experimental flow. The turbulence models are all located in the file named turbulenceProperties which lies under the constant dictionary. There are eight different, first-order turbulence models that are based on k-epsilon theory, four of which are low Re. There are two second-order turbulence models that require Reynolds stress as input.

There are also two other turbulence models that are based on other inputs than k, epsilon and the Reynolds stress tensor such as turbulent dynamic viscosity, η.

To choose one of these models type the name of the model after turbulenceModel and cheek that turbulence is marked on

Setting the numerical schemes

The fvSchemes file in the system directory sets the numerical schemes for terms involving functions such as gradients, divergences and Laplacians in equations that appear in applications being run. In this case all the calculations are done with Gaussian finite volume integration. Gaussian integration is the most common and is based on summing values on cell faces which must be interpolated from cell to cell.

Choosing solvers

The equation solvers, tolerances and algorithms are controlled from the fvSolution dictionary in the system directory. There are three types of entries in this dictionary. The first entry is for choosing which solvers that should be used for the different variables in the calculation. The second entry decides whether the calculation needs to be transient or steady state. The third entry lets the user decide the relaxation factor of the different variables. [3. s 113]

The solver entry specifies each linear-solver that is used for each discretisation for the used equations; it is emphasized that the term linear solver refers to the method of number-crunching to solve the set of linear equations. There are totally six different solvers to choose from. [3. s112]

In this case all solvers are chosen to be of the sort Incomplete-Cholesky conjugate gradient (ICCG) except for the pressure solver, which is chosen to be of algebraic multi-grid (AMG). The only alteration from the default settings is to choose AMG as solver for the pressure variable instead of ICCG. The AMG solver works by making a quick solution on a smaller amount of cells and then using this solution as an initial guess for a calculation on a finer mesh. Therefore AMG has an extra argument that needs to be set which is an approximate mesh size at the most coarse level in terms of the number of cells <nCells>, following the tolerance parameters. [3. s113]

These sparse matrix solvers are iterative, i.e., they are based on reducing the equation residual over a succession of solutions. The residual is ostensibly a measure of the error in the solution

(26)

so that the smaller it is, the more accurate the solution. More precisely, the residual is evaluated by substituting the current solution into the equation and taking the magnitude of the difference between the left and right hand sides; it is also normalized in order to make it independent of the scale of problem being analyzed. Before solving an equation for a particular field, the initial residual is evaluated based on the current values of the field. After each solver iteration, the residual is re-evaluated. The solver stops if either of the following conditions is reached:

• the residual falls below the solver tolerance, <tol>

• the ratio of current to initial residuals falls below the solver relative tolerance, <relTol>.

The solver tolerance represents the level at which the residual is small enough so that the solution can be deemed sufficiently accurate. The solvers relative tolerance limits the relative improvement from initial to final solution. [3. s 113]

The relative tolerance has been set to zero in order to totally control the solver with the requirement of tolerance.

Choosing time dependence

The second entry in fvSolution allows the user to choose whether the case shall be handled as transient or as steady state. In this work only steady-state flow is considered. Therefore the word SIMPLE (semi-implicit method for pressure-linked equations) should be typed as a name of the next function in fvSolution. This algorithm is an iterative procedure for solving equations for velocity and pressure. An additional correction to account for mesh non-orthogonality is available in SIMPE in the standard OpenFOAM solver applications. A mesh is orthogonal if, for each face within it, the face normal is parallel to the vector between the centers of the cells that the face connects. The number of non-orthogonal correctors is specified by the nNonOrthogonalCorrectors keyword. The number of non-orthogonal correctors should correspond to the mesh for the case being solved, i.e., 0 for an orthogonal mesh and increasing with the degree of orthogonality up to, say, 20 for the most non-orthogonal meshes.

Setting of controlDict

Input data relating to the control of time and reading and writing of the solution data are read from the controlDict dictionary. First the start/stop times and the time step are set for the case to run. It is desired to start running at time t = 0, which means that OpenFOAM needs to read field data from the directory named 0 that is already set. Therefore the startFrom keyword is set to startTime and the startTime keyword is set to 0. For the end time it is desired to reach the steady state solution where the flow characteristics is not changed with further running. The time step, deltaT, and end time stopTime are only set on account of how many iterations are needed for a satisfying solution. In this case deltaT is set to 0.1 and stopTime is set to 25. The writeControl is set to timeStep and writeInterval is set to 50 in order to get 5 follow-up solutions during the calculation.

3.2.3 Running the case

The case is run by typing the following line in the terminal window in the same dictionary as the case.

<application> . <case>

When a solver is executed, it reports the status of equation solution, if the level debug switch is set to 1 or 2 (default) in DebugSwitches in the $HOME/.OpenFOAM1.3/controlDict file.

(27)

An example from the beginning of a case calculation of the simpleFoam. PD112 is shown in code 3.2.3-1:

Time = 10.4

BICCG: Solving for Ux, Initial residual = 0.0104723, Final residual = 0.000440457, No Iterations 1 BICCG: Solving for Uy, Initial residual = 0.00402328, Final residual = 0.000144148, No Iterations 1 BICCG: Solving for Uz, Initial residual = 0.0087, Final residual = 0.000370788, No Iterations 1 AMG: Solving for p, Initial residual = 0.0146053, Final residual = 9.88916e-07, No Iterations 221 time step continuity errors : sum local = 1.62748e-05, global = 5.99726e-07, cumulative = 0.000181624 BICCG: Solving for epsilon, Initial residual = 5.22562e-06, Final residual = 5.22562e-06, No Iterations 0 BICCG: Solving for Rxx, Initial residual = 3.80339e-09, Final residual = 3.80339e-09, No Iterations 0 BICCG: Solving for Rxy, Initial residual = 6.19749e-06, Final residual = 6.19749e-06, No Iterations 0 BICCG: Solving for Rxz, Initial residual = 3.33768e-06, Final residual = 3.33768e-06, No Iterations 0 BICCG: Solving for Ryx, Initial residual = 6.19749e-06, Final residual = 6.19749e-06, No Iterations 0 BICCG: Solving for Ryy, Initial residual = 4.53427e-09, Final residual = 4.53427e-09, No Iterations 0 BICCG: Solving for Ryz, Initial residual = 6.51193e-06, Final residual = 6.51193e-06, No Iterations 0 BICCG: Solving for Rzx, Initial residual = 3.33768e-06, Final residual = 3.33768e-06, No Iterations 0 BICCG: Solving for Rzy, Initial residual = 6.51193e-06, Final residual = 6.51193e-06, No Iterations 0 BICCG: Solving for Rzz, Initial residual = 5.91561e-09, Final residual = 5.91561e-09, No Iterations 0 ExecutionTime = 248.9 s ClockTime = 270 s

[code 3.2.3-1] It can be seen that, for each equation that is solved, a report line is written with the solver name, the variable that is solved, its initial and final residuals and number of iterations that are needed to satisfy the tolerance conditions.

3.2.4 Display a cross section view of the results

In order to view the calculated results paraFOAM is launched as described in 3.2.6 subtopic view the mesh. In this case it is of interest to observe the results that have been saved in the latest time dictionary. The choice of time is done in the panel named Parameters together with region, vol Field and point Field. To view the results at a certain time the last time alternative box should be checked followed by pressing the accept button. Now the 3D geometry is displayed. The experimental data, which is used to validate the model, is only sampled at the outlet end of the channel. Therefore only this set of data is visualizedIn paraFoam the cut function is used to get a view of the cell values at a cross section near the end of the channel only. The user has to define at what z value of the fluid bodythe cut should be made. The z-value is chosen as the maximum length of the geometry, i.e., 11.48.

3.2.5 Make samplings of theoretical data

At this point theoretical data for all mentioned parameters are gathered for the whole geometry. It is suitable to compare the experimental data with theoretical data along the same sides. These sides are chosen to be the narrowest part of the channel, i.e., ϕ =0 and the thickest part of the channel, i.e., ϕ =30, which will be called sides 1 and 2 in the remainder of the report.

In order to get hold of the data at these sides it is necessary to use the sample function in OpenFOAM. This is done by entering the file located in OpenFOAM-1.3/applications/utilities/postProcessing/miscellaneous/sample which is named sampleDict. This file is to be copied and relocated in the system directory of the case of interest. The sampleDict file needs to be edited for purposes that fit these cases and geometries. The writeFormat is changed to raw 2 for data access for the plotting program gnuplot. The next set

(28)

for side 2. nPoints stands for the number of sampling points which have been arbitrary chosen for each case in order to get good looking plots. In some cases the number of sampling points has been altered to avoid gaps from omitted data. The last entry in sampleDict is named fields. In this entry it is possible to alter the sampling command to only sample the z-component, Uz, of the

magnitude of the fluid velocity, U . This is done by entering U(2) within the fields brackets. After the sampleDict file is altered the sampling is initiated by the command:

sample <root> <case>

This creates an extra dictionary “samples” in the case. This dictionary holds a time dictionary for every write interval. The sampling data from the latest time dictionary needs to be copied and pasted into a folder for results concerning the case. The copied file, with sampling data, consists of the measured velocity components as function of the distance from the wall.

3.2.6 Extraction of experimental data

In this section the experimental data from Ref. 2 will be studied and compared with the data from the model made in OpenFOAM by using an already existing solver named simpleFoam.

Sections of figure 3.2.6-1 has been magnified and clarified with extra help lines in order to make the iso-lines easier to read and extract numeric values from. The four sides, 1 and 2, for both PD ratios are shown together with tables of extracted data consisting of wall distance and velocity values. The extractions of the experimental values are attached in the appendix.

Figure 3.2.6-1 The distribution of time mean

velocity in a central channel with P/D =1.12 (Ub = 20.57 m s-1) and

P/D = 1.06 (Ub = 20.63 m s-1).

The numbers 1 and 2 represented in the figure shows two sides of the geometry were 2D graphs shall compare experimental data with data obtained from the calculated model.

1 1

2 2

(29)

3.2.7 Plot experimental data with theoretical data with gnuplot

In order to plot the experimental data next to the theoretical data from the OpenFAOM model the experimental data needs to be modified so that it is represented in the same manner as the theoretical data. This is done by creating a file in the same directory as the theoretical data with the wall distance in a first column, the velocity in the second column and the desired error bar width in a third column.

The plotting program gnuplot is invoked by the command

gnuplot –background white.

An example of how to make a plot similar to the ones presented in this report is shown in code 3.2.7-1. These command lines need only to be pasted after the prompt to execute the plotting command. The plot is made with an auto scale function so that no ranges need to be defined. The first lines, separated by comma signs, define which files that are requested to be plotted. Here “uz” stands for z-component of velocity and e for experimental data. “w” is a command which tells gnuplot that an extra function is requested for the plotting of the current file. In these cases it is requested to represent the data with figured points with intersected lines, which is done with the command linespoints. In order to add error bars to the experimental data it is necessary to enter the file twice, but using title on only one of them in order to get it mentioned once in the key for the figure. The set functions are very straight forward to use for altering the design of the plotting window.

The three last lines are for saving the plot. First the terminal needs to be set to a format that is compatible with the picture format. Next the plot is saved. The last command line replot is necessary to get the plot written.

plot 'uze.xy' title 'experimental data' w lines,'uze.xy' notitle w errorbars,'uznonLinear.xy'title 'T.M uznonLinear' with lines, 'uzLienCubic.xy'title 'T.M uzLienCubic' with

lines,'uzLRR.xy'title 'T.M uzLRR' with

lines,'uzRNGkepsilon.xy'title 'uzRNGkepsilon' with lines,'uzkepsilon.xy'title 'uzkepsilon' with lines set title "Velocity profiles of P/D=1.06, side 1" set ylabel "axial velocity component of flow uz (m/s)" set xlabel "distance from heated rod, y (m)"

set key top left set grid

replot

set terminal jpeg

set out "uzPD106s1.jpeg" replot

[code 3.2.7-1] In order to plot the other graphs it is convenient to give them the same file names but

separating them by appropriate dictionaries names. Then it is only necessary to change the title names when selecting plots of the different cases and sides.

(30)

3.3 Results of adiabatic flows

In this section the velocity profiles for the turbulence models (marked T.M in the plots) that have been able to generate results that resemble the data from the experiment will be presented. The turbulence models that are not represented in this section have diverged or required input data than the data that was not available. The results are presented in plots where the z-velocity component of the flow has been plotted against the wall distance. The experimental data is represented together with an error bar of 1 m/s in both positive and negative direction in order to better judge the differences between data sets. The lines are linearly interpolated between sampling points. Larger gaps are displayed where sampling data is missing.

Results from k-epsilon and RSTM models with wall functions

Figure 3.3-1 Axial velocity distribution of PD112 side 1

(31)

(32)

Results from Low Re models

Figure 3.3-5 Axial velocity distribution of PD112 side 1 Figure 3.3-6 Axial velocity distribution of PD112 side 2

Figure 3.3-7 Axial velocity distribution of PD106 side 1 Figure 3.3-8 Axial velocity distribution of PD106 side

2

In order to get an idea of how the z-component of the velocity profile is distributed between the sides 1 and 2 a cross-section graph of the outlet parameter has been presented in the appendix (cross-section plots). These graphs have not been compared to any experimental data but they are of importance in the respect that they illustrate that the models seem to work well.

(33)

3.3 Discussion

It is difficult to say which of the k - ε turbulence models that in general are the closest to the experimental data. The plots revile however that the three turbulence models kEpsilon, RNGkEpsilon and LRR are close to each other. The LRR model is however the most correct one among these three. The NonlinearKEShih and the LienCubicKE models are both often too deviated and less correlated to the experimental velocity values then the k - ε models. The plots of the two low-Re models clearly show that LienLeschzinerLowRe has the best results of all the models. The plot of PD106 side 1 is the plot that shows most deviating results. This is probably because the difficulty of calculating turbulence in the near wall region as mentioned in the theory of this section. The LRR and LienLeschzinerLowRe models are within the vicinity of the experimental data error bars in all the graphs except for the near wall region of figure PD106 side 2 where they are deviated but they approach the experimental data as the distance from the wall increases. The RNGkepsilon turbulence model shows reasonable good results except for side 1 of both PD112 and PD106 where it is more than twice the length of the error bar away from the experimental data.