Abstract

Simulation-driven design with computational fluid dynamics has been used to evaluate the flow downstream of a hydropower plant with regards to upstream migrating fish. Field measurements with an Acoustic Doppler Current Profiler were performed, and the measurements were used to validate the simulations. The measurements indicate a more unstable flow than the simulations, and the tailrace jet from the turbines is stronger in the simulations. A fishway entrance was included in the simulations, and the subsequent attraction water was evaluated for two positions and two angles of the entrance at different turbine discharges. Results show that both positions are viable and that a position where the flow from the fishway does not have to compete with the flow from the power plant will generate superior attraction water. Simulations were also performed for further downstream where the flow from the turbines meets the old river bed which is the current fish passage for upstream migrating fish. A modification of the old river bed was made in the model as one scenario to generate better attraction water. This considerably increases the attraction water although it cannot compete with the flow from the tailrace tunnel.

1. Introduction

Computational fluid dynamics (CFD) is used to simulate the flow within a tailrace channel of a hydropower plant with the purpose to scrutinize alternative positions of an entrance to a fishway. The simulations are carried out on full scale implying a length of the virtual model of 320 m, a typical width of 75 m, a typical depth of 10 m, and a maximum inlet flow rate of 1000 m3/s. A numerical challenge with the large scale is to fulfill conditions of a sufficiently resolved flow structure at locations with high gradients (e.g., at boundaries) and a good mesh overall with a decent usage of computational resources (Marjavaara and Lundström [1]). Another challenge is the validation of the simulations which is here done by measurement with an Acoustic Doppler Current Profiler (ADCP).

Studies of tagged Atlantic salmon and sea trout in the unregulated river Vindelälven in northern Sweden during 1995–2005 have shown that only a third of the upstream migrating fish find their way to their natural spawning grounds (Lundqvist et al. [2]). The main reason for this is the Stornorrfors power plant located downstream the confluence between the rivers Vindelälven and Umeälven, the latter being a regulated river. A major issue at the power plant is that the fish are attracted into the tailrace channel from the turbines rather than migrating up through the old river bed that offers a fishway around the turbines (Rivinoja et al. [3]). The flow rate from the turbines is typically 20 times larger than the flow rate from the old river bed and its entrance into the confluence is very wide. Hence, fluid flow conditions for the old river bed to attract fish are limited. The fact that migrating fish are attracted to the tailrace of the turbines instead of the weaker current from the fishway is a common problem (Arnekleiv and Kraabøl [4], Webb [5]). The difficulties of upstream migrating fish coming across in regulated rivers in northern Sweden have been documented by, for example, Rivinoja [6], Lindmark [7], and Lindmark and Gustavsson [8].

There are two major measures that are being considered for improving the upstream migration of fish at the Stornorrfors power plant. One is to construct a new fishway in the form of a fish ladder from the tailrace channel since a majority of the fish reside there for a long period of time during the migration season. The other alternative is to create better attraction water from the old river bed into the confluence area. The alternatives are here modeled with CFD, and the attraction water created using given configurations is examined. The interest in numerical simulations of flows in rivers is increasing, and, due to the rapid development of user-friendly efficient codes and computer power in recent years, more advanced models than before can be applied in areas such as fish migration, habitat modeling, sedimentation transport, erosion, and dam safety. Olsen and Stokseth [9] created a model of the Sokna River in Norway where they applied a 𝑘-𝜀 turbulence model and a porosity-based model for large roughness elements in the river bed showing good resemblance with observed data. The SSIIM model suggested by Olsen has been validated against LDA measurements in a meandering channel on a lab scale (Wilson et al. [10]) where the model showed the ability to capture secondary currents. CFD has also been applied to the River Cole, Birmingham, UK (Clifford et al. [11]), and the River Thame, Birmingham, UK (Booker [12]), where the potential for use in habitat modelling was discussed. A numerical model of a 4 km stretch of the Columbia River downstream of the Wanapum Dam has been performed and calibrated against measured data highlighting the importance of bed roughness for accurate flow predictions (Sinha et al. [13]). Dargahi [14] used the commercial code Fluent to model fluid flow and sediment transport in the River Klarälven, Sweden, and validated the results with ADCP measurements. The design of a submerged flow guiding device to increase the survivability of downstream migrating fish has been performed in the commercial code CFX-10 (Lundström et al. [15]). Simulations of flow in an ice-covered channel with the 𝑘-𝜔 turbulence model with different roughness values for river bed and ice-cover resulted in a 16% increase in mean flow depth of the channel (Yoon et al. [16]). The effects of submerged weirs in natural channels to improve the navigation conditions for barges have been investigated numerically (Jia et al. [17, 18]).

