Thermal Mixing CHT Simulations with OpenFOAM:URANS and LES

(1)

Thermal Mixing CHT Simulations with OpenFOAM:

URANS and LES

IGNACIO GALLEGO-MARCOS

Master’s Thesis at Nuclear Reactor Technology Supervisor: Prof. Henryk Anglart

TRITA-FYS 2013:41 ISSN 0280-316X ISRN KTH/FYS/–13:41–SE

(2)

(3)

iii

Abstract

In the Unit 3 of Oskarshamn Nuclear Power Plant, 2008, a control rod stem was detected to be broken. A study made by Forsmarks Kraftgrupp AB (Tinoco and Lindqvist (2009)) revealed that thermal fatigue was the origin of the cracks.

The objective of this project is to find the low frequency temperature oscillations that caused the control rod stem to break and to provide data for the experimental setup used in the THEMFE project at KTH university, Stock- holm.

The thermal mixing was simulated with OpenFOAM using the chtMultiRe- gionFoam solver. The turbulence model, mesh, numerical schemes, boundary and initial conditions were carefully chosen and adapted to our case.

The results from the URANS and LES allowed us to give a rough estimation of the length of the mixing zone of about 6 cm. Unfortunately, the thermal fatigue could not be properly studied due to the low amount of simulated seconds. Further conclusions will be made after a long term simulation.

(4)

(5)

Acknowledgements

I would like to thank all of the people who has helped me with this master thesis that Professor Henryk Anglart offered me. Starting with Roman Thiele (PhD student, KTH Stockholm) without his advices and help we would not have gone too far, Nicolas Forsberg (Forsmarks Kraftgrupp AB ) who was always there to reply our messages and David Haces (UPM Madrid), my friend and colleague at the university in Spain, it was a pleasure to work together with you.

Apart from the academical, I thank my friends for being such nice people. Toshi and Per from my corridor in Lappis, I could not have had better company than yours, Nicolas, Joshua and Riccardo from the beginning of the times in Stockholm (hope that you will come back! Otherwise, see you in our next trip) and the Tuesday dinner team with Yuki, Xueqi, Mai, Nick, Pierre, Takashi (with his Riga wine), Naoko and Yuuki, always at 9, with the best food. I will miss you all when you are gone after this summer. Thanks also to David and Claudio for our Spanish lunch time and for providing biscuits from the safety department.

Last, but not least, big thanks to my parents, brothers and grandmother, whom without I would not have been able to accomplish anything. I really enjoyed our talks via skype and your visits to Stockholm.

Special mention to my mum for going surfing with me for the first time when I was little =).

v

(6)

(7)

List of Figures

2.1 Normalized Kolmogorov energy spectrum . . . 13

4.1 Mesh of the simplified 2D model . . . 21

4.2 Low Reynolds LS k-Epsilon turbulence model validation. Left: Temper- ature profile. Right: Velocity profile . . . 22

4.3 Low Reynolds k-Omega SST turbulence model validation. Left: Tem- perature profile, Right: Velocity profile . . . 23

5.1 3D view of the physical model . . . 26

5.2 Hot inlet detail. With solid and fluid regions included . . . 29

5.3 Required cell size for the LES simulation . . . 30

5.4 Cell size used in the LES simulation . . . 31

5.5 Cell size at the hot inlet for the LES mesh . . . 32

5.6 Tetrahedral cells arrangement. Left: Tetra cells following the stream lines. Right: Inner cyl (yellow) - water (blue) inter phase . . . 33

5.7 Cell size used in the VLES simulation . . . 34

5.8 Cell size at the hot inlet for the VLES mesh. . . 34

5.9 Velocity profiles. Left: Hot inlet. Right: Cold inlet . . . 38

5.10 Turbulent properties profiles for the hot inlet: Left: k . Right: ω . . . . 38

5.11 Water propertites cp, k and µ fitted with polynomials . . . . 41

6.1 Probes location in the URANS and LES cases . . . 45

6.2 Temperature distribution from the RANS simulation . . . 46

6.3 y⁺ values for the inner cylinder . . . 47

6.4 Mean temperatures in the water and inner cylinder (URANS) . . . 48

6.5 Standard deviation in the water and inner cylinder (URANS) . . . 48

6.6 Temperautre oscillations at diferent positions of the pipe axis (URANS) 49 6.7 Mean temperatures in the water, 1mm away from the wall . . . 50

6.8 Mean temperatures in the water and inner cylinder (LES) . . . 51

6.9 Standard deviation in the water and inner cylinder (LES) . . . 51

6.10 Temperautre oscillations at diferent positions of the pipe axis (LES) . . 52

6.11 Vicinity of the hot inlet (y = 0.75m, t = 0.7s). Left: Velocity in the z direction. Right: Temperautre field . . . 53

xi

(12)

6.12 Downstream the hot inlet (y = 0.65m, t = 0.7s): Velocity in the z

direction. . . 53

6.13 Slices of the pipe at y = 0.65m, 4t = 0.4s . . . . 54

A.1 Orthogonality, 2D graphical descption . . . 63

A.2 Skewness defined by OpenFOAM, 2D graphical description. . . 64

B.1 Problem with O-Grid. Left: O-Grid mesh. Right: Section of the velocity profile: . . . 65

B.2 Solution to the O-Grid. Left: After Linking edges. Middle: After a second O-Grid. Right: Turbulent kinetic energy profile. . . 66

B.3 T-Junction mesh . . . 66

