SIMULATION OF TURBULENT FLOW THROUGH TARBELA DAM TUNNEL 3

Tarbela dam is one of the largest earth filled dam in the world. The inflow of sediments in the Tarbela reservoir has resulted in reduction in water storage capacity. During recent years, a reasonable increase of sedim ent particles in the tunnel has been observed. This is damaging to the tunnels, power ge nerating units and is a severe threat to the plant equipment. To the authors knowledge, t o-date no comprehensive simulation studies have been performed for flooding in the res ervoir or turbulent flows in the tunnels. In this paper, turbulent flow in Tunnel 3 of the Tarbela Dam is analyzed using Reynolds Stress Model with and without considering the effect of sediments particle. Results are presented for three different water hea ds in the reservoir i.e. considering summer, winter and average seasons and for one-way and two-way/full coupling for sediments particle tracking/deposition. The effect of cavitation erosion and damage to the tunnels due to erosion is investigated and resu lts are compared with the experimental erosion results for similar geometries and are foun d to be in good agreement. Sediments particulate analysis is also performed for the vali dation of the samples collected from WAPDA (water and power development authority). More ov r, pressure, velocity and erosion rate results are discussed to get complete behavior of the turbulent flow of water in the tunnel.


INTRODUCTION
Tarbela Dam comprises of six tunnels, three of which are used for power generation and three for irrigation purposes.Moreover there are four outlets of the tunnel used for power generation and four pressure relief valves.Further details of this tunnel are given in Table 1 [1,2].During recent years, a reasonable increase of sediment particles in the tunnel has been observed.This is damaging to the tunnels, power generating units and is a severe threat to the plant equipment.To the authors knowledge, to-date no comprehensive simulation studies have been performed for flooding in the reservoir or turbulent flows in the tunnels.In this paper, turbulent flow in Tunnel 3 of the Tarbela Dam is analyzed using Reynolds Stress Model with and without considering the effect of sediments particle.Results are presented for three different water heads in the reservoir i.e. considering summer, winter and average seasons and for one-way and two-way/full coupling for sediments particle tracking/deposition.The effect of cavitation erosion and damage to the tunnels due to erosion is investigated and results are compared with the experimental erosion results for similar geometries and are found to be in good agreement.The studies performed in this paper are carried out on a large structure, the tunnel lengths and diameters are in meters.It makes this research work unique and different from previous studies cited as references in the paper.In addition, Sediments particulate analysis is also performed for the validation of the samples collected from WAPDA (water and power development authority).Moreover, pressure, velocity and erosion rate results are discussed to get complete behavior of the turbulent flow of water in the tunnel.In the related work, Alamgir et al. [3] have investigated the particle deposition and suspension in a horizontal pipe flow.The deposition was studied as a function of particle diameter, density and velocity of fluid.The lighter particles were found to remain suspended with homogeneous distribution.The larger particles clearly showed deposition near the bottom of the wall.In our work the diameter of the flow passage is very large, the results show the dependance of velocity, pressure and erosion rate density on the passage diameter.Xianghui et al. [4] presented a computational fluid dynamics (CFD)-based erosion prediction model and its application to oilfield geometries specifically elbows and plugged tee geometries.This comprehensive procedure consists of three major components: flow simulation, particle tracking, and erosion calculation.The analysis procedure is taken from this study but cavitation analysis procedure is different from this.Gary [6] explains the Lagrangian approach for particle tracking, authors have used the same approach in one way coupling but Eularian approach is used for two way coupling.Hari [7] explains the role of different forces when a solid particle passes through a fluid, but in our case the rotational force is ignored.[3].
For analysis of turbulent flow in the tunnel, we divide variables into time-averaged part and fluctuating part.Governing equations of continuity and momentum for flow employed in CFX-11 are given in equations 1 and 2 respectively [4].