Rakowski et al. [19] used field-measured data with ADCP to validate their CFD simulations of a 2.7 km reach starting downstream of the Bonneville powerhouse and spillway with total river flows between 3275 m3/s and 11328 m3/s. The velocities were measured and averaged over a 10 minutes period to get adequate representation of the mean velocity. When comparing the CFD simulations (steady state, 𝑘-𝜀 turbulence model) to ADCP data, the modeled velocity was slightly lower than the measured, but within the standard deviation of the field velocity. Viscardi et al. [20] also used ADCP measurements to validate CFD simulations in a 3 km stretch of the Paraná de las Palmas River with flow rates ranging from 2200 to 5000 m3/s (steady state, 𝑘-𝜀 turbulence model, rigid lid, bed roughness Manning 𝑛=0.025). In their case, the velocities were averaged over 2 seconds in each vertical sample in order to minimize the effect of the tidal change and the velocities correspond reasonable accurate. To summarize, two-equation turbulence models are in most cases used to simulate the flow in rivers and no one is considering how attraction water from a fishway competes with the flow in the river.

2. Geometry

The actual geometry in the present study consists of four parts, the tunnel from the turbines, the tailrace channel, the old river bed, and the confluence area, see Figure 1 where the tailrace channel and the old river bed are defined. The confluence area is located where the water from the old river bed and the tailrace channel meet, while the tunnel from the turbines is located upstream the tailrace channel. CFD calculations are performed on all parts except the old river bed, while velocity measurements are only reported for a couple of transects within the tailrace channel.

3. Experimental

To measure topology and water velocity downstream Stornorrfors power plant, an ADCP was used. The ADCP has four transducers directed into the water. The transducers send out sound waves that reflect on small particles traveling with the water, and the transducers detect the Doppler frequency of the reflected sound waves. These frequencies are proportional to the velocity of the water (the particles). ADCP is a relatively fast way of measuring velocities in field and to calculate river discharge. The ADCP used in this case is a RiverBoat RioGrande, and the data processing was performed with the software Winriver II, both from RD Instruments.

The bathymetry in the area was measured using two setups. The ADCP was dragged besides a motorboat with a pole and rope, which enabled measurements close to the shoreline. By combining the bottom-tracking feature of the ADCP with GPS data, a point cloud consisting of ADCP provided depths at specific satellite coordinates was obtained. The ADCP however fails to find the bottom of the deepest area in the tailrace channel; hence, a SIMRAD EY60, GPT 200 kHz, split beam echo sounder with the transducer mounted vertically on the boat was used near the tailrace tunnel outlet. There is a small shallow part of the tailrace that is located above and behind the tunnel outlet. The GPS reception was very low this far into the spillway due to the surrounding terrain which caused large uncertainties in the acquired coordinates. The length of this region is approximately 50 m, and observations suggest that it consists of a slow circulating flow and the assumption was made that this part of the spillway does not have any significant effect on the flow in the remaining channel. It was thus omitted from the numerical model, and the entrance to this innermost part was excluded from the geometry. The points of measurements are shown in Figure 1(a).