B.4 Bad behaviour of splitMeshRegions. Left blue: Water, left yellow: In- nercyl. Right: Explanation . . . 67

(13)

List of Tables

2.1 Closure coefficients for the Launder and Sharma k − model . . . . 7

2.2 Closure coefficients for the k − ω model . . . . 8

2.3 Constants used in Kolmogorov -5/3 Law . . . 13

4.1 Geometrical parameters of the 2D model . . . 22

4.2 Mean errors form the validation of LS kEpsilon and kOmegaSST Low Reynolds models . . . 23

5.1 Geometrical parameters of the physical model . . . 26

5.2 Geometrical parameters of the inlets and outlet . . . 26

5.3 Cell size needed for a y⁺= 1 at different locations . . . 28

5.4 Cell size needed for a y⁺= 30 at different locations . . . 28

5.5 URANS mesh properties . . . 29

5.6 LES mesh properties for the WATER region . . . 32

5.7 VLES mesh properties for the WATER region . . . 35

5.8 Discretization schemes for the URANS and LES simulations . . . 36

5.9 Solver schemes for the URANS simulation . . . 37

5.10 u, p, and T patch definition for the water region . . . 37

5.11 Wall boundary conditions for the turbulent properties . . . 39

5.12 Inlet and outlet boundary conditions for the turbulent properties . . . . 40

5.13 RANS residuals after 13000 iterations. . . 40

5.14 URANS run control parameters . . . 40

5.15 SS 316 thermo physical properties . . . 42

5.16 LES boundary conditions for µ^SGS_t . . . 42

5.17 LES and VLES run control parameters . . . 43

xiii

(14)

(15)

Nomenclature

Dimensionless parameters Co Courant number Pr Prandtl number Re Reynolds number

ReL Reynolds number associated with large eddies u⁺ Velocity at the near wall region

y⁺ Wall distance

Greek letters

αt Thermal diffusivity

β Thermal expansion coefficient 4 Cell size

Dissipation rate εⁿ Error after n iterations η Kolmogorov length scale

κ Von Karman constant; wave number κ^∗ Wave number for the 80-20% division λ Thermal conductivity

µ Dynamic viscosity ν Kinematic viscosity νt Turbulent eddy viscosity νef f ν + νt

ξ Space coordinates

π Pi number

ρ Density

ρ0 Reference density

σT Standard deviation of the temperautre τ^R Reynolds-Stress tensor

τ^SGS Sub-grid scale stress tensor τ_ω Wall shear stress

xv

(16)

ω Specific dissipation rate; weighting factor

Roman letters

A Matrix

B Matrix

cp Heat capacity

Cf Skin friction coefficient

E Energy

fL Production function fη Disspation function

g Gravity

G Filter function h Enthalpy; height i Complex number

k Specific turbulent kinetic energy l Eddy length scale

l₀ Maximum eddy length scale

l^∗ Eddy length scale for the 20-80% division

M Matrix

n Iteration number

N Matrix

prgh p − rgh

P Preconditioner matrix p⁰ Fluctuating pressure q⁰⁰ Heat flux

rⁿ Residual after n iterations

t Time

T Temperature

T Averaged temperature T0 Reference temperautre u Velocity

U Mean velocity u⁰ Fluctuating velocity u Mean velocity uτ Friction velocity U^∗ Intermediate velocity

Abbreviations

CHT Conjugate Heat Transfer DIC Diagonal Incomplete Cholesky DILU Diagonal Incomplete LU

FOAM Field Operation and Manipulation

(17)

List of Tables xvii

GAMG Geometric Algebraic Multi Grid KTH Kungliga Tekniska Hogskolan LES Large Eddy Simulation NWM Near Wall Modelling

PBiCG Preconditioned Bi Conjugate Gradient PCG Preconditioned Conjugate Gradient PISO Pressure Implicit Split Operator RANS Reynolds Averaged Navier Stokes SIMPLE Semi Implicit Pressure Linked Equations THEMFE Thermal Mixing and Fatigue Experiment URANS Unsteady Reynolds Averaged Navier Stokes VLES Very Large Eddy Simulation

(18)

(19)

Chapter 1

Introduction

BWR reactors introduce the control rods from the bottom of the reactor vessel.

These control rods are fixed to a solid cylinder called stem that allows the up and down movement.

In the Unit 3 of Oskarshamn Nuclear Power Plant, 2008, a control rod stem was detected to be broken. A further examination in the rest of the units revealed multiple cracks at the upper part of the stem of many control rods and also a broken stem of a control rod at Forsmark 3. Thermal fatigue due to the temperature difference between the lower and upper bypass (60 and 276^◦C respectively) was thought to be the source of the problem. A extensive study carried out by Forsmarks Kraftgrupp AB (Tinoco and Lindqvist (2009)) found that the origin of the cracks were low frequency temperature oscillations of about 0.2 Hz in the mixing region.

This master thesis was done at the Nuclear Reactor Technology department of KTH, Stockholm. It is part of the THEMFE project (Thermal Mixing and Fatigue Experiment) initiated by SSM, Forsmarks Kraftgrupp AB, Ringhals AB, Oskar- shamns Kraftgrupp AB and Teollisuuden Voima Oyj (TVO). The aim of THEMFE is to obtain experimental data to predict heat transfer.

The project continues the work done by Reijo Pegonen at his master thesis in 2012 (Pegonen (2012)). The objective is to simulate the thermal mixing and analyze how the low frequency temperature oscillations affect the control rod stem. Comparing with Reijo’s work, the main points of innovation are: Fluid and solid parts are all included in the CFD simulation, temperature dependent water properties c_p, k and µ are used, the turbulence models for RANS are Low Reynolds, the cell size needed for LES is theoretically derived from Kolmogorov’s energy spectrum.