( )
. 0 Here, u u ′ ′ ρ ⊗ is the Reynolds stress; and stress tensor σ is given by: ( ) For sediments particle deposition in turbulent flow, we use Lagrangian particle transport and Eulerian-Eulerian multiphase approaches [5], using a RANS framework.
One-way and two-way/full coupling options are used depending upon the value of β, which is defined as the ratio of the particulate mass per unit volume flow to the fluid mass per unit volume flow and is given by; where, β = 0.2 is the threshold value [6].One-way coupling is valid for volume concentration up to 14.86% and simply predicts particle paths during post-processing based on the flow field without affecting the flow field (i.e.particles are assumed not to interact with each other).In two way/full coupling particles exchange momentum with continuous phase, allowing the continuous flow to affect the particles and vice versa.
Generally, a small, rigid spherical particle entrained in the turbulent pipe flow encounters many forces [7], such as drag and lift forces, buoyancy force, rotation force and turbulence force etc and has the following general form (5) where, F x represents additional forces in the particle force balance that can be used under special circumstances.Moreover, g x is the force of gravity.For present study we have ignored the force of gravity on the particle, it is important to note that for most cases the gravitational acceleration is taken zero.One can include the gravity force, but should define the magnitude and direction of the gravity vector.Virtual mass force which is used to accelerate the fluid surrounding the particle is also ignored alongwith additional forces arise due to the pressure gradient in the fluid and due to rotation of the reference frame.Later forces are important only when fluid is flowing in a rotating frame.Another force which is not related for our model is produced because of temperature gradient effecting small particles suspended in a gas.This phenomenon is known as thermophoresis.For micro particles, the effects of Brownian motion is optional and if required can be included in the additional force term.Therefore, the governing particle equation of motion takes the following form: ( ) where drag force per unit mass is expressed as: with τ p as the particle response time and is given by: Furthermore Re p is the particle Reynolds number based on the relative velocity between the particles and carrier phase given by ( ) Moreover C D is the drag coefficient and is used in CFX by the Schiller Naumann correlation.Schiller Naumann correlation is derived for flow of a single spherical particle and is valid in the dilute limit of small solid phase volume fractions for one-way coupling with a coefficient of 0.1.Particle transport drag coefficient correlation is used in twoway/full coupling phenomena with a coefficient of 0.44 [8].

Erosion Estimation and Surface Damage
Sediment erosion phenomenon is highly complicated and a wide range of factors contribute to erosion severity [9].Analysis is done using ANSYS CFX to determine the erosion rate under turbulent flow in the tunnel for different heads, and variation of particles concentration.In ANSYS CFX, only Finnie and Tabakoff erosion models are available.Tabakoff model provides more scope for customization with its larger number of input parameters but presence of large number of coefficients in the model can become a source of errors in the numerical scheme [10].We therefore use simplified erosion model of Finnie with Langrangian particle tracking and Eulerian-Eulerian multiphase approaches.
We have used the following expression for erosion rate density: In ANSYS CFX, implementation of overall erosion rate at each point on the surface is calculated by multiplying E with mass flow carried by the Langrangian particle impacting the surface, and then summing for all particles.This ultimately leads to an erosion rate density variable.Deformation wear occurs when repeated particle impacts at high impact angles and plastically deform surface layers of the material, eventually causing material loss through surface fragmentation.Cutting wear occurs due to particle impacts at small angles, with a scratch or cut being formed on the surface if the shear strength of the material is exceeded.The other critical factor affecting wear is the particle impact velocity, with both cutting and deformation wear being proportional to impact velocity raised to a power n determined through physical tests.In general n varies between 2 and 3 depending on both the surface and particle materials.The value of k is set to 1 which has been determined experimentally [11].
The total erosion rate at a particular point on a surface is measured by summing the contributions due to the deformation and cutting mechanisms and depends on the properties of the material, with deformation wear being more significant for hard or brittle materials and cutting wear being more significant for softer or ductile materials.For standard commercial grade steels, as used in most of the bend surfaces, peak erosion rates have been measured to occur at impact angles of 25-30° [12], indicating that cutting wear dominates.

