Aeration and risk mitigation for flood discharge tunnel in Zipingpu water conservancy project

(1)

IN

DEGREE PROJECT THE BUILT ENVIRONMENT, SECOND CYCLE, 30 CREDITS

STOCKHOLM SWEDEN 2020,

Aeration and Risk Mitigation for Flood Discharge Tunnel in

Zipingpu Water Conservancy Project

JORGE CONTRERAS MORENO

KIBRET DAWIT GHEBREIGZIABHER

KTH ROYAL INSTITUTE OF TECHNOLOGY

SCHOOL OF ARCHITECTURE AND THE BUILT ENVIRONMENT

(2)

(3)

Aeration and Risk Mitigation for Flood Discharge Tunnel in

Zipingpu Water Conservancy Project

J ORGE C ONTRERAS M ORENO

K IBRET D AWIT G HEBREIGZIABHER

Master of Science Thesis

Stockholm, Sweden 2020

(4)

TRITA-ABE-MBT- 20191 ISBN: 978-91-7873-574-7

KTH School of ABE SE-100 44 Stockholm SWEDEN

Department of Civil and Architectural Engineering

(5)

i

Abstract

The importance of hydraulic structures has become an essential mitigating mean for floods that occur more often due to climate change. Thus, the importance and safety of flood discharge tunnels has promoted further studies and experiments on the topic to mitigate damages, such as cavitation that arise because of high speed flows.

After an experimental study on a physical model was carried out on the flood discharge tunnel in Zipingpu Water Conservancy project, a CFD model was designed and simulated in the commercial software ANSYS Fluent. The simulations aimed to evaluate and examine the risk for cavitation in the tunnel, examine the design problems of the structure and analyse the installed aerators for the mitigation of cavitation. Moreover, using CFD models as a complementary form to physical models was analyzed.

A three dimensional geometry of the discharge tunnel was built in ANSYS Spaceclaim and the mesh conducted with ANSYS mesh generator. The known boundary condition such as the design flow conditions, velocity inlet, pressure inlets and pressure outlet were set. For the model a multiphase VOF scheme with RANS approach, k-ϵ turbulence model and a standard wall function was set.

The results from the initial simulations showed that the discharge tunnel was under cavitation risk, since the recorded cavitation index in the tunnel was below 1.8. After having revised the layout of the aerators in order to mitigate cavitation risk, the results from the simulations with added aerators were sufficient to mitigate the risk as the cavitation index was still below 1.8.

The results for the cavitation index remained unchanged even in the simulated models with a different solver setup that were used in the comparison with the experimental data in order to validate them.

As a conclusion, it was recommended that the tunnel design has to be revised and improved by adding more aerators and air vents to mitigate the cavitation risk. Furthermore, more studies on the discharge tunnel or similar tunnels with similar conditions should be carried out in order to validate the results of this study and determine if numerical models are preferable to physical models.

Keywords: Flood discharge tunnel, cavitation risk, cavitation index, aerators, CFD model.

(6)

(7)

iii

Sammanfattning

Betydelsen av hydrauliska strukturer har blivit ett väsentligt förebyggande medel för översvämningar som förekommer oftare på grund av klimatförändringar. Således har vikten och säkerheten för översvämning tunnlar främjat ytterligare studier och experiment om ämnet för att förebygga skador, såsom kavitation som uppstår på grund av hög hastighets flöden.

Efter att en experimentell studie av en fysisk modell genomfördes på avrinningstunneln i Zipingpu Water Conservancy projekt, genomfördes en CFD-modell i den kommersiella programvaran ANSYS Fluent. Simuleringarna syftade till att utvärdera och undersöka risken för kavitation i tunneln, undersöka strukturens konstruktion problem och analysera de installerade luftnngsmekanismer för att minska kavitation. Dessutom analyserades användning av CFD-modeller som komplement till fysiska modeller.

En tredimensionell geometri för avrinningsstunneln byggdes i ANSYS Spaceclaim och nätet genomfördes med ANSYS nätgenerator. Det kända gränstillståndet såsom designflödesbetingelserna, hastighetsinloppet, tryckinloppen och tryckutloppet inställdes. För modellen sattes ett flerfasigt VOF-schema med RANS-tillvägagångssätt, k-ϵ turbulensmodell och en standard vägg funktion.

Resultaten från de initiala simuleringarna visade att urladdningstunneln var under kavitationsrisk, eftersom det registrerade kavitationsindexet i tunneln var under 1,8. Därefter redigerades beläggningen av luftningsmekanismerna för att minska risken för kavitation och resultaten från simuleringarna med tillagda luftningsmekanismer var inte tillräckliga för att förebygga kavitationsrisken eftersom kavitationsindexet fortfarande låg under 1,8. Resultaten för kavitationsindexet förblev oförändrade även i de simulerade modellerna med en annan lösningsuppsättning som användes i jämförelsen med experimentella data för validerings skull.

Som en slutsats rekommenderades att tunnel designen måste revideras och förbättras genom att lägga till flera luftningsmekanismer och luftventiler för att minska kavitationsrisken. Vidare bör fler studier på urladdningstunneln eller liknande tunnlar med liknande förhållanden genomföras för att validera resultaten av denna studie och bestämma om numeriska modeller är att föredra framför fysiska modeller.

Nyckelord: Avrinningstunnel, kavitationsrisk, kavitationsindex, luftningsmekanism, CFD- modell.

(8)

(9)

v

Preface

The diploma project reported in this Master Degree thesis was carried out during the first half of 2020. Due to the Covid-19, we were unable to travel to China as originally planned. Instead, the work was performed at Royal Institute of Technology (KTH) and still in close cooperation with Tsinghua University, Beijing. In the light of the situation, necessary adjustments and arrangements were made in terms of project topic, layout, supervision, means of communications etc.

We would like to thank Prof. Yongliang Zhang, Tsinghua University, and Prof. James Yang, Vattenfall and KTH, and for the supervision, access to data, advice and discussions. We would also like to devote our thanks to Ph.D. student Shicheng Li and Dr. Penghua Teng at Department of Civil and Architectural Engineering, KTH, for all the help with numerical simulations including program learning during the performance of the project.

We are grateful to our examiner Anders Ansell at the Division of Concrete Structures, Royal Institute of Technology for the coordination. Jorge Contreras Moreno would like to give special thanks and gratitude to the National Council of Science and Technology of Mexico (CONACyT) for providing financial support for living and studying at the Royal Institute of Technology in Stockholm.

The diploma project is funded by Energiforsk AB within the frame of dam safety (www.energiforsk.se), with James Yang as coordinator. The project has been going since 2004, with more than 115 university students who have done their diploma work in different universities of technology in China.

Stockholm, June 2020 Jorge Contreras Moreno Kibret Dawit Ghebreigziabher

(10)

(11)

vii

List of acronyms

CFD Computational Fluid Dynamics CFRD Concrete Faced Rockfill Dam

CV Control Volume

GCI Grid Convergence Index FVM Finite Volume Method

RANS Reynolds Averaged Navier Stokes VOF Volume of Fluid

(16)

(17)

1

1 Introduction