For the velocity measurements, a steel wire was stretched across the tailrace channel and the ADCP was tethered to it. A manual winch enabled the ADCP to travel across the channel and capture the velocities in the entire cross-section. The transect T1 in Figure 1(a) was measured on several occasions at different turbine discharges and a minimum of four times at each flow. Three vertical profiles in the T1 transect and three in the T2 transect were measured during a minimum of 600 s. The profiles were collected when the flow rate through the power plant was ~500 m3/s (according to the discharge calculation in WinRiver and data from the hydropower company). Profiles were measured with a time difference of 0.95 s between ensembles. During measurements, the distance to the shore was measured with a laser distance meter. The total width of the T1 section was measured to 40 m, and the profiles were located at 16, 23, and 32 m from the south shore. Transect T2 had a measured width of 85 m, and the verticals were located 30, 44, and 59 m from the north shore.

The accuracy of the ADCP depends on many factors, such as side-lobe interference, ringing, and ADCP-flow interaction that exclude the ADCP from doing any measurements near the water surface or close to the bottom of the river (Simpson and Oltmann [21]). Nystrom et al. [22] compared ADCP accuracy with an Acoustic Doppler Velocimeter (ADV) in a lab flume within which the turbulence intensity of the flow was 0.1. The ADCP measurements were carried out during 15 min, and the error was less than 3% in the areas away from the boundaries not being affected by ringing, side lobe interference, and flow disturbance.

4. Numerical Setup

The point cloud collected with ADCP and SIMRAD seen in Figure 1 was converted to a bottom surface in the software Imageware 13. The surface was imported to Ansys Icem Cfd 11 where a solid model was created. The formed geometry was divided in two parts, the tailrace channel and the confluence area between the channel and the old river bed. The simulation volumes were discretized as tetrahedral elements in the CFD model. Local refinements of the grid were carried out in areas of simulated attraction water to increase the resolution in the most interesting parts of the flow. All simulations were carried out with the commercial software CFX11 from Ansys Inc. A mesh sensitivity study of the tailrace channel was performed with different numerical grids, ranging from 239 k to 7389 k nodes. The velocity in the east direction was evaluated at T1 and T2 for the different grids and the conclusion was that the coarsest mesh did not capture the flow field with sufficient accuracy. A mesh with 526 k nodes however produced a very similar flow field to that of the 7389 k mesh with significantly lower computational cost. Hence, the final grids for the tailrace channel consisted of ~500 k nodes and the confluence area of ~600 k nodes.

In reality, the water from the power plant goes through an approximately 4 km long tunnel before entering the tailrace channel. To create a realistic inlet boundary condition for the simulations, this tunnel was modeled separately and the velocity profile at the end of the tunnel was used at the inlet of the tailrace channel simulations. The tunnel was given a sufficient length to give a fully developed velocity profile, and the tunnel walls were given a wall roughness of a typical excavated rock.

The high-resolution advection scheme was used for solving the equations of fluid motion and turbulence closure. The high-resolution scheme uses a close to second-order solution in areas with low variable gradients, and, in areas where the gradients change sharply, it will be close to a first-order solution (ANSYS [23]). The incompressible Reynolds-Averaged Navier Stokes equation and the continuity equation are expressed as𝜕𝑈𝑖𝜕𝑡+𝑈𝑗𝜕𝑈𝑖𝜕𝑥𝑗1=𝜌𝜕𝑃𝜕𝑥𝑖+𝜈2𝑈𝑖𝜕𝜕𝑥𝑗𝑢𝑗𝑢𝑖,𝜕𝑈𝑖𝜕𝑥𝑖=0,(1) where 𝑈 is the mean part of the velocity component, 𝑃 is the pressure, 𝜈 is the viscosity of the fluid, 𝜌 is the fluid density, and 𝑢𝑗𝑢𝑖 are the Reynolds stresses. All simulations were run with the 𝑘-𝜀 turbulence model with scalable wall functions. In the 𝑘-𝜀 turbulence model, Reynolds stresses are linearly related to the strain:𝑢𝑗𝑢𝑖=2𝜈𝑇𝑆𝑖𝑗23𝑘𝛿𝑖𝑗,(2) where 𝑘 is the turbulent kinetic energy, 𝜈𝑇 is the eddy viscosity, and 𝑆𝑖𝑗 is the mean strain tensor defined as𝑆𝑖𝑗=12𝜕𝑈𝑗𝜕𝑥𝑖𝜕𝑈𝑖𝜕𝑥𝑗.(3)