The case setup was carried out using Solid Edge for the geometry, ANSYS-ICEM for the mesh and OpenFOAM for the CFD calculations.

The geometry of the control rod and casing was taken from the experimental setup at KTH. The thermal mixing was studied between two streams at different temperatures and max fluxes: Cold inlet at 60^◦C and 0.07kg/s and hot inlet at 276^◦C and 0.8kg/s, both of them at 72bar. For RANS the Low Reynolds k-Omega SST tur-

1

(20)

bulence model was chosen for its accuracy at the near wall region, for LES we used Smagorinsky sub grid scale model, not the Dynamic version, as this one was not properly implemented in OpenFOAM 2.2. The results were analyzed with the data given by the probes placed at the problematic areas where cracks were suspected to appear.

These results will be useful for the experimental setup of the THEMFE project in order to know where to place the thermocouples and for the PhD student Roman Thiele (NRT department, KTH, Stockholm), who will continue improving our work.

(21)

Chapter 2

Theoretical Background

This section shows the fluid equations that we used and gives an explanation of why we used them.

As in the simulation, water is considered to be an incompressible fluid.

Special care was taken with the turbulence models, the behavior close to the wall for the URANS and LES simulations and the Kolmogorov energy spectrum for the turbulent eddies.

2.1 Navier - Stokes Equations

Navier-Stokes equations of continuity and momentum for an incompressible fluid make up a system with 4 equations and 4 unknowns: 3 velocity components (u, v, w) and the pressure p.

∂ui

∂x_i = 0 (2.1)

ρ∂ui

∂t + ρuj

∂uj

∂xj

= − ∂p

∂xj

+ ∂

∂xj

[2µSij] + ρgi (2.2)

Strain rate tensor Sij =¹₂_∂u

i

∂xj +^∂u_∂x^j

i

The temperature field T is calculated using the energy equation (2.3) once the velocity field is known.

ρD Dt

h +1

2u_iu_i

=∂p

∂t + ∂

∂xi

[u_iS_ij] −∂q⁰⁰_i

∂xi

+ ρg_iu_i (2.3)

D

Dt represents _∂t^∂ + ui ∂

∂xI

Enthalpy h = cpT 3

(22)

Heat flux q_i⁰⁰using Fourier’s law becomes −λ_∂x^∂T

i

The system of equations is complete, it seems that we could directly discretize the equations and calculate the fluid properties at every grid point.

But solving Navier-Stokes equations numerically is not a simple task, it would require a very fine grid to resolve all the turbulent scales and a fine resolution in time (turbulent flow is always unsteady). This what DNS codes do in a computational time that could take many months for this particular simulation. As we are not able to wait this long, another approach must be taken.

2.2 Reynolds Averaging

RANS (Reynold-Averaged Navier Stokes) simplifies the formulation separating the variables into a mean and a fluctuating property in order to make turbulence sta- tionary (Wilcox (2010), section 2.1).

The decomposition of the velocity would be as:

u(x, t) = u(x) + u⁰(x, t) = U + u⁰ (2.4) Where the average value calculated from:

u(x) = lim 1 T

Z

u(x, t)dt (2.5)

As the fluctuating value oscillates around the mean, we can establish the following properties for any functions φ and ψ

φ⁰ = 0 φ + ψ = φ + ψ φψ = φψ + φ⁰ψ⁰ ∂φ

∂t =∂φ

∂t

When Navier-Stokes equations are averaged at both sides, the continuity (2.1) and momentum (2.2) equations become:

∂Ui

∂x_i = 0 (2.6)

ρ∂Ui

∂t + ρUj

∂Uj

∂xj

= − ∂p

∂xj

+ ∂

∂xj

µ ∂U_i

∂xj

+∂Uj

∂xi

− ρu⁰_iu⁰_j

+ ρgi (2.7)

τ_ij^R= −u⁰_iu⁰_j (2.8)

We can see that we have just created 6 more unknowns u⁰_iu⁰_j in symmetric tensor τ_ij^R called Reynolds Stress Tensor. The next step is to find a solution for τ_ij^R. This is where turbulence modeling starts.

(23)

2.2. REYNOLDS AVERAGING 5

2.2.1 Turbulence Modeling

As we saw in the last section, the system of equations needs to be closed. Prandtl introduced the turbulent kinetic energy in 1945, the kinetic energy contained in the velocity fluctuations, as the basic variable for turbulence modeling.

• Specific turbulent kinetic energy k:

k = 1

2u⁰_iu⁰_i (2.9)

We can see from equations (2.8) and (2.9) that k is proportional to the trace of the Reynolds-Stress Tensor. This means that if k is calculated we reduce the number of unknowns by 3.

τ_ii^R= −u⁰_iu⁰_i= −2k (2.10) The developing of the k equation follows a complicated process that can be read at Wilcox (2010) section 2.4. In essence, it introduces relation (2.10) in the Reynolds- Stress equation.

∂k

∂t + U_j ∂k

∂x_j = τ_ij^R∂Ui

∂x_j − + ∂

∂x_j

ν∂k

∂t −1

2u⁰_iu⁰_iu⁰_j−1 ρp⁰u⁰_j

(2.11) The term in Equation (2.11) is called dissipation rate. A turbulent model that has as a variable in the k equation is known as k-Epsilon turbulent model. If is replaced by ω, the specific dissipation rate, the turbulent model is called k-Omega.