With the occurrence of changing climatic and temporal characteristics taking place on earth due to environmental degradation and pollution, heavy rains and flooding are becoming more common. China is a country hit by the effects of climate change, which apart from flooding have also created water scarcity and drought problems ((Ministry of water resources, 2006)).

To tackle these problems, the country has built numerous water conservancy projects that serve as flood control, water supply and hydropower generation.

In this study the discharge tunnel at the Zipingpu Water Conservancy project, which is used for flood control, has been analysed. A numerical model will be established using a computational fluid dynamics (CFD) software called ANSYS. The model will then be validated in respect to the results procured from the physical scale model that has been constructed at Tsinghua University in Beijing, China. Numerical models are preferable because they are more cost and time effective in comparison to physical models. But, nevertheless, physical scale models are used to validate numerical models.

Aims and objectives

To complete this study various aims and objectives were set and are presented here.

The aim of the study is to:

- Analyse if there are design problems in regard to cavitation and optimize the design of the tunnel to minimize the risks of cavitation

- Analyse whether the shape and the installed aerators in the discharge tunnel are suitable to mitigate cavitation risks

- Analyse whether a numerical model is preferable to a physical scale model.

To achieve this aims the following objectives had been set:

- Establish a numerical model in ANSYS Fluent for 3D complex fluid flow, to simulate the flow characteristics such as aeration and cavitation risk mechanisms

- Calibrate the numerical model through a comparison with the results acquired from the experimental study on the physical scale model.

(18)

Limitations

This study was to be conducted at Tsinghua University in Beijing where a physical scale model was conducted, and other studies of the Zipingpu Water Conservancy project were conducted.

However due to complications that arose in connection with the Covid-19 pandemic, the study was conducted in Stockholm instead. This meant that limitations had to be considered in the study and especially in the results. The encountered limitations were as follows:

- The ANSYS R3 2019 version that was used had an educational license which means that the mesh element number was limited to a maximum of 512 000 mesh elements and therefore simulations with finer mesh elements were not able to run.

- The experimental data that was used in the post-processing section and validation to compare the calculated results was not the data for the tunnel of this study, but instead of a similar tunnel that is a hydraulic structure in the Zipingpu Water Conservancy project. The experimental data was extracted from an older study that was conducted on the other tunnel since it was said to be similar after deliberation with the supervisors (Hamberg and Dahlin, 2019). Therefore, the comparisons that were made between the calculated results and experimental data might deviate slightly or could have been better if the right experimental data for the tunnel in hand was available.

- Time constraint was a limitation. The study took off behind schedule but was concluded within the deadline that was set beforehand. This affected the process of the grid independence study and possibility of running different simulations with different types of interpolation and discretization methods.

(19)

3

2 Background

Water scarcity and water conservancy projects

China is home to the largest population on earth. With its ever-growing economic development, population and living standards, must meet the needs of its citizens. One of the biggest challenges that the country is facing is water scarcity and meeting the consumption needs (Jiang, 2009). This crisis has given the Chinese government the incentive to establish a water conservancy plan (Liu and Yang, 2012). However, this plan although aimed to achieve water sustainability may also cause environmental and socio-economic repercussions if correct assessments are not carried out beforehand.

China’s water scarcity can be connected to different causes: natural characteristics, where the southern part of China holds most of the water resources, whereas the northern part China that accounts for 45.2% of the population has only 19.1% of the water resources; economic development, where the hasty industrialization resulted in a 15% consumption of world natural resources for water; population growth, where China is home to over 13 billion people (20% of the world population) but only has 6.5% of the world’s total freshwater resources; water resource management, where the water conservancy plans are carried out poorly; water quality, where poor water quality in the rivers and basins add up to the scarcity and threaten China’s economic development, food security and life quality (Jiang, 2009). Moreover, due to the climatic and temporal characteristics that result in an annual precipitation of 60-70% in summer season add to the water scarcity problem (Liu et al., 2013).

China has been tackling water scarcity by raising funds for investments to water conservancy projects and to achieve a sustainable water management. China has more than 87000 dams, the established plan and investments are aimed to repair and sustain more than half of them and moreover, as mentioned, to manage them sustainably. However, the introduced water management plans set up by the government focus mainly on the quantity and not quality of the water. The quality of water is to be taken in account, considering that almost 40% of the rivers are polluted and 80% of the lakes suffer from eutrophication (Liu and Yang, 2012).

As stated earlier, China being one of the fastest growing and developing economies have dedicated resources and investments to counteract the water scarcity present in the country.

Which lead to a boom of construction of hydraulic structures such as water conservancy projects, that serve different purposes such as hydropower generation, urban water supply, irrigation supply and flood control (Liu et al., 2013). Hence, China is home to the largest number of dams and the largest hydropower plant, i.e. the Three Gorges Hydropower Project (TGHP) (Jiazhu, 2002).

(20)

The construction of dams and reservoirs has had a positive impact in flood control measures (Liu et al., 2013). Moreover, China has a total reservoir capacity of 932.312 billion m³ (Ministry of water resources, 2006) and an installed hydropower generation of 352 million kW year 2018 that accounts for almost 20% of China’s electricity generation and 27% of the world’s capacity (International Hydropower Asssociation, 2019, Jia, 2016). Furthermore, the construction of dams and reservoirs has facilitated the possibility of farming areas that previously was not possible due to the temporal and climatic distribution of precipitation. This resulted in a 60.35 million ha of irrigated area by 2010, that is pathing the way to food security for the country (Liu et al., 2013).

China has executed numerous water conservancy projects to alleviate their water scarcity. To ensure quality, sustainability and financial support various policies and plans have been carried out by the government. Amongst the recent plans is the “Tenth Five-Year Plan” that consist of crucial water conservancy projects such as Baise, Linhuaigang, Shapotou, Ni’erji and Zipingpu (MINISTRY OF COMMERCE, 2002).

The Zipingpu water conservancy project, which this study is focused on, is located in Sichuan province that is deemed to be one of the most important bases for hydropower and water resource hence it has an ample water resource, and therefore, many hydraulic plants are built in this area. The Zipingpu water conservancy project started construction in 2001 and concluded in 2006. It is mainly designed for the purpose of irrigation and urban water supply, but it also functions as a hydropower plant with a capacity of 760 MW and flood control (Xinhua News, 2002). The Zipingpu water conservancy project consists of different hydraulic structures and is located in the upper reaches of the Minjiang river, 60 km northwest of the capital of Sichuan province, Chengdu, and 9 km west of Dujiangyan City (Tanchev, 2014). To name some of the main structures: a concrete faced rockfill dam (CFRD) with a height of 156m, a spillway, a sand blasting hole, a hydropower system and two flood discharge tunnels that act as the main structures for flood control and sand discharge in the occurrence of a flood . One of the flood discharge tunnels will be the focus of this study.

Flood discharge tunnel

Apart from the water scarcity that is present in China, flooding is also another problem that the country tackles with. Due to its climatic characteristics and geographical location, the country has been subject to the Eastern Asian monsoon that results in heavy flooding (Zhang and Liu, 2006). China's investment in water conservancy projects and the adherent policies also covers the flood control and management in the country, and therefore, China has had achievements in their flood control programs (Ministry of water resources, 2006).