The eddy viscosity is modeled as𝜈𝑇=𝐶𝜇𝑘2𝜀,(4) where 𝐶𝜇 is a model constant and 𝜀 is the turbulent dissipation rate. For more information, see Launder and Spalding [24]. The RMS residual target for all simulations was set to 10−6. This convergence target could not be achieved with a steady-state solver due to initial fluctuations of the flow. Simulations were instead run on a transient solver until it approached a steady solution and the final values from the simulation were used. A physical time step of 0.5–2 s was selected depending on grid size, and the solution was considered steady when the velocity in 18 monitored points throughout the domain had been virtually constant for at least 1000 time steps.

The water surface was modeled as a rigid lid with zero friction. This approximation is viable when the surface level variation is smaller than 10% of the total channel depth (Rodriguez et al. [25]). The outlets are given a pressure type boundary condition. The bottom surface of the numerical model was defined as a rough wall. A scalable wall function that uses an extension of the method suggested in [24] was selected for near wall modeling. The dimensionless velocity 𝑢+ in the logarithmic layer close to the rough wall is typically written as𝑢+=𝑈𝑡𝑢𝜏=1𝜅ln𝑦++𝐵Δ𝐵,(5) where𝑦+=𝑢𝜏Δ𝑦𝜈,𝑢𝜏=𝜏𝜔𝜌,(6) and where 𝑢𝜏 is the friction velocity, 𝑈𝑡 is the velocity tangent to the wall at a distance Δ𝑦 from the wall, 𝑦+ is the nondimensional wall unit, 𝜏𝜔 is the wall shear stress, 𝜅 is the von Karman constant, 𝐵 is a constant, and Δ𝐵 is the so-called roughness characterization function. Since the roughness of the channel is not well documented, a global representation of roughness was selected. The wall roughness can be described as an equivalent sand-grain roughness, 𝑘𝑠 (ANSYS [23]). With this formulation, the roughness characterization function can be described as (White [26])1Δ𝐵=𝜅ln1+0.3𝑘+𝑠,(7) where the dimensionless roughness height 𝑘+𝑠 is defined as𝑘+𝑠=𝑢𝜏𝑘𝑠𝜈.(8)

To obtain a realistic value for the equivalent sand-grain roughness of the channel, an empirical Gauckler-Manning coefficient 𝑛 that describes the channel is selected (Arcement and Schneider [27]). The advantage of using the Manning’s 𝑛 instead of other coefficients is that 𝑛 is nearly constant regardless of flow depth, Reynolds number (Re=4𝑈𝑅/𝜈), or 𝑘𝑠/𝑅 for fully developed turbulent flow over a rough surface (Yen [28]). The selected 𝑛 is used to calculate a Darcy friction factor 𝑓 given by𝑓=8𝑔𝑛2𝑅1/3.(9)

The friction factor obtained is then used to find 𝑘𝑠 from the Colebrook-White formula (Colebrook [29])1𝑓=2log10𝑘𝑠3.714𝑅+2.51Re𝑓,(10) and the derived 𝑘𝑠 is finally used for input into the simulations. In the present case, the second term in (10) can be neglected due to the high Reynolds number of the flow.

The wall function approach is common in river simulations since the scales of roughness are very costly to model physically in problems of such large scales. The limitations of this method are discussed by Patel [30]. Since the roughness is only an approximation and is difficult to measure in reality, a parameter study was here performed in the numerical model with a flow rate of 350 m3/s from the turbines.

