Juni 2019
Numerical Study on Hydrodynamic Characteristics of Flood Discharge
Tunnel in Zipingpu Water Conservancy Project
Using RANS equations and the VOF model Micaela Hamberg
Signe Dahlin
Teknisk- naturvetenskaplig fakultet UTH-enheten
Besöksadress:
Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0
Postadress:
Box 536 751 21 Uppsala
Telefon:
018 – 471 30 03
Telefax:
018 – 471 30 00
Hemsida:
http://www.teknat.uu.se/student
Numerical Study on Hydrodynamic Characteristics of Flood Discharge Tunnel in Zipingpu Water
Conservancy Project
Micaela Hamberg och Signe Dahlin
To avoid the large amount of damage that floods can cause, spillway tunnels are used to control water levels. To ensure the safety of water transportation through spillway tunnels, the behaviour of the water throughout the tunnel is important to know. Physical experiments are time consuming and expensive, hence CFD simulations are a profitable option for investigating the performance of the spillway tunnel. In this project, simulations of water flow in a spillway tunnel were executed.
A three dimensional model of the spillway tunnel in Zipingpu Water Conservancy Project was created in the software ANSYS Gambit. A coarse, middle and fine mesh with both hexahedral- and tetrahedral elements were also created for the model in ANSYS Gambit. The meshes were imported to ANSYS Fluent where the simulations, and a convergence analysis were made. The water flow was set to be described by the Reynolds-Averaged Navier-Stokes model, using the pressure solver, k-epsilon model and the VOF model. Physical experiments had previously been performed, and the simulated results were compared to these, in an attempt to find the parameters to replicate the experimental results to the greatest extent possible. The inlet velocity of the tunnel was known and the inlet boundary was set as a velocity inlet. The ceiling of the tunnel was set as a pressure inlet, the floor and walls were set as wall, and the outlet was set as pressure outlet. The simulated results showed similar behavior as the experimental results, but all differed from the experimental results. The grid convergence index, estimating the results' dependency on the mesh was 6.044 %. The flow was analyzed, and where the flow had unfavorable characteristics, such as a high cavitation number, the geometry of the spillway was altered in ANSYS Gambit to investigate if an improved geometry for the spillway tunnel could be found. The water flow in the revised geometry was simulated in ANSYS Fluent, and results showing flow with lower cavitation numbers was found.
Tryckt av: Uppsala
ISSN: 1401-5757, UPTEC F 19018 Examinator: Doc. Tomas Nyberg Ämnesgranskare: Dr. Per Norrlund
Handledare: Prof. James Yang, Prof. Yongliang Zhang och Prof. Ling Li
Populärvetenskaplig sammanfattning
På grund av Kinas geografiska placering och långa regnperioder påverkas landet ofta av översvämningar vilket bland annat kan orsaka stor skada på landets vattenkraftverk och omkringliggande städer. Således behövs tekniska lösningar för att leda bort överflödigt vatten och förebygga skada under översvämningar. En teknisk lösning som kan används i vattenreservarer är så kallade avrinningstunnlar. Designen av avrinningstunnlarna är avgörande för att upprätthålla säkerheten under översvämningar. Problem som kan uppstå och som kan leda till skador på avrinningstunneln är bland annat stora variationer av vattenytan, stor tryckvariation i tunneln samt att flödet blir för turbulent.
I detta projekt har syftet varit att undersöka vattenflödet genom en av vattenreservoaren Zipingpus avrinningstunnlar med hjälp av simuleringar samt undersöka geometrins påverkan på vattenflödet och därefter föreslå en förbättrad geometri. Då fysiska experi- ment ofta är svåra att implementera och kostsamma är en mer effektiv och ekonomisk lösning att använda sig av beräkningar. Ekvationerna som kan används för att erhålla en matematiskt korrekt beskrivning av ett vattenflöde blev härledda för mer än hundra år sedan och kallas för Navier-Stokes ekvationer. Trots att de har funnits i mer än ett århundrande har fortfarande ingen lyckats hitta en analytisk lösning till ekvationerna för majoriteten av vattenflöden. På grund av detta måste numeriska metoder och modeller ofta användas för att kunna beskriva majoriteten av vattenflöden korrekt utan allt för stora approximationer och felkällor.
I detta projekt har Reynolds-Averaged Navier-Stokes, RANS, ekvationer använts som approximerar Navier-Stokes ekvationer med ekvationer som enklare går att lösa nu- meriskt. När RANS ekvationer används måste även en turbulensmodell användas då RANS ekvationer inte är tillräckliga för att beskriva turbulens och i detta projekt har turbulensmodellen k-e används. Stegen i detta projekt har varit att först konstruera ge- ometrin av tunneln och därefter sätta upp en matematisk modell för att sedan genomföra beräkningar vars resultat jämfördes med experimentiellt uppmätt data. I den numer- sika metoden som använts i detta projekt har geometrin delats upp i många mindre volymer och i varje delvolym har beräkningarna genomförst. Oftast är det inte känt i första beräkningen hur många delvolymer geometrin behöver delas upp i, med för få delvolymer kan noggranheten bli väldigt liten samtidigt som för många delvolymer inte heller är att föredra då det ökar tiden det tar att genomföra beräkningarna markant.
Därför har beräkningar med tre olika storlekar på delvolymerna genomförts för att undersöka vilken storlek som är bäst lämpad. Därefter har tre olika variationer av
turbulensmodellen k-e (standard, RNG och Realizable) använts för att undersöka om de beräknade resulaten påverkas av valet av turbulensmodell. Tillsist ändrades tunnelns geometri för att undersöka hur det påverkar vattenflödet.
För alla simuleringar skiljde sig det resultat som erhölls från den experimentiella datan.
Dock följde de beräknade resultaten till stor del samma trend som den experimen- tiella datan. Störst likhet mellan det simulerade resultatet och den experminetiella datan erhölls då minsta storleken på delvolymerna användes och därför genomfördes resterande beräkningar med denna storlek (trots att simuleringen med denna storleken tog längst tid att genomföra). När de olika turbulensmodellerna jämfördes visade det sig att resultatet från simuleringar med standardmodellen stämde bäst överens med den experimentiella datan. Det visade sig även att vattenflödet genom den uppdaterade geometrin påvisade lägre mängd av vissa oönskade effekter än den ordinarie geometrin.
Numerical Study on Hydrodynamic Characteristics of Flood Discharge Tunnel in Zipingpu Water Conservancy Project
Vattenfall Uppsala University Tsinghua University
Department of Hydrualic Engineering
Approved by
Supervisor, Prof. Yongliang Zhang at Tsinghua University
Supervisor, Prof. Ling Li at Tsinghua University
Supervisor, Prof. James Yang at Vattenfall
Subject reviewer, Dr. Per Norrlund at Uppsala University
Examiner, Ass. Prof. Tomas Nyberg at Uppsala University
June 28, 2019
Acknowledgements
Our deepest gratitude goes to professor Yongliang Zhang for always guiding and en- couraging us before and throughout the project, and to professor Li Ling for helping us with any difficulties we encountered with the simulations.
Our sincere appreciation also goes to Li Zhengwen, Wenchuang Chen, Guo Peng and all team members in the office for ensuring that we had all help and possibility to do the project.
We also want to thank James Yang and Energiforsk for allowing us to do this project, complete our master degree and gain memories for life.
Table of Contents
Abstract II
Populärvetenskaplig sammanfattning III
Acknowledgements V
List of Acronyms VIII
1 Introduction 1
1.1 Background and aim . . . 1
1.2 Simulations . . . 2
2 Theory, the mathematical model and discharge tunnel 3 2.1 Mathematical description of fluids and fluid flows . . . 3
2.2 Numerical methods, CFD . . . 6
2.2.1 Discretization method . . . 6
2.2.2 Pressure or density based numerical solver . . . 9
2.2.3 Reynolds-Averaged Navier-Stokes Model . . . 10
2.2.4 Turbulence models . . . 12
2.2.5 Boundary conditions . . . 13
2.2.6 Multiphase flow modeling . . . 14
2.2.7 Mesh . . . 16
2.3 Discharge tunnel . . . 18
2.3.1 Design of spillway tunnel . . . 20
3 Methodology 21 3.1 Work process . . . 21
3.2 Building geometry . . . 21
3.3 Creating mesh . . . 24
3.4 Setting boundary conditions . . . 25
3.5 Establish solver model . . . 25
3.6 Running simulation . . . 26
3.7 Post-processing of the calculations . . . 26
3.8 Previously measured data from physical experiments . . . 27
3.9 Convergence analysis . . . 27
3.10 Validation of model . . . 27
3.11 Revision of spillway tunnel geometry . . . 28
4 Numerical results 29
4.1 Geometry . . . 29
4.2 Mesh . . . 30
4.3 Running simulation . . . 32
4.4 Results from the convergence analysis . . . 35
4.5 Validation of model . . . 40
4.6 Results from revision of spillway tunnel geometry . . . 45
5 Analysis 59 5.1 Mesh and convergence analysis . . . 59
5.2 Comparison of turbulence models . . . 60
5.3 Geometry and revision . . . 60
5.4 Sources of error . . . 62
6 Conclusion and outlook 63 7 References 64 A Appendices 66 A.1 Experimental data . . . 66
A.2 List of Responsibilities . . . 72
List of Acronyms
CFD - Computational fluid dynamics CV - Control volume
DNS - Direct Numerical Simulation FDM - Finite Difference Method FEM - Finite Element Method FVM - Finite Volume Method GCI - Grid convergence index
RANS - Reynolds-Averaged Navier-Stokes RNG - Renormalization group
VOF - Volume of Fluid
1 Introduction
1.1 Background and aim
The influence of nature is indisputably eminent, and controlling it is practically im- possible. Heavy rains have caused immeasurable amount of damage by flooding, and China in particular has been severely affected by this due to its geographic and climatic conditions. Therefore, the ability to control flooding is crucial to study and improve1. A common way to control flooding is using water reservoirs.
The Zipingpu water reservoir is located in the upper part of the Minijang River, in the Sichuan province, China. The reservoir was built in 2006 and is mainly used for irrigation but also for power generation and flood control. The site contains four hydro-electric generators and has a total generating capacity of 760 MW2. In a hydraulic reservoir, flood control is needed to prevent damage of the reservoir and the downstream riverbed during flooding. To enable flood control in a reservoir, flood discharge tunnels can be used. By using discharge tunnels large amounts of water can be discharged through the reservoir during flooding, and thus ensure safety. The water flow through the tunnel is affected by the geometry of the tunnel and therefore the geometry needs to be considered carefully3.
In this project the hydraulic characteristics of a flood discharge tunnel in the Zipingpu Water reservoir were studied by numerical simulations of water flow through the tunnel.
This was done by establishing a mathematical 3D model of partial free surface fluid flow and by simulations of the fluid flow using the computer program ANSYS Fluent. The velocity field, pressure distribution and position of the water surface were calculated.
The numerical results obtained were compared to previously experimental and pub- lished results. The main aim of this project was to study and design an optimized flood discharge tunnel. Other goals were to master both the engineering aspect of designing a flood discharge tunnel as well as computational fluid dynamics (CFD). The CFD parts included modelling, such as model set-up, turbulence models, meshing technique, model validation and result analysis.
1Tan Xuming, Wang Yinghua, Zhou Kuiyi, Flood Control and Management in China, Beijing: China Water- Power Press, 2005, p.19
2Xinhua News Agency, New Water Control Project Under Construction, China Through a Lens, 2002
3Charles C. S. Song, Fayi Zhou, Simluation of Free Surface Flow Over Spillway, Journal of Hydraulic Engineering, Vol. 125, 1999, p.959
1.2 Simulations
One of the most important things to consider when building a flood discharge tunnel, or any kind of hydraulic structure, is the safety. To ensure that the flood discharge tunnel is safe to use for operation during different circumstances hydraulic model tests can be done. However they are expensive and time consuming. A way to avoid this cumbersome problem is to use CFD technology to simulate the hydraulic problem4. If it is known that the results are reliable, this method can save both a lot of time, money and makes it possible to examine and decide between different designs before build- ing them5. To validate the model used in the CFD simulation, the results need to be compared to experimental data to ensure that the numerical methods and equations are applicable6. The tool used for the simulations in this project was ANSYS Fluent. ANSYS Fluent is a CFD software by ANSYS that is used to model and simulate various fields, such as turbulent flow and heat transfer7. ANSYS Fluent is based on mathematical equa- tions such as the Navier-Stokes equations and different turbulence models. Different boundary- and initial conditions can also be used.
It has been proven by experiments that Navier-Stokes equations can describe the flow of a Newtonian fluid, a fluid that obeys Newton’s law of viscosity, accurately. Since it is not possible, with some exceptions, to solve Navier-Stokes equations analytically, numerical methods must be used and these methods are called CFD8. Fluid flows are described by several partial differential equations and in order to obtain an approximate solution, by a numerical method, the partial differential equations must be approximated by a system of algebraic equations that can be solved by a computer. There exist many different approaches to approximate the differential equations and some of the most commonly used are; Finite element -, Finite volume- and the Finite difference methods9.
4Nils R. B. Olsen, Hilde M. Kjellesvig, Three-dimensional numerical flow modelling for estimation of spillway capacity, Journal of Hydraulic Research, Vol. 36, 1998, p.775
5Song, p.959
6Birjan Dargahi, Experimental Study and 3D Numerical Smulations for a Free-Overflow Spillway, Journal of Hydraulic Engineering, Vol. 132, 2006, p.899
7ANSYS, Inc, ANSYS Fluent, ANSYS, 2019
8Joel H. Ferziger, Milovan Peric, Computational Methods for Fluid Dynamics, 3 ed, Berlin: Springer, 2002, p.12
9Ferziger, p.23
2 Theory, the mathematical model and discharge tunnel
2.1 Mathematical description of fluids and fluid flows
Even though the basic equations of fluid dynamics have been known for over a century several of the fundamentals in this field are yet to be elucidated1. To be able to describe the theory of fluid dynamics, it is important to start by explaining what a fluid is. Every- body has a intuition of what a fluid is and how it behaves however the mathematical description is complicated.
A fluid is a substance that has no resistance to external shear forces, meaning a fluid particle will be influenced and deformed even by the smallest force. A fluid appears to be continuous on a macroscopic level but on a microscopic level there exists discontinu- ities, however it is convenient, in most cases, to treat a fluid as a continuous medium.
A continuous medium is a substance whose sub-parts are evenly distributed in the space it occupies and that can be divided into infinitesimal elements that still have the same properties as the original volume had. A fluid fulfills the above statements to some extent on a macroscopic level and can therefore be approximated as a continuous substance when the macroscopic behaviour of the fluid is of interest. The appearance of fluids are not always the same. Both liquids and gases are considered to be fluids and obey the same laws of motion and act similarly under influence of external forces, even though there is a significant difference between them. Macroscopic properties of fluids, such as density or viscosity, can differ a lot between different substances2. For instance, a liquid is not easily compressible meaning that the density of a liquid is not strongly depending on pressure and temperature. A gas on the other hand is highly compressible3.
Viscosity can be described as a fluid’s resistance to flow. For example it is often said that syrup is a viscous fluid while water is not. However, this expression is not always correct since the viscosity actually is not a property of a fluid but rather a property of a fluid flow. Thus, the same fluid substance can sometimes behave more viscous and
1Claude Godréche, Paul Manneville, Hydrodynamics and nonlinear instabilities, New York: Cambridge University Press, 1998, p.25
2Ferziger, p.1
3O.Kochukhov, V. Pavelenko, L. Rosenqvist, Fluid Dynamics Compendium, 12 ed, Uppsala: Uppsala University department of physics and astronomy, 2017, p.9
sometimes less viscous depending on the situation. Dynamic or absolute viscosity is the ratio between the shear stress and the shear rate and is a measurement of the internal friction. The kinematic viscosity is the ratio between the dynamic viscosity and the density4.
Fluids can be classified into two subgroups; Newtonian- and non-Newtonian fluids. The classification of a Newtonian fluid is that it obeys Newton’s law of viscosity while a non-Newtonian fluid does not. The following equation is Newton’s law of viscosity for one dimensional flows,
τ=µdv
dy, (2.1)
where τ [N/m2] is the shear stress, µ [Ns/m2] is the dynamic viscosity and dv/dy [s−1] is the shear rate5. The behaviour of Newtonian fluids is important to study since most of the familiar fluids are Newtonian, for example air, water and oils. Non-Newtonian fluids are also an important aspect to study but with fewer applications. Some examples of non-Newtonian fluids are blood, honey and liquid plastics. The fluid simulated in this project is a Newtonian fluid whose flow can be described by Navier-Stokes equations which includes a time dependent continuity equation and a time-dependent momentum conservation equation. In the following equations the continuity equation and the momentum conservation equations are presented,
∂ρ
∂t + ∇(ρv) =0 (2.2)
and
ρdv dt ≡ρ
h∂v
∂t + (v· ∇)vi
= −∇p+ρg+µ∇2v+ (µ+b)∇(∇ ·v), (2.3) where ρ [kg/m3] is the density, v [m/s] is the velocity, t [s] is the time, p [N/m2] is the pressure, g [m/s2] is the external acceleration, µ [Ns/m2] is the dynamic viscosity (first coefficient of viscosity) and b [Ns/m2] is the the second coefficient of viscosity.
The value of the b is often negligible and approximated to zero even for compressible flows6. However for some specific flows the value of b can be of importance and then it is commonly approximated to - 23 µ. The Navier-Stokes equations are derived from the mass- and momentum conservation equations, however the derivation is above the scope of this thesis and thus it will not be included. It is not possible to obtain an analytic solution, with some exceptions, to the Navier-Stokes equations without any simplifications. In the limited number of cases where the Navier-Stokes equations can be solved analytically, it is for specific fluid flows where several terms in the equations are zero. For other fluid flows, several of the terms in the equations are minor and can be neglected, however these kinds of simplifications give rise to errors in the final
4Godréche, p.27
5Kochukhov, p.25
6Kochukhov, p.115
solution. Also, for some flows even the simplified equations cannot be solved analyti- cally. Because of this a numerical method must be used for the majority of cases in order to obtain a description of a fluid flow7. This will be discussed in more detail in section 2.2.
Some commonly used simplifications of fluid flows are incompressible flow and steady state flow. For many fluid flows the density of each fluid particle and the density of the unified fluid can be assumed to be constant. This approximation can be made for many liquids but also for some gases. The easiest way to determine if a fluid is incompressible or not is to compute its Mach number which is the ratio between the flow speed and the sound speed. If the Mach number is below 0.3 the fluid can be approximated as an incompressible fluid. When a fluid is incompressible the divergence of the velocity becomes zero,
∇ ·v=0, (2.4)
which simplifies the Navier-Stokes equations8. From equation 2.3 and 2.4 the following simplified Navier-Stokes equations, valid for incompressible fluids, can be derived,
ρdv dt ≡ρ
h∂v
∂t + (v· ∇)vi
= −∇p+ρg+µ∇2v. (2.5) A fluid is never truly incompressible, however in this project the fluid can be approx- imated as an incompressible fluid because the fluid velocity is much smaller than the speed of the sound. This equation can much more easily be solved, especially when other approximations, such as the steady state flow are valid as well. Steady flow means that the velocity of the flow is constant in time at each point the fluid occupies and can be written mathematically as,
∂v
∂t =0. (2.6)
When calculating a flow it is important to identify the flow as either laminar or turbulent.
Turbulence has baffled engineers for decades and still there is not a distinct separation between laminar and turbulent flows. Almost every flow encountered in real life is tur- bulent, for example everything from raising smoke from a cigarette and water spiraling downwards the pipe in a sink to air flow around cars and planes. Words often used to describe turbulence are random and chaotic and it is well known that small disturbances in a laminar flow can eventually lead to the flow becoming turbulent9. Some of the main characteristics of turbulence are that the flow appears to be highly unsteady and the behavior of the velocity as a function of time could be described as random. Furthermore turbulent flows contain rotational flow structures called eddies and the transportation of these structures in the flow is called vortex stretching. Large eddies are formed which then cause new smaller instabilities and eddies to form hence energy is transferred from larger eddies to smaller ones. This process is called dissipation which provides energy
7Ferziger, p.12
8Ferziger, p.12
9Ferziger, p.266
to maintain the turbulence behavior of the flow until the eddies become too small to sustain the process. It is more complicated to mathematically describe a turbulent flow than describing a laminar flow since turbulent flows contain larger fluctuations in both space and time10.
The fluid flow in this project is Newtonian, incompressible, turbulent and unsteady until an equilibrium is reached. The flow is a multiphase flow consisting of two phases, air and water, with a free surface.
2.2 Numerical methods, CFD
As stated earlier, in order to calculate a fluid flow it is often necessary to use a numerical method. There exist many different approaches and models that can be implemented and in this chapter the ones used in this project will be described thoroughly.
2.2.1 Discretization method
When choosing a numerical approach to solve the Navier-Stokes equations the partial derivatives are approximated by algebraic equations which can be solved by a math- ematical algorithm, using a computer. This step requires a discretization model. For a transient flow, unsteady flow, the equations must be discretized in both space and time. There exist several discretization models and the most commonly used are; Finite element-(FEM), Finite volume-(FVM), and Finite difference method (FDM)11. In this project the computer software ANSYS Fluent has been used and it uses the FVM ap- proach that discretizes the domain into control volumes, CVs, using meshing technique.
Then in each CV the momentum and continuity equations are integrated to construct algebraic equations for the unknown variables, i.e pressure, velocity etc12. In each CV the unsteady transport equation for a scalar φ is applied,
Z
V
∂ρφ
∂t dV+
I
ρφ~v·dA~ =
I Γφ∇φ·dA~ +
Z
VSφdV, (2.7)
where A is the surface area vector,~ Γφ is the diffusion coefficient and Sφ is the source term. When using the FVM, equation 2.7 can be written as,
∂ρφ
∂t +
Nf aces
∑
fρf~vfφf · ~Af =
Nf aces
∑
fΓφf∇φf · ~Af +SφV, (2.8)
where Nf aces is the number of faces of the given CV, φf is the value of φ on the face f , A~f is the area of the face, ∇φf is the gradient of the scalar at face f and V is the
10Jiyuan Tu, Guan Heng Yeoh, Chaoqun Liu, Computional fluid dynamics - a practical approach, 2 ed, Oxford: Elsevier Ltd, 2008 p.96
11Ferziger
12Ansys Inc, Overview of flow solvers, Ansys Fluent Documentation, 2009
volume of the CV13. By default, ANSYS Fluent stores the values of the scalar φ at the center of each control volume, however face values φf appear in equation 2.8 and must therefore be calculated separately. Fortunately they can be interpolated from the CV center values by using an upwind scheme. By upwind scheme it is implied that the face scalar value is calculated for quantities upstream relative to the fluid velocity, it is calculated from cell centred values from where the wave is propagating from. There are several upstream schemes available in ANSYS Fluent; first-order upwind, second-order upwind and QUICK are some among several.
The first order upwind scheme is the simplest scheme. When using the first-order upwind scheme only first order of accuracy is obtained and the scalar value at the face φf is approximated to the same value of the scalar value at the CV center φ. When higher order accuracy is required one option is to use the second order upwind scheme which yield second order accuracy14. Furthermore when using the first order upwind scheme a large amount of diffusion occur which can be avoided when using the second order upwind scheme15. This scheme uses Taylor expansion of the CV centred scalar values to calculate the corresponding face values. This can be done using the the following equation,
φf =φ+ ∇φ·~r, (2.9)
where~r is the vector from the CV center to the face center. It is not always an advantage to calculate with a higher order of accuracy scheme since this can disturb the conver- gence due to local flow fluctuations. For some meshes the solution can be of higher accuracy if the QUICK scheme is used, which is based on the second order upwind scheme. The QUICK scheme is commonly implemented when a hexahedral mesh is used, however it can be implemented also when other meshes are used. When the QUICK scheme is used with hybrid meshes the second order upwind scheme will be used for the volumes consisting of other mesh than hexahedral16.
To solve the gradient, ∇φ, see equation 2.9, it is possible to choose between differ- ent approaches in ANSYS Fluent. The three approaches available are; Green-Gauss Cell-Based, Green-Gauss Node-Based and Least Squares Cell-Based. When the Green- Gauss-based method is used the gradient is calculated in cell center c0 by the following equation17,
(∇φ)c0= 1 ν
∑
f
φ¯fA~f, (2.10)
13Ansys Inc, General Scalar Transport Equation: Discretization and Solution, Ansys Fluent Documentation, 2009
14Ansys Inc, General Scalar Transport Equation: Spatial Discretization, Ansys Fluent Documentation, 2009
15Nils R. B. Olsen, Hilde M. Kjellesvig, Three-dimensional numerical flow modelling for estimation of spillway capacity, Journal of Hydraulic Research, VOL 36, 1998, p.782
16Ansys Inc, General Scalar Transport Equation: Spatial Discretization, Ansys Fluent Documentation, 2009
17Ansys Inc, Evaluation of Gradients and Derivatives, Ansys Fluent documentation, 2009
where
φ¯f = φc0+φc1
2 (2.11)
and ν [m2/s] is the kinematic viscosity. In equation 2.11 the face value ¯φf is calculated from averaging the two neighborhood cell center values, φc0and φc1. In ANSYS Fluent the option of an implicit- and explicit scheme is available for the time discretization and the following partial differential term is discretized,
∂φ
∂t =F(φ). (2.12)
Also the chioce of 1st- and 2nd-order implicit time scheme is available. For both the implicit- and explicit scheme the equations are solved at each time step n. When using the explicit formulation the partial differential term F(φ)is evaluated by the following equation,
φn+1−φn
∆t =F(φn), (2.13)
leading to the explicit formulation of φn+1,
φn+1=φn+∆tF(φn), (2.14) where φn+1 is calculated using previous values φn. However when the implicit for- mulation is used the partial differential term F(φ) can be defined by the following equation,
φn+1−φn
∆t =F(φn+1). (2.15)
From equation 2.15 the implicit formulation of φn+1can be derived,
φn+1=φn+∆tF(φn+1). (2.16) One advantage by using the implicit scheme is that it is unconditionally stable with respect to the size of the time step while the explicit scheme can be advantageous for flows with transient moving waves, i.e shock waves, since it is then less computationally expensive and more accurate. However the explicit time scheme is more restricted than the implicit time scheme and can only be used in combination with the density implicit based solver. Also the explicit time scheme can not be used for incompressible flows18. Due to restrictions of model options in this project explained in following sections 2.2.2 and 2.2.6, the pressure based solver must be used, hence the implicit time discretization scheme must be used.
18ANSYS Inc, Temporal Discretization, Ansys Fluent documentation, 2009
2.2.2 Pressure or density based numerical solver
In ANSYS Fluent the choice between pressure and density based solver is available.
Decades ago the pressure based solver was introduced to solve low speed incompress- ible flows while the density based solver was used to solve high speed compressible flows. However over the years both models have been evolved and can now be used to solve a larger variety of flows than what they initially were intended to solve. The choice of solver for a specific fluid flow is not as clear nowadays as it was decades ago. Both solvers use the momentum equation to derive the velocity field and both solvers use the same discretization model mentioned in section 2.2.1. The most distinct separation between the solvers is that the density based solver uses the continuity equation to cal- culate the density field and the equation of state to calculate the pressure field while the pressure based solver uses a pressure correction equation extracted from the continuity and momentum equations to obtain the pressure field19. Due to certain restrictions on the multiphase model used in this project, which will be explained in section 2.2.6, it was only possible to use the pressure based solver.
The pressure based solver uses an algorithm where the constraint of the velocity field is obtained by solving a pressure correction equation. The pressure correction equation is derived under the condition that the velocity field, with pressure correction, satisfies the continuity equations. Since the set of equations are coupled and non-linear the algebraic equations must be iterated and solved several times until convergence is reached. In AN- SYS Fluent the option of two different pressure based solver algorithms is available, the segregated algorithm and the coupled algorithm. When using the segregated algorithm the governing equations are solved separately and since the equations are non-linear and coupled the equations are solved by several iterations20, see figure 2.1.
19Ansys Inc, Overview of flow solvers, Ansys Fluent Documentation, 2009
20Ansys Inc,pressure based solver, Ansys Fluent documentation, 2009
Figure 2.1: The computational steps when using the segregated algorithm.
When implementing the coupled algorithm a system of coupled equations of the mo- mentum equations and pressure continuity equation are solved simultaneously. One advantage by using the coupled algorithm is that the rate of convergence becomes much faster but on the other hand one disadvantage is that the memory requirement can get twice as large21. Because of the simplicity of the segregated algorithm it was chosen for the simulations in this project. In ANSYS Fluent there are four different options for pressure velocity coupling available in combination with the pressure segregated algorithm named SIMPLE, SIMPLEC, PISO and Fractional Step. In this project the SIMPLE velocity pressure coupling has been implemented since it is the simplest and least computationally expensive22.
2.2.3 Reynolds-Averaged Navier-Stokes Model
The most accurate method to calculate a fluid flow is to directly solve the Navier-stokes equations numerically without any approximation or averaging other than using one of the discretization methods mentioned in section 2.2.1. This method is called Direct Numerical Simulation, DNS. When implementing this approach the errors of the final solution can both be estimated and restrained. When using the DNS method the solution contains a substantial amount of information about the flow with a high degree of
21Ansys Inc, pressure based solver, Ansys Fluent documentation, 2009
22Ansys Inc, Pressure-Velocity Coupling, Ansys Fluent documentation, 2009
accuracy. However from an engineering point of view the information retained from the DNS method is excessive. Additionally because of the method’s complexity it is thoroughly computationally expensive and the time it takes to obtain a solution can be far too long. For these reasons this method can not be deployed too often, hence it is generally not suitable to use as a design tool and in consequence engineers tend to seek for a less complex method23.
One method that is eminently less complex, and that has been used in this project, is the Reynolds-Averaged Navier-Stokes, RANS, model, which is based on ideas intro- duced by Osborne Reynolds over a century ago. Engineers often use this model when computing specific attributes of a flow, for example the average force distribution. By using the RANS model the Navier-Stokes equations are time-averaged and fluctuations are regarded as a part of the turbulence. Due to the complexity of turbulence a single RANS model can not represent all turbulent flows, thus supplementary turbulence models must be used in combination with the RANS model24. This will be discussed further in section 2.2.4.
For an unsteady flow every variable can be written as,
φ(xi, t) =φ¯(xi) +φ0(xi, t), (2.17) where
φ¯(xi) = lim
N→∞
1 N
∑
N i=1φ(xi, t). (2.18)
In equations 2.17 and 2.18, φ(xi, t)is the variable that is going to be calculated, ¯φ(xi)is the ensemble averaged function of the same variable, φ0(xi, t)is the fluctuation and xiis the position vector. The time average of fluctuations are defined to be zero, ¯φ0(xi, t)= 0.
Using equation 2.17 the velocity, ui,- and pressure field, p can be described by,
p(xi, t) = ¯p(xi) +p0(xi, t), (2.19) and
ui(xi, t) =u¯i(xi) +u0i(xi, t), (2.20) hence both pressure and velocity are described by an average component and a fluctua- tion component25. By using the equations for the average velocity 2.20 and pressure 2.19 as well as the continuity equation 2.2 and the incompressible Navier-Stokes equation 2.5 the incompressible RANS equations in Cartesian tensor form can be derived,
∂ρ ¯ui
∂xi =0 (2.21)
23Ferziger, p.268
24Ferziger, p.292
25Ferziger, p.293
∂(ρ ¯ui)
∂t
+∂(ρ ¯uiu¯j)
∂xj
= −∂ ¯p
∂xi + ∂
∂xj h
µ
∂ ¯ui
∂xj + ∂ ¯uj
∂xi
i
+∂(−ρu0iu0j)
∂xj
+ρgi, (2.22) where u0iu0j is called Reynold’s stress term which needs to be modeled in order to find a solution to the equations26.
2.2.4 Turbulence models
Since the three-dimensional, time-dependent Navier-Stokes equations are not suitable to solve analytically, the Reynolds-averaged equations are used. However, simply using this method looses information, so additional estimations are needed27. There are several models available to approximate the flows, such as the CS-model, that is used for incompressible flows in two dimensions. Two widely used models for flows in three dimensions, that can be chosen in ANSYS Fluent are the k-e model and k-ω model.
Another model for three-dimensional flows is the SST model that uses features from both the k-e model and the k-ω model. A fifth turbulence model, that can also be used in ANSYS Fluent, is the SA model, that uses a semi-empirical transport equation28.
k-e model
The k-e model is a turbulence model that is based on eddy viscosity concept. The k-e model uses the Boussinesq hypothesis, which is used to relate Reynolds stresses to mean velocity gradients. It can be presented as
−ρu0iu0j =µt
∂ui
∂xj + ∂uj
∂xi
−2
3ρδijk, (2.23)
and
k= 1
2u0iu0i (2.24)
where µt[Pas] is the computation of the turbulent viscosity and k [J/kg] is transport of turbulence energy29. In this eddy viscosity concept based model, the eddy viscosity is written as
em= cµk
2
e , (2.25)
where em[Pasm3/kg] is the eddy viscosity, cµ[s−1] (default as 0.09) is a constant and e [J/kgs] is rate of dissipation30. In addition to equation 2.23 there are two other equations for the k-e model, the equations for k and e, respectively. There are different equations to
26Ansys Fluent, 2009
27Tunger Cebeci, Turbulence Models and Their Application, Germany: Horizons Publishing Inc, 2004, p.1
28Cebeci p.2-13
29Ferziger p.294
30Cebeci p.2
describe k and e, which depend on the choice of k-e model. The transport of turbulence energy in accordance with the k-e standard model is defined as
Dk Dt = ∂
∂xk h
ν+em σk
∂k
∂xk i
+em
∂ui
∂xj + ∂ui
∂xi
∂ui
∂xj −e, (2.26) where t [s] is the time, xk, xjand xiare coordinates, ν [m2/s] is the kinematic viscosity of the fluid, σk [m−1] (default as 1.0) is a constant and ui[m/s] is the velocity field of the flow31. The rate of dissipation in the k-e model is defined as
De Dt = ∂
∂xk h
ν+em σe
∂e
∂xk i
+ce1e kem
∂ui
∂xj + ∂ui
∂xi
∂ui
∂xj −ce2e2
k , (2.27) where σe[m−1] (default as 1.3), ce1[m−1] (default as 1.44) and ce2[m−1] (default as 1.92) are constants32.
In ANSYS Fluent there are three variants of the k-e model; the standard model, the renormalization group (RNG) model and the realizable model. The RNG- and realizable model are newer and modified versions of the standard k-e model and have shown significant enhancement over the standard k-e model33. The equations for k and e are slightly different from the standard model. The disadvantage of the RNG- and realizable model are that they are more computationally expensive than the standard model.
2.2.5 Boundary conditions
To be able to calculate the fluid flow it is important to specify the boundary conditions.
The No-Slip boundary condition has been used for the walls inside the discharge tunnel and assumes that the velocity, both tangential and normal to the walls, of the fluid particles in contact with the wall is zero. This phenomenon occurs due to the viscosity of the fluid. In the main parts of the spillway tunnel in this study, the water will not fill up the whole volume, hence it is important to consider the free water surface that works as a barrier between the water and air inside the spillway tunnel.
Fluid flows with free surfaces are specifically complicated to describe mathematically since they give rise to moving boundaries. In order to be able to treat these kinds of flows it is crucial to keep track of the free surfaces. However, the fluid surface is only known at the initial time and therefore it must be calculated as a part of the solution at each time step. Most commonly the free surface is a water-air mixture, which is also the case in this thesis34. If phase change can be neglected the following kinematic and dynamic boundary conditions can be applied.
31Cebeci p.5
32Cebeci p.5
33SAS IP, Inc, ANSYS Fluent, Sharcnet, 2015
34Ferziger, p.381
The kinematic condition implies that the free surface is a sharp boundary separating water and air and does not allow any flow through it,
(v−vb) ·n]f s=0 (2.28) where f s stands for the free surface, v [m/s] is the velocity of the fluid at the surface, vb [m/s] is the velocity of the free surface and n is the normal to the surface. The dynamic condition implies that the forces acting on the fluid particles at the surface must be in equilibrium, which indicate that the normal forces must be equally large and in opposite direction on each side of the boundary. However the tangential forces must be equally large on each side but in the same direction, see equations 2.29, 2.30 and 2.31.
(n·T)l·n+σK= −(n·T)g·n, (2.29)
(n·T)l·t− ∂σ
∂t = (n·T)g·t, (2.30)
(n·T)l·s− ∂σ
∂s = (n·T)g·s. (2.31)
Here T [K] is the temperature, σ is the surface tension [N/m], K [m−1] is the curvature of the surface,(n, t, s)are unit vectors in a orthogonal coordinate system where n is the normal to the free surface while t and s are the tangents to the surface. Additionally l and g stand for liquid and gas respectively.
The shape and location of the free surface is only known at the initial time and must therefore be calculated iteratively, which increases the computational cost and complex- ity of the calculations greatly. Hence it is not as simple as to implementing the earlier stated boundary conditions directly, see equations 2.28 - 2.31. There exists different methods to calculate the shape of the free surface and they can be categorized into two main groups; the ones that define the surface to be a sharp boundary and the ones that does not define it to be a sharp boundary35. The method used in this thesis belongs to the second group, hence does not define the free surface to be a sharp boundary. This will be explained more thoroughly in section 2.2.6.
2.2.6 Multiphase flow modeling
In nature it is common that a flow consists of several different phases, i.e gases, liquids and solids. When numerically calculating a flow not only the physical phases mentioned above are considered as a phase, also particles of different sizes in the same material can be defined as different phases if their interaction with the flow is different. When a flow consists of two or more phases it is called a multiphase flow36. Different multiphase flows can be simulated in ANSYS Fluent and in this project a gas-liquid flow has been
35Ferziger, p.383
36Ansys Inc, Multiphase Flow Regimes, Ansys Fluent Documentation, 2009
simulated which consists of a moving liquid flow with a free surface surrounded by gas. Currently there exist two different types of multiphase flow approaches, the Euler- Lagrangian approach and the Euler-Euler approach. In this project the Euler-Euler approach has been used. When using this approach it is possible to choose between three different models; The Volume of Fluid, VOF, model, the Mixture model and the Eulerian model37. The most complex model out of these three is the Eulerian model and both the Mixture model and Euler model are suitable to use for flows with phase transition. On the other hand the VOF model is a surface tracking model and is therefore suitable to use when calculating fluid flows with free surfaces. Since the flow in this project have a free surface the best choice of multiphase regime is the VOF model38.
Volume of Fluid model
The VOF model can be used for fluid flows containing two or more phases. Two important restrictions of the VOF model are that the phases can not penetrate each other and that the pressure based solver must be used39. When using the VOF scheme the shape of the free surface is derived by dividing the geometry into a finite number of cells and then calculating the volume fraction of each cell that is partially filled, i.e that has both water and air inside. In this way the free surface can be distinguished from the rest of the fluid flow. For each control volume, CV, the following is true,
∑
q 0αq=1, (2.32)
where αqis the filled fraction of the qthfluid. Only the following three cases are allowed, αq= 1, αq= 0 and 0<αq<1. When αqis equal to 1 it means that the control volume only is filled with the qthfluid and when αqis equal to zero it means that the control volume does not contain that fluid. Lastly when the filled fraction has a value larger than zero but smaller than one it implies that the CV contains a mixture of that fluid and other fluid(s). For this project only two phases are considered, the primary phase air and the secondary phase water. It is common to select the primary phase to be the fluid with lowest density. For each CV, an equation for the filled fraction, αq, is solved. It is common practice to define the free water surface to have a volume fraction of 0.5 meaning that the CVs constructing the free surface contains 50 % water and 50 % air.
Also on the free surface the kinematic and dynamic boundary conditions apply40. In this project the filled fraction of the water phase can be calculated by solving following continuity equation,
1 ρw
∂ρwαw
∂t + ∇ · (αwρwvw) =Sαw+m˙aw−m˙wa], (2.33)
37Ansys Inc, Approaches to Multiphase Modeling, Ansys Fluent Documentation, 2009
38Ansys Inc, Approaches to Multiphase Modeling, Ansys Fluent Documentation, 2009
39Ansys Inc, Overview and Limitations of the VOF Model, Ansys Fluent Documentation, 2009
40Ferziger, p.384
where ρw [kg/m3] is the density of the water, αw is the volume fraction of water, vw [m/s] is the velocity of the water fluid, Sαw [kg/s] is the source term and is by default set to zero, ˙maw [kg/s] is the mass flow from phase air to phase water and ˙mwa [kg/s] is the mass flow from phase water to phase air. It is possible to solve this equation with both an implicit and explicit scheme. The volume fraction of the air can then be calculated from the result from this equation, 2.33, combined with equation 2.32. A limitation of the VOF approach is that it is not possible to use the second order implicit time discretization scheme in combination with the explicit VOF scheme41.
Wall function
When working with fluids of high Reynolds number near walls, the existence of the wall needs to be taken into consideration to be able to use boundary conditions at the very thin viscous sublayer since it is troublesome to operate a sufficient amount of gridpoints in this layer to solve it42. These adjustments are called the wall function43. The wall function is a logarithmic law, relating a particle’s velocity to the logarithm of the distance between the particle and the wall. As an option to the standard wall function, ANSYS Fluent offers a function called non-equilibrium wall function, that is more complicated since it, in contrast to the standard wall function, involves two equations instead of one.
The non-equilibrium wall function has shown to be advantageous over the standard wall function for spillway tunnels44.
2.2.7 Mesh
Choice of mesh
As prevously mentioned, the Navier-Stokes equations can usually not be solved analyti- cally for fluid flows. Therefore a continuity of subdomains of the geometry, i.e a mesh, has to be used to analyze the fluid flows. When choosing a suitable mesh, the accuracy of the solution, as well as the convergence rate (computing effort), has to be taken into account. In general a finer grid gives a finer solution, but is also more computationally expensive45. There are several types of meshes, for example triangular and quadrilateral meshes in 2D as well as tetrahedral, hexahedral and polyhedral meshes in 3D. Different meshes are preferable for different calculations. It is also possible to use different meshes in different parts of the geometry, for example a finer grid size in sections where it is needed, like where the flow is more turbulent, without making the computation more expensive than necessary, if the grid size is small enough for the areas where the flow is not very turbulent.
41Ansys Inc,Volume Fraction Equation, Ansys Fluent Documentation, 2009
42Cebeci p.5
43Ferziger p.298
44Dargahi p.903
45Ferziger, p.344
For simulations in 3D, tetrahedral meshes are commonly used, however they can cause unanticipated issues when oscillations can make convergence unattainable46. Another difficulty with tetrahedral meshes are that the angle between the cell face and the line that connects the centers of two neighbouring cells are not allowed to have a value far from 90 degrees, since that can cause extensive flaws and convergence issues47. Hexahe- dral and polyhedral meshes have more nodes in each cell as well as more neighbouring cells, which allows improved gradient approximation, for example studies have shown that a polyhedral mesh can provide the same accuracy as a tetrahedral mesh when using one fourth of the number of cells as the tetrahedral mesh, and as little as a tenth of the computational time48. The choice between using a polyhedral mesh or a hexahedral mesh can depend of the type of flow, where the hexahedral mesh is beneficial for a flow whose streamlines have similar direction as the grid lines, since it causes less smearing than the polyhedral mesh. However, if the flow pattern is more turbulent, a polyhedral mesh is more suitable49. Therefore, if the flow varies throughout the model, sections of hexahedral meshes, together with sections of polyhedral meshes could be a good solution.
Even though tetrahedral meshes are less desirable than hexahedral and polyhedral meshes, it is sometimes not possible to implement a hexahedral and polyhedral mesh due to the geometry of the model that is being meshed. Then a tetrahedral mesh could be a sufficient choice of mesh.
Also using an adaptive grid can be beneficial when simulating flows with a moving surface since it gives flexibility, and can help avoiding false diffusion, that smear out gradients and degenerate the results50. Having a finer grid size near the walls is bene- ficial to avoid incorrect pressure distribution, and at smaller flow depths is crucial for a spillway tunnel to retrieve reliable results51. There are several softwares that create meshes. One is GAMBIT by ANSYS, where two- and threee-dimensional meshes can be created and imported to ANSYS Fluent.
Grid convergence index and quality control
Since the results of simulations are obtained with discretizations, errors can occur, which depend on the mesh. Therefore it should be examined if the results changes distinctly when the size of the mesh is changed, to analyze if the mesh size should be decreased,
46Ferzinger, p.343
47Ferzinger, p.342
48Milovan Peric, Steven Fergusson The advantage of polyhedral meshes, Dynamics 24:45, 2012
49Gus Nasif Ronald Barron Ram Balachandar Jet Impingement Heat Transfer: Stationary Disc, International Journal of Surface Engineering Materials Technology. 4. p34-38, 2014
50Olsen p.782
51Dargahi p.905
and to examine whether the mesh gives results that converge. To examine if the results generated with the mesh converges, an approximation of error bands, known as grid convergence index (GCI) can be checked. This can be done by using the formulas
e= f2− f1
f1 , (2.34)
where f is the concerned quantity, such as water depth that is obtained with the fine, middle and coarse grid size, and
p= ln(ff3−f2
2−f1)
ln(r) , (2.35)
where r is the grid refinement ratio (how much the grid size is decreased between the different meshes, in this project the grid refinement ratio 2 was used, meaning the grid size was halved when the mesh was refined), to find the relative error and order of convergence, respectively. After obtaining these values, they can be inserted into
GCIf ine= |e|Fs
rp−1, (2.36)
where Fsis a safety factor that should span between 1.25 and 3, where 3 is reccomended52 and was used in this study, to find the GCI for the finest mesh used in the simulations.
The grid convergence index is based on Richardsson Extrapolation and by using a second- order method it relates the results from an arbitrary grid refinement test approximately to the results that are expected from increasing the grid size by two53. A low grid convergence index implies that the solution does not depend considerably on the grid size. How low the GCI should be to be "good" depends on the case and depends on what the safety factor is set as. To be able to reach convergence and get an accurate solution, the quality of the mesh also has to be sufficient. There are several indicators of the quality of the mesh. The three most important indexes of mesh quality are skewness, face squish index and orthogonal quality, who all vary between 0 and 1. For the skewness and face squish index, 0 is the optimal value and 1 is the opposite of the optimal value. For the orthogonal quality, 1 is the most favorable value and 0 is the least favorable value.
2.3 Discharge tunnel
Since the spillway tunnel is one of the structures in hydropower schemes that is used to affirm security during floods, the design of the spillway tunnel is of high importance54.
52P. J. Roache Quantification of uncertainty in computational fluid dynamics, Annual Reviews Inc, 29:123–60, 1997, p.136
53P. J. Roache Perspective: a method for uniform reporting of grid refinement studies, Journal of Fluids Engineering, 116(3):405–413, 1994
54MU Zhenwei, Zhang Zhiyan, ZHAO Tao, Numerical simulation of 3-D Flow Field of Spillway based on VOF Method 2012 International Conference on Modern Hydraulic Engineering, 2012 p.808
Tunnels are a good way to divert water since with them it can be avoided to disturb the natural landscape and because of their structural stability55. The material of the spillway tunnel is usually rock. Spillway tunnels made from soft ground can also be used, but are refrained when possible56. On the inner wall of the spillway tunnel, the tunnel can be lined with concrete, roller compacting concrete or steel to protect and enhance the strength of the tunnel. If the spillway tunnel is made from good sound rock, usually no lining is needed57. The necessity of lining is also dependent on the water pressure in the spillway tunnel and the desire of low hydraulic resistance58. During the constructing of a dam, a tunnel is often used to divert water flowing towards the dam away from the construction site, and this tunnel can after the completion of construction be converted into a spillway tunnel, to cut expenses of building two tunnels59.
The shape of the cross section of the spillway tunnel varies. A structurally beneficial shape is a circular cross section, since it has high stability during high pressures. Other commonly used shapes that are non-circular have a flat bottom, arched ceiling and walls that are gently flaring or almost vertical. The most commonly used cross section is the
"horse-shoe" shape, that has a flat bottom and arched ceiling and walls, in the shape of a horseshoe60. The cross section can also vary throughout the spillway tunnel, as in the spillway tunnel used in this study. Depending on the alignment, the spillway tunnel can be classified differently. A spillway tunnel with a small slope is simply called tunnel.
If the spillway tunnel has vertical alignment it is called a shaft, and if it has a steep inclination it is called an inclined shaft61but is still a spillway tunnel. When designing a spillway tunnel, some of the aspects that need to be considered to optimize the spillway tunnel is to make it as short as possible, align it with the intake and surge shaft, as well as away from critical fault planes in the structure of the surrounding mountain62. Before this study a "flood discharge sluicing sediment tunnel was transformed from the diversion tunnel into a “dragon-up” type pressurized short tunnel whose inlet is non-pressure flow, including imported open channel, inlet tower, inclined shaft section, diversion tunnel with slow slope section, and exit tunnel section"63, which was used in this study.
55M M Dandekar, K N Sharma, Water Power engineering, Dehli: Vikas Publishing House PVT LTD, 1979, p.287
56Dandekar p.287
57Dandekar p.289
58Dandekar p.290
59John A. Roberson, John J. Cassidy, M. Hanif Chaudhry Hydraulic Engineering, 2nd ed, United States of America: John Wiley Sons, Inc, 1998, p.275
60Dandekar, p.288
61Dandekar p.288
62Dandekar p.296
63James Yang, Swedish student diploma work in China, Vattenfall, 2019