Among flood control hydraulic structures are flood discharge tunnels. A flood discharge tunnel is commonly structured by a sloping section, a toe curve and an approximately sloping or almost horizontal section that joins into an energy dissipator or tailwater at its end (Khatsuria, 2005).

(21)

5

The performance and stability of discharge tunnels is dependent on the design schemes used to tackle different problematics regarding shape, size and flow properties. Regarding the shape and size, the toe curve is of importance for the performance of a tunnel. In theory, the flow in a tunnel is three dimensional and therefore, given the high velocity of the flow, the pressure on the tunnel surface would be highly positive as the flow fastens to the surface and cavitation damage occurs. Thus, the toe radius must be designed accordingly, too large of a radius is economically unsustainable and too small might cause unwanted and risky flow conditions, therefore, the toe radius should be larger than the tunnel radius before the toe start (i.e. 2.5 - 10.5 times the radius of a tunnel) (Khatsuria, 2005). Another important aspect is the cross- section shape of the tunnel. Given that it is a tunnel in question, a circular cross-section is the obvious choice for its stability and uniformity under high pressure. However, the most frequently used shape is in the form of a “horse shoe”, where the ceiling and walls of the tunnel are arched and the bottom is flat to resemble a horse shoe(U.S. Army Corps of Engineers, 1980, Dandekar and Sharma, 1979).

In order to obtain an economically sustainable dam construction, the diversion tunnel used for deviating the river flow away from the construction can be then used as a permanent flood discharge tunnel/tunnel spillway upon completion of the project by adding an energy dissipator to it (Tian et al., 2009). The energy dissipator at the outlet of a discharge tunnel is commonly a flip bucket or a stilling basin, and is critical to achieve a smooth transition from the tunnels cross-section shape into a flat bottom to ensure the stability and integrity of the structure (Khatsuria, 2005). Furthermore, also for economic reasons, in the design of the cross-section the diameter of the tunnel is held to a minimum without undermining the purpose to be fulfilled by it, hence a tunnel is never allowed to flow with full capacity in order to leave space for air flow (Khatsuria, 2005). Therefore, it is advised to design a spillway tunnel with a ratio of ¾ or

⅞ of the full flow capacity to accomplish and allow a balanced air-water flow to avoid unsafe flow conditions (Khatsuria, 2005). The impacts of these unsafe flow conditions can be seen due to the flow's high velocity and low pressure, which in return can result in cavitation damages, hence the need for aeration.

For the purpose of this study one of the discharge tunnels in the Zipingpu water conservancy project was considered. According to the design information of the Zipingpu Water Conservancy project, the flood discharge sluicing sediment tunnel was transformed from a 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.

This tunnel has a length of 720.55 m, the inlets bottom elevation and exit elevation are 800 m and 745.156 m, respectively. The given discharge rates are 1530.87 m³/s and 1666.74 m³/s at the design water level and the check level respectively, with a maximum flow velocity of 45 m/s and a flow rate of 212.9 m³/s.m at the outlet of the discharge tunnel.

(22)

General description of Cavitation

As mentioned above, cavitation is a severe damage that affects the stability and integrity not only of a discharge tunnel but other hydraulic structures such as spillways and chutes (Yazdandoost and Attari, 2004). In this chapter the cavitation phenomena will be covered and explained.

Cavitation is when voids are formed in a liquid medium and they can be classified in two groups: vaporous cavitation, when the void is filled with water vapor and gaseous cavitation when the void is filled with a gaseous medium (Falvey, 1990). Often, it is compared and examined with the aspect of boiling water at a local atmospheric pressure. The vapour pressure increases in boiling water. When the vapour pressure of the boiling water reaches and is equal to the local pressure, bubbles are formed due to the transformation phase of the water to vapour at the boiling point (Khatsuria, 2005, Falvey, 1990).

The pressure obtained in the water can be seen as a governing factor for the phenomena of cavitation. Therefore, boiling can be obtained at lower pressures consequently as the pressure decreases even up to a level where it can be obtained at room temperature. Nevertheless, bubbles are still formed even in this case and are defined as vaporous cavitation (Khatsuria, 2005, Falvey, 1990). Boiling and cavitation, although intertwined with each other, are not the same thing. Both illustrate a phase change from liquid to vapour, where boiling changes in temperature but holds the local pressure constant, and the opposite for cavitation where the temperature is held constant while the local pressure changes (Falvey, 1990).

To describe gaseous cavitation, the example of bubbles formed in a bottle of carbonated water is used. Before a bottle is opened, the carbonated water is still because it’s kept sealed under high pressure. When opened, the pressure in the bottle decreases, and bubbles form as the liquid becomes saturated and the carbon dioxide diffuses (Khatsuria, 2005, Falvey, 1990).

The damage caused by cavitation occurs when the bubbles filled with vapor pressure collapse.

When the bubbles collapse or explode along a solid structure, the high pressure due to the collapse generates a force capable of damaging hydraulic structures such as a discharge tunnel (Khatsuria, 2005). The collapse mechanism of a single bubble is described as a process of phases where the bubble diameter decreases until it reaches a minimum and then increases again or rebounds. This process is repeated several times, and in every cycle the bubble diameter decreases eventually down to a microscopic size, as shown in Figure 2.1. A shock wave with a velocity equivalent to the speed of sound in water is formed in the rebound phase. It is stated that the shock wave emits a pressure 200 times the ambient pressure at the collapse site (Khatsuria, 2005, Falvey, 1990).

(23)

7

Figure 2.1. Collapse mechanism of a single bubble near a solid surface (as shown in Khatsuria 2005 and Falvey 1990)

Cavitation damages

Cavitation damages are more prone to be found in the structures where there are surface irregularities, such as rough spots and joints between structures (Falvey, 1990, Chanson, 1989).

Basically, in the regions where the surface irregularities are present, a flow separation might be formed, and the pressure lowered. And in the presence of high flow velocities, it results in bubble formation, due to the low pressure (below the local pressure), which will later on collapse when it reaches a region of higher pressure and be liable to cavitation damages (Chanson, 1989). To be noted is that the cavitation bubble consists of several bubbles and is referred to as a cavitation cloud (Falvey, 1990).

The surface irregularities on spillway are directly connected to the surface roughness, of which they are classified as singular or isolated roughness and uniformly distributed roughness.

Singular roughness may also be specified as: offset into the flow, offset away from the flow, abrupt curvature or slope away from the flow, voids or grooves, roughened surface and protruding joint (Khatsuria, 2005, Falvey, 1990). As a result of the unusual change in the flow at the irregularities, cavitation occurs due to the turbulence in the shear region in all the listed cases of singular roughness. As cavitation is formed, cavitation damage will occur as the bubbles collapse either close to the flow boundary or in the flow itself depending on the roughness shape (Falvey, 1990).

In difference to singular roughness, the uniformly distributed roughness develops cavitation in the flow due to fluctuations that take place over a larger area. The fluctuations that arise, can be a result of the concrete surface erosion due to abrasions or to poor lining of the concrete surface finish (Khatsuria, 2005).