• Dissipation rate : Rate at which turbulent kinetic energy is converted into thermal internal energy.

 = ν∂u⁰_i

∂xk

∂u⁰_i

∂xk

(2.12)

• Specific dissipation rate ω: Rate of dissipation of energy in unit volume and time. The variable l is the turbulent length scale.

ω = ck^1/2

l (2.13)

It seems that the system of equations is even bigger after the introduction of the k equation. The last two tensors on the left side of equation (2.11) just gener- ated another 16 new unknowns. These unknowns were just replaced by closure approximations, which are as follows.

(24)

2.2.2 Reynolds - Stress Tensor

Notice that in equation (2.11) we still have the Reynolds stress tensor τ_ij^R. Boussi- nesq stated that this tensor can be approximated to the product of the turbulent eddy viscosity and the mean strain-rate tensor (Wilcox (2010), section 4.1).

τ_ij^R= 2νtSij−2

3kδij, Sij = 1 2

 ∂u_i

∂xj

+∂uj

∂xi

(2.14) Since Sii= 0 for incompressible flow and δii= δ11+ δ22+ δ33= 3; equation (2.14) is in accordance with (2.9)

Using dimensional arguments, the momentum transfer by turbulent eddies, the turbulent eddy viscosity ν_t, can be written as:

νt= Cµk^1/2/l (2.15)

Where l is a physical quantity describing the size of the large energy-containing eddies and Cµa closure constant.

2.2.3 Turbulent Transport and Pressure Difussion

The turbulent transport term, which can be regarded as the turbulent energy trans- ported through the fluid by turbulent fluctuations, and the pressure diffusion term, which characterizes the turbulent transport by pressure fluctuations (most of the time this term is so small that can be neglected), are approximated to (Wilcox (2010), section 4.1):

1

2u⁰_iu⁰_iu⁰_j+1

ρp⁰u⁰_j= −ν_t σk

∂k

∂xj

(2.16) Where σk is a closure coefficient.

2.2.4 System of Equations

The result of the two previous approximations gives us a new version of the k equation (2.11)

∂k

∂t + Uj

∂k

∂x_j = τ_ij^R∂Ui

∂x_j − + ∂

∂x_j

ν + νt

σ_k

∂k

∂x_j

(2.17) We now have a system composed by equations (2.6), (2.7), (2.14) and (2.17); 10 equations in total

With 13 unknown variables:

• Three velocity components (u, v, w)

• Pressure p

• 6 components of the Reynolds stress tensor τ_ij

(25)

• Turbulent kinetic energy k, dissipation and turbulent length scale l One more equation for epsilon or omega and we would have closed the system.

2.2.5 k-Epsilon Turbulence Model

The idea of this model is to derive the exact equation for . Using the definition of

 presented at equation (2.12) Chou arrived to a formula in 1954 which after big simplification turns to appear as:

∂

∂t+ Uj

∂

∂xj

= C1

kτij

∂U_i

∂xj

− C2

² k + ∂

∂xj

ν + ν_t σ

∂

∂xj

(2.18) The turbulent eddy viscocity is calculated as:

ν_t= C_µk²

(2.19)

The corresponding closure coefficients are:

Table 2.1: Closure coefficients for the Launder and Sharma k − model

C₁ C₂ C_µ σ_k σ 1.44 1.92 0.09 1 1.3

This k-Epsilon model preforms reasonably well for 2D and free shear flows, but the results are quite inaccurate when dealing with 3D flows which involve impinging jets or adverse pressure gradients at the boundary layer (Pope (2000), section 11.10).

The model also fails to predict a satisfactory value of the constant C in the law of the wall and cannot be integrated through the viscous sub-layer.

2.2.6 k-Omega Turbulence Model

This model was first postulated by Kolmogorov in 1942. At first, he derived an equation for ω combining physical reasoning and dimensional analysis. After some improvements the model became as follows (Wilcox (2010), section 4.3.1):

∂k

∂t + Uj

∂k

∂xj

= τ_ij^R∂U_i

∂xj

− βkω + ∂

∂xj

ν +ν_t

ω

∂k

∂xj

(2.20)

∂ω

∂t + Uj

∂ω

∂x_j = −αω kτ_ij^R∂Ui

∂x_j − βkω +σd

ω

∂k

∂x_j

∂ω

∂x_j + ∂

∂x_j

νt+ σk ω

∂ω

∂x_j

(2.21) Turbulent eddy viscosity is calculated as:

(26)

νt= k

eω, ω = maxe (

ω, 7 8

s 2SijSij

β^∗ )

(2.22) The corresponding closure coefficients are:

Table 2.2: Closure coefficients for the k − ω model

α β β^∗ σ σ^∗ σdo σd

13/25 Function 9/100 1/2 3/5 1/8 Function

The k-Omega model can be directly used to treat the near wall viscous sub-low layer.

It also accounts for effects in of the stream wise pressure gradient and provides good accuracy for free shear and separated flows. However, it does not work well at free stream conditions and it is very sensitive to a change in the inlet value ω (Pope (2000), section 10.5.1 & ANSYS (2012), section 2.2.2.4.1).

2.2.7 k-Omega SST Turbulence Model

Developed by Menter in 1994. This model combines the strengths of the previously explained turbulence models: It works as k-Omega in the near wall region and as k-Epsilon in the free shear flow (ANSYS (2012), section 2.2.2.5).

The equation for the k is the same as in the k-Epsilon model, the ω equation is written as:

Uj

∂ω

∂xj

= γ 1 νT

Pk− βω²+ ∂

∂xj

(ν + σωνt) ∂ω

∂xj

+ 2(1 − F1)σω2

1 ω

∂k

∂xj

∂ω

∂xj

(2.23) Turbulent eddy viscosity is calculated as:

ν_t= k

ω (2.24)

F₁ is a blending function that switches from one model to the other.

This model is recommended for high accuracy boundary layer simulations.

2.2.8 Behaviour Close to the Wall

The previous turbulent models are called High Reynolds models, they state equations that can only be applied at fully turbulent conditions. If we use them we will not take Low Reynolds effects into account, such as the sharp peak of the turbulent kinetic energy close to the wall.

There are two ways of dealing with this:

(27)

Wall functions

Low Reynolds number effects are not computed. A cell at the wall boundary is taken thick enough to ensure that we can assume fully developed turbulent boundary layer conditions (30 < y⁺ < 200 for the first cell centroid). This first grid point of the cell is known as “matching point” and it uses the law of the wall with a suitable value of C to relate the velocity field to the surface shear stress and to achieve a proper matching between the viscous (Eq.(16)) and logarithmic velocity profiles (Eq.(17)) (Wilcox (2010), section 4.7.1).

U⁺= y⁺ (2.25)

U⁺= 1

kln y⁺ + C (2.26)

The dimensionless parameters U⁺and y⁺ are defined as _u^U

τ and ^u^τ_ν^y respectively.

Once the friction velocity uτ is calculated we use the wall functions to compute the values for k, or ω.

k = u²_τ

√β^∗, ω = k^1/2

β^∗1/4κy,  = β^∗3/4k^3/2

κy (2.27)

Where β^∗ is a constant and κ the Von Karman constant.

Keep in mind that this procedure is only applied to the first cell near the wall. The calculated values of k, and ω are used as boundary conditions.

Wall functions are only an approximation, they work quite well simulating internal flows but the results would not be very accurate if they were applied to external flows or heat transfer calculations. Furthermore, they are very sensitive to the location of the first grid point near the wall as it can be seen in the inverse dependence on the coordinate y at the equations above.

Low-Re turbulence models

Low Reynolds number effects are resolved. This new approach is based on inserting enough number of grid lines close to the wall (y⁺= 1 for the first cell centroid) so that the boundary layer can be adequately resolved. To do so, some modifications have to be made on the equations:

High Reynolds turbulent models predict that for y → 0 (Wilcox (2010), section 4.6.3)

k ∼ yⁿ (2.28)

βy²ω/ν = m (2.29)

Applying the non-slip boundary condition and the continuity equation close to the wall, we can approximate the order of magnitude of the fluctuating velocities u⁰,

(28)

v⁰and w⁰. Then, using the k equation (2.17) we get that n = m = 2 (Wilcox (2010), section 4.6.3).

k-Epsilong and k-Omega turbulence models fail to predict both of these parameters.

• Launder-Sharma k-Epsilon: n = 1.39 and m = 0.53

• Wilcox k-Omega: n = 3.23 and m = 7.63

The Low Reynolds models include a set of damping functions fµ, f1, f2, 0, E with some modified coefficients C1, C2, Cµ, σk , σ whose purpose is to make a correction on then and m values. The resulting modified equations for k and can be seen at Wilcox (2010), section 4.9.1. Notice that the variable for the dissipation rate is no longer bute. is calculated as +e, where the damping function 0is the value of at y → 0, defined differently for each particular model.

The eddy viscosity is defined as:

ν_t= C_µf_µk²

(2.30)

2.3 Large Eddy Simulation

Hinze (1975) wrote:

“Turbulent fluid motion is an irregular condition of flow in which the various quantities show a random variation with time and space coordi- nates”

This turbulent process starts by the creation of large scale turbulent structures (large eddies) that release their energy transferring it to smaller eddies. The energy cascade continues until the Reynolds number is sufficiently small that the eddy motion is stable and molecular viscosity is effective in dissipating the kinetic energy.

Large eddy simulations simulate the large eddies and model the small ones. The possibility to simulate eddies is what makes LES a much powerful tool than RANS, never less, requiring more computational power.

2.3.1 Filtering

The filtering operation decomposes the velocity into a sum of filtered or resolved component and a residual or sub-grid scale, SGS, component (Wilcox (2010), section 8.3.1).

ui= ui+ u⁰_i (2.31)

This is done by applying a filtering function G:

ui(x, t) = Z Z Z

G(x − ξ; 4)ui(ξ, t)d³ξ (2.32)

(29)

2.3. LARGE EDDY SIMULATION 11

Where 4 = (4x4y4z)^1/3 and G is expressed as:

G(x − ξ; 4) =