With no surface roughness, the jet leaving the tunnel barely leaves the bottom of the channel which does not seem likely with regards to the characteristics of free surface channel flow, see Figure 2. With a 𝑘𝑠 value of 0.3 m (Manning 𝑛0.033) which can be considered typical for a rock excavated channel such as the tailrace channel, the flow characteristics change considerably. The jet emerging from the tunnel now moves towards the free surface of the channel. Increasing 𝑘𝑠 to 0.5 m (Manning 𝑛0.037) did not affect the solution in any major way, and all following simulations on the tailrace channel, were run with 𝑘𝑠=0.3 m.

Two ways of improving the upstream fish migration around the power plant were studied: a new fishway in the tailrace channel and higher attraction to the old river bed. In the tailrace channel, two positions and two angles of a new fishway entrance were evaluated. The positions were selected from previous observations of fish during the migration season. The dimensions of the entrance were 2×2.7 m2, and the flow rate used was 10 m3/s. The two inlet angles of the fishway entrance were perpendicular and 45° to the main flow.

The modification of the confluence area to improve the attraction to the old river bed was realized by adding a wall at a distance from the river bank directing nearly all the flow in the old river bed to a narrow open channel between the wall and the shoreline. In such opening the water may be accelerated with a ramp, for instance, as suggested in Lindmark and Gustavsson [8] and Green et al. [31]. The flow in the old river bed was set to 20 m3/s and that from the tailrace tunnel to 500, 750, and 1000 m3/s, representing a low flow, a normal flow, and a flow close to the maximum flow, respectively.

5. Results and Discussion

The characteristics of the flow will be described followed by a comparison to experimental data, and finally the results from simulation of attraction water will be presented.

5.1. Characterization of the Flow in the Simulations

The flow exiting the tunnel takes the form of a jet that gradually develops into an open channel flow profile. Approximately after 2/3 of the channel length, the jet maximum velocity is at the surface of the water as can be seen in Figure 3.

The jet however influences the surface orientation flow much earlier in the channel as revealed by plots of vorticity and turbulence intensity in a plane at 1 m depth, see Figures 4(a) and 5(a). High-vorticity areas are found at the bottom of the channel near the edges of the channel and at one large area of recirculation after the expansion, near the north shore, but there is also a noticeable rotation of the flow close to the surface near the inlet of the tailrace channel, see Figures 4(a) and 4(b). This is also reflected by relatively high turbulence intensity in this area, see Figures 5(a) and 5(b).

5.2. Comparison to Experiments

The results from ADCP measurements in the tailrace channel yield an unstable behavior of the flow, see Figure 6 showing a 12×12 m2 section in the middle of the T1 transect at five different times where the raw data from the ADCP has been averaged to 1×1 m2 cells. The measurements were taken in succession, and the velocities have been normalized with transect average velocity to account for minor differences in total flow from the tunnel. The jet stemming from the tunnel is apparent in all transects but not as well defined as in the simulations, compare Figures 3 and 6. To examine the time dependence of the flow, the ADCP was kept in the same position and the velocity was measured during a longer period of time. Three vertical profiles at 15, 22, 31 m from the south shore were measured. The standard deviation from the mean distance was 0.01-0.02 m. The results from the measurements show a highly fluctuating flow. Initial frequency analysis does not indicate any periodicity; however, it cannot be excluded that fluctuations are influenced by large-scale structures of the flow, originating from upstream instabilities. How the RMS velocity (east) stabilizes with time is shown for the profile at 22 m in Figure 7. From the results, it is concluded that, to measure representative mean velocities, each profile must be measured during at least 600 s. The measurements over a complete transect presented in Figure 6 took about 120 s which means that these measurements by no means represent the mean velocity in that transect which explains the different velocity patterns.