Hydraulic structures always run a risk of high probability for cavitation damages when they are in contact with high velocity flows. Therefore, factors have been established in order to know if a surface will be damaged or not (Falvey, 1990), and are as follows:

(24)

- Determining the cause of cavitation - Determining the location of the damage - Determining the intensity of the cavitation - Determining the flow velocity

- Determining the air concentration in the water - Determining the surface resistance to damage - Determining the exposure time of the surface

The cavitation index is also used to assess whether a flow is prone to cavitate. The index is dimensionless and is calculated according to the following equation:

𝜎𝑖 =𝑃_𝑜− 𝑃_𝑣 1 2 𝜌𝑉^𝑜²

where ꝍi is the cavitation index, P0 (Pa) is the static fluid pressure, Pv (Pa) is the vapour pressure, 𝜌 (kg/m³)is the density of the fluid and V0 (m/s) the velocity of the flow (Khatsuria, 2005). For the cavitation index ꝍi = 3 there is no cavitation, for ꝍi = 1.8 the cavitation is incipient, for 0.3

< ꝍi < 1.8 it’s a developed cavitation and for ꝍi < 0.3 it’s a supercavitation (Falvey, 1990).

As for the location of the cavitation damage, it always takes place downstream of the cavitation source, which is the collapsing cavitation cloud. Near the end of the cavitation cloud, it has been shown to be the region for the maximum damage. Moreover, with the increase of the discharge and surface irregularities’ height, the maximum damage also results in an increase.

However, for a cylinder, the cavitation damage takes place when the cavitation cloud’s length equals the cylinder diameter (Falvey, 1990).

Good examples of cavitation damages in spillway tunnels are the ones occurred in the spillway tunnels of the Hoover Dam and Glen Canyon Dam, both located in the USA. For the Hoover dam, reports show that cavitation occurred due to misalignment upstream of the damage location. The velocity of the flow at the time of the damage was recorded to be 45 m/s and the damage size was a hole 14 m deep, 35 m long and 9 m wide (Khatsuria, 2005). The Glen Canyon dam also experienced the same damages as the Hoover dam after the flooding of the Colorado river in 1983. The damage left the spillway with a hole as big as its diameter. The damage was reported to be due to residing cavitation damages on the concrete lining in several locations (Falvey, 1990).

(25)

9

Figure 2.1. Cavitation damage on the Hoover dam to the left and Glen Canyon dam to the right (Falvey, 1990).

To mitigate and avoid cavitation damages, aerators have been used in spillways that serve the purpose of inserting air into the water, especially close to the surface to minimize the damage risk (Yazdandoost and Attari, 2004, Falvey, 1990, Chanson, 1989). Moreover, studies have shown that for velocity of 12 m/s - 20 m/s cavitation damages may be avoided by concrete- lining, improving or eliminating irregularities in the surface and using better material for flows between the (Ruan et al., 2007, Chanson, 1989).

Air entrainment

As mentioned previously, aerators are commonly used to prevent cavitation damages. The use of aerators falls in the category of air-entrainment. Air entrainment is one of the best ways to avoid cavitation damages especially when the air is put as close as possible to the flow boundary (Khatsuria, 2005).

The phenomena of air entrainment is defined as the exchange of air contained within the atmosphere and water. Some other synonyms for air entrainment are air bubble entrainment or aeration. In addition, the air entrainment may occur from natural or artificial origins (Chanson, 1996). Natural air entrainment is referred to as self-aeration, which occurs when turbulence starts in a spillway and the turbulent boundary intersects the water surface and air entrains the bubbles in the turbulent boundary. For the artificial air entrainment is meant forced aeration by means of modification to the design, i.e. installing an aerator. Thorough studies should be carried out to decide the type and location of an aerator. The aerator is fed air through different mechanisms of air supply systems, such as the commonly used air-intake conduit or a duct

(26)

system. The air entrainment process can be seen when the water surface appearance turns from clear and glossy to irregular, white and bubbly (Khatsuria, 2005).

Air entrainment is also defined as the entrapment of un-dissolved air bubbles and air pockets through the flowing liquid. These concentrations can be classified locally (local aeration) or continuously (interfacial aeration) along the air-water flow (Chanson, 1996). The local air entrainment is a concentration of air bubbles located at the intersection of the impinging jet and at hydraulic jump. On the other hand, the interfacial aeration is defined as the air entrapped along an air-water interface, for example a chute (Yazdandoost and Attari, 2004).

In tunnel spillways air entrainment mechanism has differing characteristics depending on the flow condition, mainly partly full or pressurized. When the tunnel is partly full, the air entrainment mechanism is considered as an open channel flow (as described earlier). Here the main parameter to be considered in the design is the total air discharge, which consists of both the air flowing freely on the water surface and the air added to the flow through external means such as air vents (Khatsuria, 2005, Chanson, 1997). For pressurized or full flow tunnels air entrainment can cause serious damages of conveyance, unstable flow conditions and damage to the concrete lining. Moreover, depending on the design of the tunnel, the air might flow upwards instead and create air pockets on the water surface that will need to be released through air diffusing systems. Moreover, the downstream flow conditions in the shear layer are to be considered since they affect the air transport within the tunnel, and therefore, plays the ratio of length and diameter of the tunnel a role in the transport (Khatsuria, 2005).

Spillway tunnel aerators

As mentioned earlier, aerators are measures taken to induce forced aeration in order to prevent and mitigate cavitation damages. It was also mentioned that cavitation damage can be avoided by eliminating surface irregularities for flow with a velocity up to 20m/s, but for flow with a velocity over 30m/s aerators must be used to avoid cavitation damages (Ruan et al., 2007, Chanson, 1989).

There are various types of aerators, and the most basic devices are illustrated as: steps, deflectors, grooves, offset or a combination of them (Khatsuria, 2005, Falvey, 1990, Chanson, 1989, Volkart and Rutschmann, 1984). Figure 2.3 shows the illustrated aerators and their combinations.

(27)

11

Figure 2.3. Basic types of aerators. (Khatsuria, 2005, Falvey, 1990, Volkart and Rutschmann, 1984)

Some basic criteria to consider in design an aerator system are (Khatsuria, 2005):

- Positioning of the first aerator - Form of the aerator (type and size) - Quantity of entrained air by the aerator

- Form of the air supply mechanism (type and size) - Spacing between the aerators

To establish the location of the aerator, firstly the area where cavitation is deemed to be possible is the most likely primary parameter and is based on the cavitation index with a value of 0.2 or smaller. Afterwards the flow characteristics, depth and velocity, at the cavitation region and the curvature of the boundary should be taken in consideration (Khatsuria, 2005).

In tunnel spillways it’s hard to apply the above mentioned aerator types due to the provided space in tunnels because of their shape and the depth of the tunnel. And if applied the size of the aerators is usually limited because of the available space, but also because, in some examples such as the Glen Canyon dam, the criteria of the trajectory not hitting the tunnel roof to avoid blocking the flow is set (Khatsuria, 2005). Nevertheless, with careful and thorough design they all can be used.

Grooves are usually used in tunnels since they make the air supply through air vents easier.