Cavitation Erosion
Cavitation is a phenomenon of formation of vapor bubbles in low pressure regions and collapse in high pressure regions.High pressure is produced and metallic surfaces are subjected to high local stresses.It is a phenomenon which manifests itself in the pitting of the metallic surfaces of tunnel parts because of the formation of cavities.The most critical location is the S-bend at the middle of the tunnel where pressure falls well below the vapour pressure i.e; 750 Pa in one way coupling and 450 Pa in two way/full coupling.The tendency for a flow to cavitate is characterized by the pressure coefficient or cavitation number, given by [13]: ( )

Cavitation Models
Different models have been proposed for cavitation phenomenon in literature, but here we consider the followings: The Rayleigh-Plesset equation provides the basis for the rate equation controlling vapour generation and condensation.The Rayleigh-Plesset equation describing the growth of a gas bubble from a liquid is given by: The Homogeneous model can be viewed as a limiting case of Eulerian-Eulerian multiphase flow in which the inter phase transfer rate is very large.This model is used in highly dispersed flows with one gas phase and the other dispersed liquid droplet phase.
For the present study we use Rayleigh Plesset cavitation model which requires the model parameters given in Table 3 [14].

MODELING AND ANALYSIS
Modeling of tunnel is done in Pro-Engineer software as per dimensions taken from [15] and is shown in Fig. 1.Tunnel model is then meshed in ICEM CFX with free mesh option using 1843803 tetrahedral elements as shown in Fig. 2.  Meshed model is imported into ANSYS CFX for detailed analysis.The particles are assumed to be randomly distributed at the inlet.Particles injected at the inlet are proportional to the mass flow rate of water flowing into the tunnel.The sediment particles volume fraction is only 0.007% at the high head during months of July till September, which increases to 6.1% at the minimum head level during March till June.These fractions fall in the case of one-way coupling.The samples collected and analyzed in February 2010 show that the concentration of sediment particle has increased to 19.6% which lies in the two-way/full coupling region as it crosses the threshold value of 14.86% [6].The analysis reveals that the continuous flow changes to dispersed flow and flow field is now affected by the sediment particles, which is more critical for the damage caused by erosion.Standard no-slip wall functions were applied at all solid surfaces for the fluid phase.Furthermore, coefficient of restitution for the particles were taken as 0.9 for parallel flow direction and 1.0 for perpendicular flow direction.
The effect of cavitation erosion is analyzed at the S-bend of tunnel of 200 m length.The velocity boundary conditions of 38.51 m/s and 32.09 m/s are specified at the inlet and outlet of the S-bend portion.
During analysis, pressure relief valves are excluded in the geometry which might affect the water velocities and pressures at different locations of the tunnel.Boundary and initial conditions applied are listed in Table 4. Zero pressure is specified at the tunnel outlet being exposed to the atmosphere.The list of other input parameters is given in Table 5. Navier-Stokes equations for mass, momentum and fluid turbulence are solved with the commercial code CFX-11 using the finite volume technique.Convection terms in momentum equations are discretized using a simple high resolution advection scheme.

Mesh Sensitivity Analysis
The purpose of mesh senstivity analysis is to determine the minimum grid resolution required to generate a solution that is independent of the grid used.Starting with a coarse grid the number of cells was increased in the region of interest until the solution from each grid was unchanged for successive grid refinements.
Mesh sensitivity analysis was performed using the following methodology: 1. Meshes of four different element sizes (2, 3, 5 and 7) were generated using ANSYS ICEM. 2. Boundary conditions were applied to all these four meshes.3. Transient analysis was performed for all the four meshes.4. Final residuals were compared for differences in comparison for changes in mesh sizes.
All the results were compared with a mesh size of 3 as the standard.We have evaluated the effect of further refining the mesh to a mesh size of 2 and also coarsening it to a mesh size of 5.The mesh senstivity statistics are tabulated in Table 6.
Results of important variables such as velocity, pressure and erosion rate density were also taken at a particular point in the tunnel defined by the coordinates (793, 1, and -44) and are given in Table 7.