To validate the simulations, time-averaged velocities of fixed-point measurements at both T1 and T2 are derived. The agreement between simulation and experiment at T1 is rather poor, see Figure 8(a) where normalized velocity profiles are compared. The velocity is normalized with the bulk velocity 𝑈bulk=𝑄/𝐴𝑇1, where 𝑄 is the flow rate and 𝐴𝑇1 is the area of the T1 transect being 516 m2 as derived from the virtual model. The jet that exits the tunnel appears closer to the water surface in the measurements than in the simulations, and it is much more diffuse in the measurements. This is most apparent for the measurements at 31 m at T1 where measurements indicate a plug flow while the simulations yield a sinus-shaped profile. Hence, there is a discrepancy at the surface and at the bottom and the jet penetrates the surface much earlier in reality as compared to the simulations. For T2, the agreement between simulations and experiments is better especially close to the free surface, see Figure 8(b). The maximum velocity of the flow in the middle of the channel is lower in the experiments than in reality, while it is actually higher towards the shores indicating a more diffusive flow also in this transect. One reason for the differences, especially apparent in T1, might be the inlet boundary condition in the simulations, which is described as a stationary velocity profile where in reality effects of the turbines, larger discrete wall roughness elements or sudden changes in discharge may come into play. Other contributing factors may be difference between model geometry and real geometry as to surface roughness, for instance, and oversimplified modeling of turbulence or that the rigid lid assumption creates unphysical behavior when the jet from the tailrace tunnel approaches the water surface of the tailrace. It is also likely that the flow field is smeared out by the method to measure the velocity field. The discrepancy between simulations and measurements is a subject for future research as to turbulence intensity, for instance. When later on discussing the results from the simulations with the fishway entrances, it should be remembered that the jet is more diffuse and surface orientated in reality as compared to the simulations.

5.3. Simulation of Attraction Water

For position 1 in the tail race channel, the perpendicular entrance gives a noticeable jet that stretches to the center of the tail race channel while the angled inlet gives a jet that aligns with the flow from the tailrace tunnel and reaches further downstream, see Figure 9. Even better attraction water is created at the second position as shown in Figure 10. Since the small jet from the fishway does not collide with the large jet from the tailrace tunnel, the generated attraction water stretches further out in the channel, see Figure 11. Noticeable attraction water was created even at the highest flow (1000 m/s) from the turbines, see Figure 12. The relatively high-vorticity levels and turbulence intensities are thus too weak to influence the attraction water to any larger extent. Hence, position 2 is, as to generation of attraction water, a better choice than position 1. This conclusion is strengthening by the fact that, in reality, the jet from the turbines is more surface orientated which probably will make the attraction water created at position 1 less prominent than obtained in the simulations and from the turbines stressing the fact that the fishway should be placed as long into the tail race channel as possible for optimum generation of attraction water.

When scrutinizing possible improvement of the attraction water from the old river bed, the simulations yield a rather different result. The attraction water cannot compete with the flow from the tail race channel except in an area quite close to the shore, see Figure 13 where the simulated attraction water competes with two flow rates from the turbines (500 and 750 m3/s). As compared with the current situation, this modification of the confluence would still provide considerable improved attraction water along the north side (the right-hand side in the simulated results in Figure 13 and see Figure 1 for cardinal directions). This should improve the probability that fish migrating upstream on the north side of the river or fish exiting the tailrace tunnel on the north side find the fish passage in the old river bed.

6. Conclusion

The measurements show that the flow is considerably more unstable in reality as compared to the simulations. The flow fields in the simulations are therefore less diffuse as to time-averaged quantities, and the tailrace jet from the tunnel outlet is stronger but less surface orientated in the simulations as compared to reality. Keeping this in mind, a number of additional conclusions can be made from the work here presented. A fishway in the tailrace channel can generate noticeable attraction water for all relevant flows from the turbines. Of the cases studied, the simulations show that position two gives considerably stronger attraction water as compared to position one. It is likely that this difference is enlarged by the diffusivity of the tailrace jet existing in reality. By a concentration of the flow from the old river bed, noticeable attraction water can be created at the confluence area. In this case, the attraction water only stretches a short distance into the tailrace flow since this flow is completely surface orientated in the confluence area. However, if the fish migrate along the north shore, they will sense the attraction water.

Acknowledgments

This work was financed by Vattenfall AB and Umeå Kommun. The authors also acknowledge Vindeln Utveckling that always kept their administration in phase with the real EU time and the National Board of Fisheries for their suggestions on the design of this study.