Grooves are commonly combined with either an offset or a deflector to make the aeration mechanism more efficient (Volkart and Rutschmann, 1984). In this study, the discharge tunnel

(28)

at Zipingpu water conservancy has a shape of horseshoe and the installed aerator is as shown in Figure 2.3.

Numerical method

In hydrodynamics the study of fluid and flow characteristics and dynamics is an important step for the design of a hydraulic structure. To tackle the laws of physics and nature that surround fluid dynamics different approaches have been taken. Amongst these approaches, there are the practical and theoretical ones, where the first ones are dependent on experiments and the latter on relations between nature and mathematical equations (Griebel et al., 1997). In recent decades, thanks to the advancement in technology, a more efficient and cost-effective approach has evolved, the numerical method and simulation. The numerical method is a complement to mathematical analysis where complex equations of fluid dynamics cannot be solved. The detailed explanations of the equations and the theory behind computational fluid dynamics (CFD) will be discussed in chapter 3.

The numerical method is used to obtain an approximate numerical solution by discretization of the differential equation governing fluid dynamics. The numerical method is part of the computational fluid dynamics (CFD). CFD run methods have the capability of solving fluid equations in 2D and 3D (Ferziger et al., 2020). For the discharge tunnel at Zipingpu Water Conservancy project a 3D model will be simulated using CFD. The CFD software used for the simulation of the fluid flow is a software called ANSYS Fluent developed by ANSYS.

(29)

13

3 Theory

Mathematical model

In Newtonian fluid flows, Navier-Stokes equations and the governing equations of fluid mechanics are commonly used in mathematical modelling (Griebel et al., 1998). This is one of the mathematical fundamentals for the Computational Fluid Dynamics modelling.

For centuries, scientists have tried to explain the mechanism of the fluids through physics and mathematics. The Navier-Stokes equations describe the flow of Newtonian fluids accurately at a mathematical level analysis. For very low Reynolds numbers and simple geometries, it is possible to obtain complete explicit results (Wolfram, 2002).

These equations were derived from the conservation mass laws and momentum equations in the 1840’s (Wolfram, 2002) to a system of nonlinear partial equations with independent variables (Griebel et al., 1998). The equations include three different parameters described by the following concepts (Griebel et al., 1998) :

● 𝜐̅:Velocity field,

● 𝜌:pressure,

● 𝜚:density

The system of equations basically describes the advective-convective forces of the fluids (velocity field) and the external forces (i.e. pressure and viscosity) that oppose resistance. The term Re is the dimensionless Reynolds (Re) number, the 𝑔 term belongs to body forces like gravity acting throughout the bulk of the fluid and Δ represents a gradient differential operator (Griebel et al., 1998). The equations that are dimensionless can be given in their simplest form (where 𝑑𝑖𝑣 𝜐̅ is the divergence of velocity), as follows:

𝜕𝜌𝜐̅

𝜕𝑡 𝑑𝜐̅ + (𝜐̅ ∗ Δ)𝜐̅ + Δ𝜌 = 1

𝑅𝑒∗ Δ𝜐 + 𝑔 → 𝑀𝑜𝑚𝑒𝑛𝑡𝑢𝑚 𝑒𝑞𝑢𝑎𝑡𝑖𝑜𝑛 𝑑𝑖𝑣 𝜐̅: = 0 → 𝐶𝑜𝑛𝑡𝑖𝑛𝑢𝑖𝑡𝑦 𝑒𝑞𝑢𝑎𝑡𝑖𝑜𝑛

On the other hand, the Navier-Stokes equations are not suitable for this project due to its complexity regarding boundary conditions (model setup). Therefore, this mathematical model has to be adapted parallelly with numerical methods in order to solve the CFD model in this study case.

(30)

Discretization method (Finite Volume Method)

After choosing the mathematical model, a discretization approach is commonly used in order to calculate a fluid flow. This process consists of approximation of the differential equations by a system of algebraic equations (Ferziger et al., 2020). In this case study, ANSYS Fluent will be used for CFD modelling and it is based on the Finite Volume Method (FVM) approach (Jeong and Seong, 2014).

3.2.1 Finite Volume Method

The FVM subdivides into a finite number of control volumes (CVs) that makes it different from the Finite Difference Method. The FVM uses the integral form of the conservation equations as the starting point (Equation 3.3) (Ferziger et al., 2020). FVM solution depends of the integral from with respect to the conservation equations (Chakraverty, 2019):

∫ 𝜌𝜙 ∗ 𝑛 𝑑𝑆

𝑠

= ∫ ΓΔ𝜙 ∗ 𝑛 𝑑𝑆

𝑠

+ ∫ 𝑞_𝜙𝑑𝑉

𝑣

(𝐸𝑞. 3.3)

The surface integrals represent the convection (𝜌𝛷𝑉 ∗ 𝑛), diffusion (𝛤𝛥𝛷 ∗ 𝑛) and the flux vector (𝑞_𝜙𝑑𝑉) in each CV face (Ferziger et al., 2020). However, the velocity and the fluid properties are commonly known, but not the value for 𝛷 (i.e. dimensionless scalar value). In this project, the Φ value will be calculated by using an interpolation method explained in the next chapter. Also, in terms of transportation equations, an integration over the volume of a CV is required (𝑞_𝛷𝑑𝑉) (Ferziger et al., 2020). The converted expression (Equation 3.4) commonly used in the FV method can be written as:

𝜕𝜌𝜙

𝜕𝑡 + ∑ 𝜌_𝑓∗ 𝜐̅_𝑓∗ 𝜙_𝑓∗ 𝑆̅_𝑓 = ∑ Γ_𝜙𝑓 ∗ Δ𝜙_𝑓∗ 𝑆̅_𝑓

𝑁 𝑓𝑎𝑐𝑒𝑠

𝑓 𝑁 𝑓𝑎𝑐𝑒𝑠

𝑓

+ S_𝜙∗ 𝑉 (𝐸𝑞. 3.4)

One of the FVM approaches usually used for CFD to define CVs is by using a grid system (meshing) and computational nodes at the centre of the CVs. On the other hand, the computational node locations can be defined before the CVs (Ferziger et al., 2020). The choice between both approaches depends on the geometry of the domain.

3.2.2 Discretization methods

Even though there are approximations to the integrals, an interpolation method is needed to calculate the 𝛷 value. Usually, the commercial codes use different schemes, but also advice to choose an appropriate method for a situation in specific. As mentioned before, ANSYS Fluent is used to run the CFD model where the Upwind Interpolation (First order and Second Order)

(31)

15

and the Quadratic Upwind Interpolation (QUICK) methods can be used among others (Hamberg and Dahlin, 2019).

When a first-order accuracy is required, the CV faces are determined by assuming that the nodal center values represent a mean value throughout the CV (ANSYS, 2009). Thus, the value 𝛷f is equal to the CV center value of 𝛷 in the upstream CV. On the other hand, in the second order scheme the CV faces are calculated using a multidimensional linear reconstruction scheme.

Moreover, in the second order approach a higher accuracy is achieved at the CV faces through an expansion of Taylor series of the nodal-centred solution on the nodal center (ANSYS, 2009).