( 1/4³, |xi− ξi| < 4xi/2

0 otherwise (2.33)

Equation (2.32) will filter out scales smaller than the mesh size, the mesh size is the cut-off length l for the smallest eddie simulated.

2.3.2 Simulated Eddies

When Navier-Stokes equations (2.1) and (2.2) are filtered we get:

∂ui

∂xi

= 0 (2.34)

ρ∂ui

∂t + ρ∂(ujui)

∂xj

= − ∂p

∂xj

+ ∂

∂xj

µ ∂ui

∂xj

+∂uj

∂xi

(2.35) The nonlinear filtered term u_ju_i will be the main problem in LES modeling. The deomposition of this term into a fluctuating and filtered variables is (Wilcox (2010), section 8.3.1):

uiuj = uiuj+ uiu⁰_j+ u⁰_iuj+ u⁰_iu⁰_j (2.36) The terms in the left side of equation (2.36) can be physically regarded as the interactions between the large scale filtered eddies u and the sub grid scale modeled eddies u⁰. These terms are separated into three variables known as:

Leonard Stress L_ij= u_iu_j− u_iu_j

Cross Term Stress Cij = uiu⁰_j+ u⁰_iu_j

SGS Reynolds Stress Rij = u⁰_iu⁰_j

The summation of the three terms is known as the SGS stress tensor τ_ij^SGS = L_ij+ R_ij+ C_ij = u_iu_j− u_iu_j. With this, equation (2.35) is rewritten to:

ρ∂u_i

∂t + ρ∂(u_iu_j)

∂xj

= − ∂p

∂xj

+ ∂

∂xj

µ ∂u_i

∂xj

+∂u_j

∂xi

+ τ_ij^SGS

(2.37)

To close the system of equations we need a turbulence model that can compute τ_ij^SGS.

(30)

2.3.3 Sub-Grid Scale Models

Smagorinsky (1963) postulated a model based on the Boussinesq approximation for the Reyndols srtess tensor (2.14):

τ_ij^SGS = 2νtSij, Sij =1 2

 ∂u_i

∂xj

+∂uj

∂xi

(2.38) Here, the turbulent eddy viscosity is no longer calculated as a combination of k,  and ω as we saw in equations (2.19) and (2.22). Instead, it is directly computed form:

ν_t= (C_s4)²pS_ijS_ij (2.39) Where C_s is the Smagorinsky coefficient that can be chosen between 0.1and 0.24 and 4 the grid scale (4_x4_y4_z)^1/3.

The problem is that if we use a fixed value of C_s the method will predict a non- zero sub-grid eddy viscosity ν^SGS close to the wall where the flow is laminar.

To overcome this problem, an improved version of this method called Dynamic Smagorinsky was postulated by Germano in 1997. In this model the Smagorinsky constant Csis no longer a fixed value; instead, it is calculated at every cell for each time step. The process it follows is (Wilcox (2010), section 8.3.2):

• The resolved stress tensor τij is used to compute a value for Cs.

• This C_svalue is used for the sub-gird scale tensor τ_ij^SGSin the next time step.

Resulting in a better performance at the boundary layer.

2.3.4 Kolmogorov Energy Spectrum

The eddy length scale in the energy cascade varies in a wide range, Kolmogorov -5/3 Law relates the spectral distribution of the energy E in terms of the wave number κ = 2π/l, equation (2.40) (Pope (2000), section 6.5.3).

Kolmogorov’s hypothesis was:

“In a turbulent flow at sufficiently high Reynolds number, the statistics of the small scale motions have a universal form that is uniquely defined by ν and ”.

E(κ) = C^2/3κ^−5/3fLfη (2.40) The variable l represents the eddy length scale and the dissipation rate, which in a k-Omega model can be computed as ₁₀₀⁹ kω.

The production range is governed by fL:

f_L=







κl₀ h

(κl0)²+ cL

i1/2







p₀+5/3

(2.41)

(31)

2.3. LARGE EDDY SIMULATION 13

The dissipation range is governed by fη: f_η= exp(−β

h

(κη)⁴+ c⁴_Li1/4

− c_η

(2.42) The length scales are calculated as:

Large eddies l₀=

ν³

1/4

Re^3/4_L

Small eddies (Kolmogorov length scale) η =

ν³

1/4

Reynolds number associated with large eddies:

ReL =k^1/2l₀

ν (2.43)

Constant values used:

Table 2.3: Constants used in Kolmogorov -5/3 Law

cL cη C p0 β 6.78 0.4 1.5 2 5.2

An example of the Kolmogorov spectrum is displayed in Figure 2.1.

10⁻³ 10⁻² 10⁻¹ 10⁰ 10¹

10⁻¹⁵ 10⁻¹⁰ 10⁻⁵ 10⁰ 10⁵

k*eta

E/(eta*u2)

Kolmogorov energy spectrum

Figure 2.1: Normalized Kolmogorov energy spectrum

This spectrum is divided into three zones:

(32)

• Bulk of the energy 6l₀> l > ¹₆l₀: Where the most energetic eddies are.

• Inertial range (¹₆l0> l > 60η): Kinetic energy is transferred to smaller eddies,

.

• Dissipation range (60η > l > η): Kinetic energy is dissipated by viscosity, ν.

A good LES simulation should resolve the eddies that contain 80% of the energy and model the remaining 20% (Pope (2000), section 13.1). Kolmogorov equations will be very useful to determine the mesh size needed to achieve this.

2.3.5 Behaviour Close to the Wall

In the last paragraph of the previous section we said that we should simulate 80% of the energy containing eddies of the flow to have a proper LES. But what happens at the boundary layer? Below y⁺ = 20 the turbulent kinetic energy and other terms have a sharp peak, meaning that if we want to solve 80% of the energy we would need to have a very fine mesh resolution which will make the simulation as expensive as a DNS (Pope (2000), section 13.4.5).

To reduce the computational cost, instead of resolving 80% of the energy everywhere we will just resolve 80% of it in the remote part of the wall and the boundary layer will be modeled. This is called LES-NWM, LES with Near Wall Modeling.

LES-NWM uses wall functions that act as a boundary condition at the first cell point. We used Spalding’s wall function, which provides a turbulent viscosity ν_t value based on the u⁺. The main advantage of this formula is that it can be applied at any y⁺ value at the boundary layer (de Villiers (2006), section 3.3).

y⁺= u⁺+ 1 E

exp(κu⁺) − 1 − (κu⁺) −1

2(κu⁺)²−1 6(κu⁺)³

(2.44)

Von Karmann Constant κ = 0.41 Thermal diffusivity αt is computed from νt

2.4 Buoyancy

Apart from the turbulence modeling in URANS and the sub-grid scales in LES there is still one section that has to be taken into account for a proper CFD calculation.

This section is the buoyancy, the force experienced in the fluid when density changes happen.

There are two ways of including buoyancy in our calculations:

Boussinesq Approximation

In a fluid where the temperature changes are not so big (2^◦C for water), density can be treated as a constant in every term of the momentum equation excepting the gravitational term where it will be treated as variable. When temperature

(33)

2.4. BUOYANCY 15

differences are larger, the use of the Boussinesq approximation can lead to big errors (Ferziger and Peric (2001), section 1.7.5).

The modification of the gravity term would be:

ρgi= ρ0gi[1 − β(T − T0)] (2.45) Where T0 is a reference temperature, ρ0 the density at T0 and β the thermal expansion coefficient.

Compressible flow

Using compressible fluid flow equations with a temperature dependent density will automatically compute the buoyancy forces, no need to use any approximation here.

(34)

(35)

Chapter 3

CDF tool OpenFOAM

OpenFOAM is a open source CFD software written in C++. In this project, we used OpenFOAM version 2.2.

OpenFOAM uses finite volume discretization schemes. The finite volume method uses the integral form of the conservation equations to calculate the properties at centroid of the cells. Interpolation between the centroids is used to calculate fluxes at the cell boundaries, the faces (Ferziger and Peric (2001), section 4.2).

As an example, the finite volume form of the contnuity equation (2.1) is solved by stating that the net mass flux through the cell faces must be equal to zero.

3.1 Solver

The only available solver to simulate a multi region heat transfer problem is chtMultiRegionFoam.

This solver can run steady state calculations using a SIMPLE loop and transient doing PIMPLE. As we are dealing with a transient problem, I will only explain how PIMPLE works.

3.1.1 File Structure

The solver uses different files to solve the equations. The order in which they care called is:

chtMultiRegionFoam.c











1 SolveFluid.H







1.1 UEqn.H → Momentum Eq (2.7).

1.2 EEqn.H →Energy Eq (2.3) 1.3 pEqn.H → Pressure Eq.

2 SolveSolid.h → Conduction Eq.

3.1.2 Calculation Method

There is no independent pressure equation in Navier-Stokes, this means that it is difficult to find a pressure field that satisfies the continuity equation. To overcome

17

(36)

this problem, OpenFOAM and most of the CFD solvers take these steps to compute U and p (Ferziger and Peric (2001), section 3.1).

1. Solve the momentum equation to compute an intermediate velocity field U^∗: This equation does not take pressure into account.

∂U_i^∗

∂t + U∂U_j^∗

∂xj

= ∂

∂xj

νef f

 ∂U_i^∗

∂xj

+∂U_j^∗

∂xi

+ gi (3.1)

2. Compute the max fluxes at the cell faces: Needed for the finite volume discretization.

3. Solve the pressure equation: Using the previously non used terms of the momentum equation.

∂Ui

∂t + ∂p

∂xj

= 0 (3.2)

Taking the divergence of equation (3.2):

∂

∂xi

 ∂U_i

∂t + ∂p

∂xj

= 0 (3.3)

The continuity equation states that that ^∂U

t+4t i

∂xi = 0 → _∂x^∂

i

h_U^t+4t

i

4t

i

= 0:

The unknown U_i^t+4t is simplified and p becomes the only unknown in the laplacian equation (3.3).

4. Correct the mass fluxes at the cell faces.

5. Use p and U^∗ to solve equation (3.2) and compute the real velocity field U_i^t+4t.

6. Update the boundary conditions.

7. Repeat the process from step 2 as many times as the user specified in the nOuterCoor parameter for the PISO loop.

chtMultiRegionFoam solver does not use pressure p to do the calculations but a modified prgh = p − ρgh. This is done to include the buoyancy term ρgi in the pressure correction shown at step 2.

As we explained before, pressure is hard to converge and it needs more iterations than the velocity or energy. The user can define the number of iterations via these two parameters:

• nNonOrthCorr: Inner loops for pressure accounting for non-orthogonalities of the mesh.

• nCorr: Outer correction loops for pressure.

(37)

3.2. CASE STRUCTURE 19

3.2 Case Structure

The basic file structure of a non-multi region case in OpenFOAM is:











Init & Bound Cond 0









 U T p p_rgh k omega mut

Constants constant











M esh polyMesh

T urbulence model turbulenceProperties ρ, cp, k, µ thermoPhysicalProperties Bouyancy transportProperties

N umerics system







N umercal schemes fvSchemes

Solvers fvSolution

Run control controlDict

In a multi region case every region needs to its own parameters defined. The directories 0, constant and system include as many folders as regions and the structure changes to:

(38)









 0









 WATER

n

U , T , p, , p_rghk, ω, , ν_t INNERCYL

n T OUTERCYL

n

”

constant









 WATER











polyMesh

turbulenceProperties thermoPhysicalProperties transportProperties INNERCYL

(

polyMesh

thermoPhysicalProperties OUTERCYL

n

” polyMesh

regionProperties

system









 WATER

(

fvShcemes fvSolution INNERCYL

(

fvSchemes fvSolution OUTERCYL

n

” controlDict

In our case, the boundary and initial conditions were not specified directly in the 0 file. Instead, we included a changeDictionaryDict file in system > regionName that rewrote the files at 0 before running the simulation.

We also used many other OpenFOAM utilities such as topoSet, createPatch, mapFields, mirrorMesh, etc. All these utilities need a dictionary named utilityDict which is usually located in the system folder, we did not include them in the graph above for simplicity.

(39)

Chapter 4

Choosing the Turbulence Model

In order to fin the most suitable turbulence model for our 3D simulations we compared the results of two different Low Reynolds turbulence models (Launder- Sharma k-Epsion and k-Omega SST) to the experimental data that Kang and Pratil published in the Inernational Journal of Heat and Mass Transfer (Kang and Patil (2000)).

The physical model that they used was an annular flow channel, similar to our 3D case. In order to reduce the computational time we made a quasi 2D model.

4.0.1 2D Model

The annular flow channel is simplified to a wedge of 2.5^◦. The geometry and mesh were built in OpenFOAM using blockMesh utility (Figure 4.1 and Table 4.1).

Figure 4.1: Mesh of the simplified 2D model

21

(40)

Table 4.1: Geometrical parameters of the 2D model

Water Outer cylinder

Length [mm] 3000

φ [mm] 38.1/15.8 41.2/38.1

4.0.2 Validation

In the experimental setup used by Kand and Pratil water is pumped through an annular flow channel at different Reyndols numbers, temperatures and thermal conditions. The measured variables are the radial velocity an temperature profiles.

We choose to validate experiment number 8:

• Flow parameters: Re = 16000, P r = 7.16

• Water inlet: T = 42.2^◦C P = 269kP a

• Heat flux from the outer radius of the outer cylinder: q⁰⁰ = 9000W/m² The OpenFOAM simulations were run with the chtMultiRegionFoam solver in a transient PIMPLE loop until convergence. The two turbulence models which were validated are shown in Figures 4.2 and 4.3. In these figures, the red dots represent the experimental results, we did not include error bars because they were even smaller than the dot.

0 0.2 0.4 0.6 0.8 1

40 42 44 46 48 50 52 54 56 58 60

Dimensionless radius R*

T [C]

Temperature profile Experimental OpenFOAM

0 0.2 0.4 0.6 0.8 1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

u [m/s]

Velocity profile

Experimental OpenFOAM

Figure 4.2: Low Reynolds LS k-Epsilon turbulence model validation. Left: Temperature profile.

Right: Velocity profile

(41)

23

Both models show good accuracy. In one hand, LS k-Epsilon works better at the center of the pipe, whereas k-Omega SST has more accuracy with both temperature and velocity in the near wall region. Table 4.2 shows the mean errors of the two validations. These errors were calculated using equation (4.2)

M ean error =P errorP oint

NP oints

100 (4.1)

k-Omega SST shows a better performance in the near wall region, the mean errors for temperature and velocity are also lower than LS k-Epsilon. As the near wall region is the region of study in this project we will use k-Omega SST in the next 3D simulations.

0 0.2 0.4 0.6 0.8 1

42 44 46 48 50 52 54 56 58 60

T [C]

Temperature profile Experimental OpenFOAM

0 0.2 0.4 0.6 0.8 1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

u [m/s]

Velocity profile

Experimental OpenFOAM

Figure 4.3: Low Reynolds k-Omega SST turbulence model validation. Left: Temperature profile, Right: Velocity profile

Table 4.2: Mean errors form the validation of LS kEpsilon and kOmegaSST Low Reynolds models

4T [%] U [%]

Low Re LS k-Epsilon 0.4978 2.95 Low Re k-OmegaSST 0.933 3

(42)

(43)

Chapter 5

Methodology

The three key elements in a CFD simulation are:

1. Mesh

2. Turbulence model

3. Numerical schemes

Once these three steps are properly implemented for our specific case we can set up the program to run the simulation. While running the simulation, the stability, speed and accuracy of the results will decide if steps 1, 2 and 3 which we thought correct were really ok.

The turbulence model has already been chosen to be k-Omega SST. In this section we will deal with the mesh and the numerical schemes.

5.0.3 Physical model

The model used for the CFD calculation is equal to the experimental setup for the THEMFE project at KTH university, Stockholm. A general view of this model is shown in Figure 5.1.

The control rod is represented by the inner cylinder, the casing by the outer cylinder.

The bottom and top of the pipe are walls, there is no water flowing through them.

25

(44)

Figure 5.1: 3D view of the physical model

The dimensions of the pipe are shown in Tables 5.1 and 5.2:

Table 5.1: Geometrical parameters of the physical model

Inner Cyl Outer Cyl Water

Length [mm] 1150

φ [mm] 35/25 88.9/80 80/35

Table 5.2: Geometrical parameters of the inlets and outlet

Cold Inlet Hot Inlet Outlet

Position [mm] 150 800 1000

φ [mm] 7.5 7.5 14

5.0.4 Geometry

Solid Edge version 10.4 was used to create the geometry.

As we are dealing with a symmetric body we will only generate geometry and mesh for one quarter of the pipe. This will improve the quality of the results and suppress the possibility of having eddies due to the asymmetry of the mesh.

The different regions involved in our CFD calculation are created in different Solid Edge ISO Part files:

• Parts: Outercyl, Innercyl, InletPipe and OutletPipe

Thermal Mixing CHT Simulations with OpenFOAM:URANS and LES