Equation Residual
If the solution is "exact" then the residual is zero.Exact means that each of the relevant finite volume equations is satisfied precisely.The "residual" of an equation identifies by how much the left-hand side of the equation differs from the right-hand side at any point in space i.e; Ax-B = residual.
The RMS residual is obtained by taking all of the residuals throughout the domain, squaring them, taking the mean, and then taking the square root of the mean.The peak values of residues, final residues as summarized in Table 8, convergence plots and important variables are compared for analysis performed on different mesh sizes for Tarbela Dam Tunnels.

Convergence Plots
Turbulence, mass and momentum plots for the residuals are shown in Fig. 3 and Fig. 4. The Mass and Momentum plot of the mesh size of 2 shows smoother convergences as compared to the other mesh however none of the meshes have shown large peaks of the residuals.This must be kept in mind that the mesh size of 2 takes far more computational power than the rest and still does not substantially improves the results.
The turbulence plot of the mesh size of 2 again shows smoother convergences when compared to other meshes.This could be because of the finer mesh but the residual peak values do not justify the amount of computational time and power required for such fine mesh.

RESULTS AND DISCUSSION
The following sections are dedicated for complete analysis of results obtained using ANSYS CFX for velocity, pressure and erosion rate density for different coupling techniques and water heads.Specific attention is paid to critical region such as S-bend.

One-way Coupling at All Heads with Sediment Particles
For one way coupling maximum velocities of water 28.89, 41.73 and 51.34 m/s are measured at low, medium and high head of water respectively at the S-bend.Velocities are 7.23, 10.44 and 12.84 m/s at the inlet section.Velocity increases to 21.67, 26.40 and 32.09 m/s when the water flow is fully developed at 150 m from the vertical section at low, medium and high water head respectively.Velocity increases abruptly at the outlet branches due to reduction in their crossectional area.Moreover, maximum pressures of water 346.70,728.10 and 1103 kPa are measured at the inlet section where the velocity was minimum at low, medium and high water heads respectively.Minimum pressures 1.75, 4.50 and 0.75 kPa are measured at the S-bend where the velocity has it highest value respectively at low, medium and high water heads.Pressure decreases abruptly at the main branch and at the outlet branches due to increase in the velocity.Furthermore, higher erosion rate densities 11.75 x10 -5 , 8.13x10 -5 , 7.13x10 -5 kgs -1 m -2 due to the sediment particles are measured at high, medium and low heads of water respectively at the inlet section are measured.The highest value of the erosion rate density is measured as 23.31x10 -5 kgs -1 m -2 at the S-bend.It changes abruptly at the main branch and at the outlet branches due to the higher impact velocity and impact angle at these locations.

Effect of Coupling at High Head
In one way coupling increase in sediment particles concentration does not effect the velocity, pressure and erosion rate density, but in two way/full coupling the velocity of water deccreases to about 5% and pressure of water increases to 1% at high head.The small change in velocity and pressure values is due to the fact that the diameter of the tunnel is too large to effect the flow field completely.But the erosion rate density of tunnel wall increases to about 35% because of the higher number of sediment particles striking the tunnel walls with higher impact velocity at high head.The results are shown in Fig. 5    The velocity, pressure and erosion rate density values will be effected significantly when the sediment particles concentration reaches to the threshold value of 14.86% as the sediment delta is moving towards the tunnel.This shows that the coupling phenomena is a function of tunnel diameter.

Effect of Cavitation at High Head at S-bend
The cavitation phenomenon is observed at the S-bend for tunnel 3 due to the lower local pressure then the vapour pressure at this critical location.The local pressure at Sbend is 3.57 kPa which is less than the vapour pressure of water i.e; 4.24 kPa as shown in Fig. 11(a), therefore cavitation phenomena starts at this location.The water vapour volume fraction profile is shown in Fig. 11(b).The highest fraction of 0.795 is found at the inner periphery of the curvature of the S-bend.The low pressure at this location generates the bubble formation, when these bubbles enter the higher pressure region, they explode.As a result these bubble bursts causes the cavitaion erosion.The velocity profiles for water flow are discontinuous, i.e, water contents dominates the vapor contents at some portions and vice versa due to the mixing of water and water vapour phases (two phase flow) at S-bend as shown in Fig. 12.