Therefore, the face value 𝛷f is computed using the following expression (Eq. 3.5):

Φ_𝑓= Φ + ΔΦ ∗ 𝑟⃗ (𝐸𝑞. 3.5)

where 𝑟 is the displacement vector from the centre of the upstream CV to the centroid of the face. It is important to mention that the term ∇φ requires a formulation approach in order to be solved. The formulations to determine the gradient of φ (∇φ) can be: Green-Gauss Cell-Based, Green-Gauss Node-Based and Least Squares Cell-Based.

The QUICK scheme is usually run in quadrilateral and hexahedral meshes, where upstream and downstream faces with unique characteristics can be identified. The QUICK method usually tends to be more accurate on structured meshes aligned with the flow direction. Also, this method is based on a weighted average of second-order upwind and central interpolations of the variable (ANSYS, 2009).

For this project, the governing equations must be discretized in both space and time (ANSYS, 2009). The temporal discretization implies the integration of every differential equation over a time step. The following expression represents the first-order discretization variable φ:

𝜙^𝑛+1− 𝜙^𝑛

Δ𝑡 = 𝐹(𝜙) (𝐸𝑞. 3.6) While the second order discretization:

3𝜙^𝑛+1− 4𝜙^𝑛+ 𝜙^𝑛−1

2Δ𝑡 = 𝐹(𝜙) (𝐸𝑞. 3.7)

Therefore, there are two approaches usually used in time discretization: Implicit and Explicit Time Integration. Firstly, the implicit method is to test the function F in a future time step (Equation 3.8), which can be solved by iterations at each time step before moving to the next step.

𝜙^𝑛+1− Φ^𝑛

Δ𝑡 = 𝐹(𝜙^𝑛+1) (𝐸𝑞. 3.8)

(32)

The term “implicit integration” is referred to a given cell 𝜙^𝑛+1 that is related to 𝜙^𝑛+1 in neighbouring cells through F(𝜙^𝑛+1) (Eq 3.9). The implicit equation is solved by performing iterations at each time level moving to the next time step (ANSYS, 2009):

𝜙^𝑛+1 = 𝜙^𝑛 + Δ𝑡 ∗ 𝐹(𝜙^𝑛+1) (𝐸𝑞. 3.9)

Secondly, the explicit method performs 𝐹(𝜙) at the current time level and it is available when the model used the density-based solver. Thus, this method is not used in this project because the usage of the pressure-based solver.

3.2.3 Pressure based solver

In order to solve the governing equations a flow solver must be used. There are two types of solvers that can be used in ANSYS Fluent: the density or pressure-based solver. In this project we are analysing a multiphase flow (gas-liquid), therefore the density-based solver can be neglected due to the different densities between the fluids. Nevertheless, both approaches use a similar discretization process (i.e. FVM), but when it comes to solve the discretized equations the scheme is different (ANSYS, 2009).

In difference to the density-based solver, the pressure-based solver uses a solution algorithm where the governing equations (i.e. mass conservation equations) are not linear but coupled to one another. The process involves a series of iterations wherein the set of equations is solved until the solution converges (ANSYS, 2009). The governing equations can be solved segregated or coupled from the other equations.

In the segregated algorithm, each governing equation (e.g. u, v, T, 𝑘, 𝜖) is solved step by step, one after another. Also, this method is memory-efficient, which means that the discretized equations are stored in the memory once at a time. Therefore, the convergence solution is slower than the coupled method (ANSYS, 2009). This algorithm solution (Figure 3.1) is illustrated in the following diagram:

(33)

17

Figure 3.1 Conceptual model of coupled algorithm of the pressure-based solver (ANSYS, 2009)

The coupled algorithm approach solves coupled systems by integrating the momentum equations and the pressure-based continuity equation. Besides, since the momentum equations are solved, the rate of solution convergence shows a better performance than the segregated approach, but its memory requirement increases by 1.5 to 2 times in comparison to the segregated approach.

3.2.4 Reynolds Averaged Navier -Stokes

There are different numerical approaches to describe turbulent flows, but the Reynolds- averaged Navier-Stokes (RANS) is a practical way which does not require as much computational calculations as some others, such as the Large Eddies Simulations (LES). Firstly, the RANS method was developed by Osborne Reynolds over a century where the governing equations were averaged over the spatial volume, not time. Usually, the RANS is averaged over the time (Ferziger et al., 2020).

(34)

In RANS, the starting point is the Reynolds decomposition of flow variables into mean and fluctuating parts. Then, the Reynold decomposed variables are inserted in the Navier Stokes equations, followed by an averaging of the equations involved (Alfonsi, 2009). For the velocity components (ANSYS, 2009):

𝜐_𝑖 = 𝜐̅_𝑖 + 𝜐′_𝑖 (𝐸𝑞. 3.10)

where 𝜐̅_𝑖 and 𝜐′_𝑖 are the mean and fluctuating velocity components. The same approach is used for the scalar quantities:

𝜙_𝑖 = 𝜙̅_𝑖 + 𝜙′_𝑖 (𝐸𝑞. 3.11)

where 𝜙 denotes a scalar variable like pressure, energy, or concentration of species. Thus, after substituting the expressions (number of the equations) into the continuity and momentum equations (General Navier Stokes equations) and taking a time average yields the RANS equations (ANSYS, 2009). They can be written as:

𝜕

𝜕𝑡𝜌 + 𝜕

𝜕𝑥_𝑖(𝜌𝜐_𝑖) = 0 (𝐸𝑞. 3.12)

𝜕

𝜕𝑡(𝜌𝜐̅_𝑖) + 𝜕

𝜕𝑥_𝑗(𝜌𝜐_𝑖𝜐_𝑗)

= −𝜕𝑝

𝜕𝑥_𝑖+ 𝜕

𝜕𝑥_𝑗[𝜇 (𝜕𝜐_𝑖

𝜕𝑥_𝑗+𝜕𝜐_𝑖

𝜕𝑥_𝑖−2 3𝛿_𝑖𝑗𝜕𝜐_𝑖

𝜕𝑥_𝑖)] +𝜕𝜐_𝑖

𝜕𝑥_𝑗(−𝜐̅̅̅̅̅̅̅) (𝐸𝑞. 3.13) ^′_𝑖𝜐^′_𝑗

Multiphase flow model

Often in engineering, the multiphase flow models are used for different applications such as combustion systems. Multiphase flows with a phase change could be used to assess different parameters such as cavitation, melting, solidification, boiling and condensation (Ferziger et al., 2020). In the present model, a multiphase flow approach has to be chosen due to the water and air being the main fluxes to test cavitation within the hydraulic structure (tunnel).

In order to solve any multiphase flow, as the present project, it is important to choose an appropriate model. There are two numerical schemes for multiphase flow: the Euler-Lagrange method and Euler-Euler method, but in this case, the second one (Euler-Euler method) was chosen due to different reasons that are explained in this chapter (ANSYS, 2009).

In ANSYS Fluent, there are different Eulerian-Eulerian methods such as Volume of Fluids (VOF) model, Mixture model, Eulerian Model. The VOF model was chosen mainly due to its suitability for free-surface flows which is the case of this simulation. Moreover, VOF can only be performed with the pressure-based solver.

(35)

19

Turbulence model

To solve and close the equations (RANS) we must include a turbulence model. Most of the studied flows in engineering are turbulent and therefore require a different numerical approach compared to laminar flows. This project has a turbulent flow containing velocity fluctuations that mix transport properties such as momentum, energy, and species of concentration, but also cause the transported properties to fluctuate (ANSYS, 2009).

For many years, the studies on turbulence models were merely experimental, but the costs and time required for this approach were high in comparison with numerical simulations. Currently, the most accurate method uses the direct numerical simulation (DNS) in which the Navier- Stokes equations are solved to calculate the motions in a turbulent flow (Ferziger, 2020).

The DNS is also the simplest approach from the conceptual point of view. The obtained information of a DNS contains detailed characteristics of the flow. In other words, the results are like the experimental approaches and can be used to create statistical information or a

“numerical flow visualization”. This data can be used to acquire deeper knowledge of the physics of the flow or to build a quantitative model.

3.4.1 𝒌 − 𝝐 Turbulent model

According to ANSYS (2009), the simplest turbulence simulations are given by the two-equation models, where the solution can be independently solved in terms of the turbulent velocity and length scales. The 𝑘 − 𝜖 builds up the length and time scales from the turbulent kinetic energy 𝑘 and the turbulent dissipation rate 𝜖 (Alfonsi, 2009). This model is suitable for 3D modelling.

The CFD software (ANSYS Fluent) used for this project includes three different turbulence models: Standard, RNG and Realizable. The three models are similar, but there are some differences explained herein in terms of:

• Calculating turbulent viscosity

• The Prandtl number that governs the turbulent diffusion

• The generation and destruction in the 𝜖 equation (ANSYS, 2009)

The Standard 𝑘 − 𝜖 method is the simplest one of the three mentioned and it is the one used for this CFD model. This model assumes that the flow is fully turbulent, and the effects of molecular viscosity are negligible (ANSYS, 2009). The components (𝑘, 𝜖) in the Standard scheme are obtained from the following equations:

𝜕

𝜕𝑡(𝜌𝑘) + 𝜕

𝜕𝑥_𝑖(𝜌𝑘𝜐_𝑖) = 𝜕

𝜕𝑥_𝑗[(𝜇 +𝜇_𝑡 𝜎_𝑘) 𝜕𝑘

𝜕𝑥_𝑗] + 𝐺_𝑘+ 𝐺_𝑏− 𝜌𝜖 − 𝑌_𝑀+ 𝑆_𝑘 (𝐸𝑞. 3.14)

(36)

and

𝜕

𝜕𝑡(𝜌𝜖) + 𝜕

𝜕𝑥_𝑖(𝜌𝜖𝜐_𝑖) = 𝜕

𝜕𝑥_𝑗[(𝜇 +𝜇_𝑡 𝜎_𝜖) 𝜕𝜖

𝜕𝑥_𝑗] + 𝐺_1𝜖𝜖

𝜅(𝐺_𝑘+ 𝐺_3𝜖𝐺_𝑏) − 𝐶_2𝜖𝜌𝜖²

𝑘 + 𝑆_𝜖(𝐸𝑞. 3.15) where 𝐺_𝑘 represents the generation of turbulent kinetic energy due to the mean velocity gradients, Gb is the generation of turbulence due to buoyancy and 𝑌_𝑀 represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate. The 𝐺_𝜖’s components are constants and the 𝜎_𝜖and 𝜎_𝑘 are the Pradntl number for 𝑘 and 𝝐, respectively.

Finally, 𝑆_𝑘 and 𝑆_𝜖 are user-defined terms (ANSYS, 2020).

In parallel, the turbulent viscosity 𝜇_𝑡, is computed by the following expression:

𝜇_𝑡 = 𝜌𝐶_𝜇𝑘²

𝜖 (𝐸𝑞. 3.16) where 𝐶_𝜇 is a constant (ANSYS, 2020).

Boundary conditions

3.5.1 Inlet and outlet

For FVM analysis, the boundary conditions require that the fluxes involved must be known in terms of quantities and nodal values. CFD models have many boundary conditions that allows the fluid(s) to enter and exit the flow domain. There are different boundary conditions that are used for different CFD models, but for this project the following list shows the specific ones for the case of the discharge tunnel of Zipingpu Water Conservancy project:

● Velocity inlet value is used to define the velocity and scalar properties of the flow at inlet boundaries.

● Pressure inlet value is used to define the total pressure and other scalar quantities at flow inlets.

● Outflow boundary values, used to model flow exits where the details of the flow velocity and pressure are not known prior to solution of the flow problem (ANSYS, 2009).

3.5.2 Wall function boundar y

In CFD models, the wall boundaries are used to bound fluid and solid regions. In ANSYS Fluent, the no-slip function is given at walls by default, which was also the one chosen for this project. Moreover, it indicates that the fluid sticks to the wall and moves with the same velocity as the wall, only if it is moving.

(37)

21

According to ANSYS (2009), the following information is required for a wall boundary:

● Wall motion conditions (for moving or rotating walls)

● Shear conditions (only for slip walls, optional)

● Wall roughness (for turbulent flows, optional)

● Thermal boundary conditions (for heat transfer calculations)

● Discrete phase boundary conditions

● Wall adhesion contact angle

For this project, a stationary wall with No-slip shear condition are chosen because both comply with the real conditions of the hydraulic structure (tunnel).

Meshing

As mentioned earlier the Navier-Stokes equation is hard to solve analytically due to its complexity and therefore a numerical method will be used in this study. In solving a numerical method, the geometry of the studied structure plays a crucial role in the choice of methods, computation time and cost. Therefore, for the geometry a proper mesh/grid has to be applied in order to smooth out the model computation and get accurate results.

3.6.1 Choice of mesh

The geometry is divided into subdomains in order to carry out the calculation and these subdomains represent the mesh or grid that is chosen. Some meshing options are structured, block-structured, and unstructured. Moreover, the beforehand chosen discretization method is deciding for the shape of the mesh, e.g. if the algorithm in the discretization method is set to simulate orthogonal grids then non-orthogonal grids cannot be used; or if the CV is needed to be hexahedral or quadrilateral then triangles or tetrahedral cannot be used. To be noted is that the accuracy and quality of the mesh is higher if the CV is hexahedral in 3D and quadrilateral in 2D (Ferziger et al., 2020).

Depending on the geometry in hand, the mesh can be carried out differently. For instance, when the geometry is simple it is relatively easy to choose a mesh type. But when the geometry becomes more complicated different approaches can be carried out to ensure the accuracy and quality of the mesh. Structured or block-structured meshing is compatible with simple geometries, the problem with these techniques is that the mesh accuracy deteriorates near the wall and therefore an overlapping scheme is needed. Overlapping grids give a better accuracy and mesh quality in simulations. Moreover, the accuracy of the mesh can be obtained by having smaller mesh sizes where the geometry and the fluid characteristics are complex. Nevertheless, it is not advised to have small mesh sizes throughout the entire geometry domain, given the fact

(38)