Effect of no Sedimentation at High Head
Negligible increase of less than 2% in velocity and pressure is observed for the sediments particle flow.This demonstrates that the flow field is uneffected in one way coupling.Erosion in the tunnel for flow without the sediment particles is concluded due to the turbulent eddies [16] present at the walls of the tunnels at critical locations.

Experimental Set-up and Procedure
An experimental setup is developed to validate the numerical results discussed in section 5.2.The experimental setup is as shown in Fig. 15.The pipe loop is constructed in the horizontal plane with a valve to allow flow to be diverted to another loop as necessary.A high power stirrer is installed to help distribute the sediments in the tank.Pipe components are made of AISI 304L stainless steel with a nominal wall thickness of 3 mm.The geometric components like straight portion, bend-section and T-section are all analysed for their flow characteristics in this study.Pipe sections are prepared for weighing by firstly thoroughly rinsing with water to remove any sediments and then cleaned with warm 5% citric acid to remove calcite deposits on internal surfaces.This was necessary to allow the change in mass of the pipe sections to be attributed solely to erosion.The pipes are then allowed to air dry, usually overnight, prior to weighing.The experiment is performed for continuous flow of sediments for 80 hours with sediments concentration of about 0.4% by volume and a velocity of 0.175 m/s.The Reynolds number calculated is 2997 for this flow.The loops are then carefully dismantled for cleaning and weighing.

Numerical Solution of the Experimental Set-up
Pro-Engineer Wildfire 4.0 is used for the modeling of the components, i.e, straight portion, bend-section and T-section The ANSYS ICEM package is used for the meshing of the geometry.The number of elements used in the geometric components are 1500, 130 and 160 respectively.The analysis is done in ANSYS CFX as shown in Fig. 15.Table 12 shows the input parameters used in the software.Specifications for different regions of Tunnel 3 are summarized in Table 13.
The comparisons between the experimental and numerical results are tabulated in Table 14.The results show an error of about 8% for erosion rate density.

CONCLUSION
A CFD-based velocity, pressure and erosion prediction procedure is presented, based on one way and two way/full coupling phenomena.We have found that erosion is high at high head because of higher impact velocity of the sediment particles during July, August and September.Moreover, the results are shown for nine different locations of critical importance with and without sediment particles flowing through the tunnel where the velocities and pressures may vary causing the erosive damage.Interestingly, the erosion rate density is found to be maximum at the S-bend due to the additive effect of the cavitation erosion.Also, higher values of erosion are observed at the outlet branches due to several reasons like the higher velocity and impact angle and the production of turbulent eddies.Furthermore, the cavitation phenomenon is observed at the S-bend due to the lower local pressure then the vapor pressure at this critical location.Numerical simulations as well as experimental erosion tests are performed.Comparisons show that the CFD-based erosion prediction procedure is able to reasonably predict the erosion profile and satisfactorily capture the trend of erosion with respect to the carrier velocity with a maximum error of 8.5%.Simulation results show that particle volume fraction increases rapidly as the sediment delta moving towards the tunnels.The particle concentration is found to be 21% by mass and 21% by volume.This increase in particle concentration is changing continuous flow to dispersed flow and is damaging the tunnel walls.The particulate analysis is shown to have a good agreement with the data provided by WAPDA (water and power development authority).The error in the average diameter of the particles found to be only 7.8%.
-7 for one way coupling and Fig. 8-10 for two way/full coupling.

13. THEORY 2.1 Turbulent Flow and Particle Tracking For
flows with large, rapid, extra strains and flows with strong acceleration or retardation in the tunnel, we use second-moment closure model or Reynolds stress model (RSM) as it can reproduce anisotropy of the flow in the turbulent boundary layer.Coefficients of this model are given in Table 2 [3].

Table 4 :
Boundary conditions and Initial conditions.

Table 5 :
Input parameters used in ANSYS CFX.

Table 7 :
Variable values for Tunnel 2 at four different mesh sizes.

Table 12 :
Input parameters for the similar geometric components.

Table 13 :
Specifications for different regions of Tunnel 3.

Table 14 :
Comparison between Experimental and Numerical result.