that the smaller the mesh size the longer the computational time will be and as a result will affect the cost (Ferziger et al., 2020).

When computed profiles of a certain variable are presented, it is recommended that numerical uncertainty be indicated by error bars on the profile, analogous to the experimental uncertainty.

It is further recommended that this be done using the GCI in conjunction with an average value of 𝑝 is equal to pave as a measure of the global order of accuracy.

3.6.2 Estimation of discretization error

As mentioned earlier a discretization method is needed in order to solve the differential equation through approximations. These approximations commonly result in errors that are defined as the difference between the solution of the discrete approximation and the solution of the governing equation prior to the discretization (Ferziger et al., 2020).

The mesh size plays an important role in the weight of the obtained error in discretization. The recommended method in estimating the discretization error is the GCI (Grid Convergence Index) that is based on the RE method (Richardson Extrapolation) (Celik et al., 2008). After having set the size of the different meshes, i.e. coarse-, medium- and fine mesh (h1, h2, h3), the GCI is calculated according to the equations below. The apparent order p is calculated using Equation 3.17, as follows:

𝑝 = 1

ln(𝑟₂₁)|𝑙𝑛 |𝜀₃₂

𝜀₂₁|| + 𝑞(𝑝) (𝐸𝑞. 3.17.1)

𝑞(𝑝) = ln (𝑟₂₁^𝑃 − 𝑠

𝑟₃₂^𝑃 − 𝑠) (𝐸𝑞. 3.17.2) 𝑠 = 1 ∗ sgn (𝜀₃₂

𝜀₂₁) (𝐸𝑞. 3.17.3)

where r32 = h3 /h2, r21 = h2 /h1, r is the mesh ratio and it is preferred if the mesh refinement is structured even though the mesh itself is unstructured, i.e. same decrease or increase between the different mesh sizes. 𝜀32 = 𝜙3 - 𝜙2,𝜀21 = 𝜙2 - 𝜙1 and 𝜙 is the dimensionless scalar value obtained from the interpolation method, as mentioned earlier. Thereafter the relative error is estimated using Equation 3.18 and the GCI is calculated using Equation 3.20.

𝑒_𝑎²¹ = |𝜙₁− 𝜙₂

𝜙₁ | (𝐸𝑞. 3.18)

When the used mesh refinement ratio is 1 and it is constant other equations prevail. Therefore instead of using Equation 3.17.2, a simplified version of the equation can be used to calculate the order p for a three-grid solution (Eq. 3.19.1), even the equation for the relative error changes (Eq. 3.19.2) (Wagner et al., 2002, Roache, 1998). These are calculated using the following:

(39)

23 𝑞(𝑝) =

ln (𝑓₃ − 𝑓₂ 𝑓₂− 𝑓₁)

ln(𝑟) (𝐸𝑞. 3.19.1) 𝑒_𝑎²¹= |𝑓₁− 𝑓₂

𝑓₁ | (𝐸𝑞. 3.19.2)

where f1, f2 and f3 are the different grid solution. The GCI for fine mesh is then obtained by Equation 3.20,

𝐺𝐶𝐼_{𝑓𝑖𝑛𝑒}²¹ = 1.25𝑟_𝑎²¹

𝑟₂₁^𝑃 − 1 (𝐸𝑞. 3.20)

where 1.25 is a safety factor used for three-grid solution, but the value of the safety factor can go up to 3, which is considered moderate and mostly recommended for two-grid solutions (Roache, 1998; Wagner et al., 2000). The obtained GCI should be low in order to achieve a grid independent result.

(40)

(41)

25

4 Methodology

In this chapter, the undertaken steps in the methodology in order to conclude the study are represented. The steps will cover the building of the geometry, meshing scheme, model setup, simulation, post-processing of the results, validation of the model and lastly evaluating different scenarios of aerators and discharge rate dependent on the achieved results from the initial simulations.

The ANSYS version 2019 R3 workbench offers different tools to perform the CFD model, from the geometry construction to the post-processing stage. The used softwares in the process were as follows:

- ANSYS Spaceclaim was used in the building of the geometry - ANSYS mesh generator was used for meshing the geometry - ANSYS Fluent was used for the simulation

- Excel was used to analyze the results from the simulations in order to post-process and validate the results.

Geometry

The geometry was built as close as possible to the provided design drawings of the discharge tunnel. Some details were neglected in order to simplify the meshing and simulation process, a detailed illustration will follow. The design drawings of the tunnel are illustrated below, where Figure 4.1 represents the side view of the entire tunnel excluding the tunnel inlet, Figure 4.2 represents the side view of the tunnel inlet and Figure 4.3 the plan view of the tunnel inlet.

In Figure 4.1, all the measurements for the elevations and length distances are given in meters.

Since in this figure the inlet is excluded the tunnel length starts at 0+28.793 m and ends at 0+583.00 m. After the end point of the tunnel, an extension of the tunnel is included as it is seen in the figure. The side view of the inlet in Figure 4.2 has its measurements given in cm and the elevations are given in m. As it is the start of the tunnel the length starts from 0+0.000 m and the connection to the rest of the tunnel through a floodgate at 0+028.793 m. The floodgate is marked with a red circle since it was not included in the simulated geometry. The plan view of the tunnel inlet which is represented in Figure 4.3 has its length measurements given in cm and the red circle marks an extension that was excluded from the simulated geometry.

(42)

Figure 4.1. Side view of the entire tunnel.

(43)

27

Figure 4.2. Tunnel inlet side view

Figure 4.3. Tunnel inlet top view.

As mentioned earlier the simulated geometry was designed as identical to the provided design drawings (Figure 4.1-4.3). Some detailed characteristics of the structure have been excluded, such as the above mentioned floodgate in Figure 4.2 and the tunnel inlet extension in Figure

Aeration and risk mitigation for flood discharge tunnel in Zipingpu water conservancy project

Aeration and Risk Mitigation for Flood Discharge Tunnel in

Zipingpu Water Conservancy Project

JORGE CONTRERAS MORENO

KIBRET DAWIT GHEBREIGZIABHER

Aeration and Risk Mitigation for Flood Discharge Tunnel in

Zipingpu Water Conservancy Project

J ORGE C ONTRERAS M ORENO

K IBRET D AWIT G HEBREIGZIABHER

Master of Science Thesis

Stockholm, Sweden 2020

Abstract

Sammanfattning

Preface

Contents

List of acronyms

1 Introduction

Aims and objectives

Limitations

2 Background

Water scarcity and water conservancy projects

Flood discharge tunnel

General description of Cavitation

Cavitation damages

Air entrainment

Spillway tunnel aerators

Numerical method

3 Theory

Mathematical model

Discretization method (Finite Volume Method)

3.2.1 Finite Volume Method

3.2.2 Discretization methods

3.2.3 Pressure based solver

3.2.4 Reynolds Averaged Navier -Stokes

Multiphase flow model

Turbulence model

3.4.1 𝒌 − 𝝐 Turbulent model

Boundary conditions

3.5.1 Inlet and outlet

3.5.2 Wall function boundar y

Meshing

3.6.1 Choice of mesh

3.6.2 Estimation of discretization error

4 Methodology

Geometry