Values of parameter: *y+*, skewness, and aspect ratio for the analyzed numerical grid.

## Abstract

This chapter deals with the processes by which a single-phase 3-D CFD model of hydrodynamics in a Sulejow dam reservoir was developed, verified, and tested. A simplified volume of fluid (VOF) model of flow was elaborated to determine the effect of wind on hydrodynamics in the lake. A hexahedral mesh with over 17 million elements and a k-ω SST turbulence model were defined for single-phase simulations in steady-state conditions. The model was verified on the basis of the extensive hydrodynamic measurements. Excellent agreement (average error of less than 10%) between computed and measured velocity profiles was found. The simulation results proved a strong effect of wind on hydrodynamics, especially on the development of the water circulation pattern in the lacustrine zone in the lake.

### Keywords

- CFD model
- hydrodynamics
- dam reservoir
- Sulejow reservoir

## 1. Introduction

In recent years, numerical models of flow in water bodies are widely used to provide an accurate description of flow velocity distribution. Although all the hydrodynamic processes going on in an aquatic system are difficult to describe, some authors have reviewed the main mechanisms, from large scale [1, 2, 3, 4, 5] to small scale and mixing processes [6, 7] or both [8, 9]. Three-dimensional numerical models of flow calculations in natural water bodies have been successfully applied for rivers with single complex flow features [10, 11, 12, 13, 14] although there is a little experience with modeling of flow in large dam reservoirs [15]. Studies of flow through curved or straight, open channels of simple cross-sectional shape (rectangular or trapezoidal) have been reported by Demuren [16] and Meselhe [17]. Both authors presented a 3-D numerical model, for the calculation of turbulent flow in meandering channels. The works provided a hydrodynamic basis for the study of the mechanisms for the formation of river meanders. The authors employed a finite volume numerical method to solve the full Reynolds-averaged Navier-Stokes (RANS) equations in conjunction with the standard k-ɛ model. The case exhibits some of the features encountered in real rivers, including longitudinal curvature and varying bed topography. Demuren’s results were in reasonable overall agreement with the mean velocity measurements of Almquist and Holley [18]. Meselhe focused on proposing a simplified approach for calculating the water surface elevation as a part of the overall solution procedure.

All previous studies with both straight and curved channels had adopted a rigid-lid assumption for modeling the free surface. The rigid-lid approximation is a commonly used simplification in the study of density-stratified fluids mainly in oceanography. Assumes, that the displacements of the surface are negligible compared with interface displacements. In a research conducted by Sinha [11], a 3-D model of flow through a 4-km stretch of the Columbia river, downstream of the Wanapum Dam, was developed and validated. Authors succeeded in modeling of flow in rapidly varying bed topography and the presence of islands. In the study three-pronged strategy—composed of the field measurements, the laboratory experiments, and the numerical model—was undertaken. The model solves the (RANS) equations closed with the standard k-ɛ turbulence model.

In today’s development level of hardware and software, increasingly complex, engineering tasks can be managed, and, thus, three-dimensional hydrodynamic modeling of dam reservoirs became a potential tool for determining more accurate picture of the dependencies, which prevail in the aquatic ecosystem. Whereas in most studies one- and two-dimensional description of flows may be sufficient, in some special cases, three-dimensional approach is needed to determining, for instance, flow patterns in bends or in the vicinity of hydraulic structures (dams and weirs). These issues are evidently three-dimensional, and spatial character directly affects pollutant transport processes.

Hydrodynamic processes determine the movement of suspended and dissolved matter; heat transfer; intensity of circulation inside the ecosystem; the speed of contaminating processes; and the self-purification of the reservoirs, ultimately provide conditions of the ecosystem function.

Kennedy et al. [15] applied a 3-D model to estimate the hydraulic residence time (HRT) for the Thomas Basin (part of the Wachusett Reservoir in central Massachusetts). The basin was modeled using the FLUENT software package with particles used to track travel time in steady-state conditions. A tetrahedral mesh was used with accurate description of basin bathymetry. The model solved the transient Reynolds equations for turbulent flow with standard k-ɛ closure. Modeling was performed to simulate flow pattern during a period when conditions were isothermal and windless. HRT was estimated to 3–4 days which is about half of the HRT that would be expected based on the theoretical mean residence time. The results of the calculations show that the presence of a primary flow path, large-scale eddies, and stagnation zones contributed to the faster travel time. Reductions in inflow rates produce increased residence times and significant changes in flow patterns. However, the authors did not provide any information about model verification.

In a research conducted by Khosronejad [19], a 3-D CFD model was applied to predict the flow hydrodynamics around power intakes within the Dez dam reservoir (Iran). The study was also devoted to a qualitative analysis of sediment transport at the area around the intakes. For incorporating the effects of turbulent flow, the k-ω model was implemented. The finite volume method was used to discretize the RANS equations. The 3-D single-phase model ran for a limited 320 m long reservoir section. The results show that the flow velocity has a maximum value of 1–2 m/s near power intakes and decreases with distance. In addition, the turbulent intensity increases in the area near intake entrance resulting in increasing bed shear stress near intakes. Khosronejad did not present a test of either coarser or finer mesh resolution; moreover, the model was not verified in the field measurements.

## 2. Study area

Sulejow reservoir is situated on the middle reach of the Pilica river, which is left-hand side tributary of Vistula river, central Poland (Figure 1). One of the biggest artificial reservoirs in Poland was built by impounding the Pilica river on 138.9 km with a dam in the years 1969–1973. The reservoir is a shallow water body (mean depth 3.3 m) covering a large area (22 km^{2}).

Sulejow reservoir is a ribbon-type reservoir, which can be divided into two morphological zones each influenced by different forcing agents. The first one (consisting of a riverine zone and a transition zone) is the narrow, shallower part of the reservoir, dominated by the river inflows. The second zone, the wide lacustrine part of the reservoir, is located near the dam. This zone is quite open and behaves like a lake, and the main driving force mechanism causes the movement of water masses wind. The main axis runs from southwest to northeast which is close to the direction of winds that ripple and mix the water. A result is the formation of places with stagnant water on the southern bank of the middle and lower part of the reservoir.

## 3. Generation of three-dimensional geometry

The first step in the numerical procedure was to prepare the 3-D geometry of the 17-km-long Sulejow reservoir. Gambit Program 2.2.30 (ANSYS, USA) was employed to generate geometry of the water body.

On the basis of the most accurate, available data of the reservoir, geometry was generated using segmentation technique with 36 cross-sectional profiles of the artificial lake. Methodology adopted to determine parameters of the sections was based on the field measurements by using an integrated system: digital depth sounder RESON PC-100 and GPS Trimble 5700. The cross-sectional spacing ranged from 400 to 600 m. Figure 2 shows an example and locations of the measured sections in the water body.

In this study, the “bottom-up” approach was implemented in order to obtain a three-dimensional geometry from the set of separated surfaces. Each measured cross section was defined by the information of the distance between two banks (connection length between probing points in the section was 5 m) and elevation (normal impoundment level in the Sulejow reservoir is 166.6 m amsl). Figure 3 illustrates an example of the cross section built in the Gambit Program.

The vertices which determine the reservoir bed were connected to form the edges. Subsequently, edges were combined into faces. The last stage consisted in stitching the faces in order to obtain the volumes. After completion of the segmentation procedure, rendering process was conducted, which facilitated generation of three-dimensional geometry of the reservoir, using the faces obtained from the segmentation. Figure 4 presents the complete geometry that consists of 36 volumes totally.

Terminal berm of the reservoir consists of earth dam and weir, with integrated hydroelectric power station (Figure 5). The length of the dam with a weir is 1200 m, the maximum height is 16 m, and the total volume is 567,000 m^{3}. The jazz is concrete, in the overflow riffle span, middle and left of the weir, drain pipes have been built. Weirs are closed with the oval valves. Cross section of earth dam and weir is shown in Figure 6a and b. Geometry and computational mesh of the dam are shown in Figure 7.

The final verification of the accuracy of the reservoir shape was feasible on the basis of satellite photographs, which confirm the correct approximation of the geometry. Figure 8 shows the satellite image of the Sulejow reservoir with marked cross sections, according to which the 3-D geometry was prepared. Due to the fact that sections were determined every ∼500 m, it could not be possible to precisely map the shape of the shoreline; however, this approximation does not affect the nature of the flow throughout the 17-km-long dam reservoir. The generated 3-D geometry of the artificial lake was discretized to perform the numerical solutions.

### 3.1. Computational mesh

Three-dimensional grid generation for a complex geometry is not a straightforward task, as the grid type has a significant influence upon the quality of the simulations [20]. The crucial issues controlling grid quality are the type of grid, that is, structured and unstructured, grid spacing, and skewness. To apply a finite volume method in CFD model of the reservoir, the 3-D computational domain must be subdivided into a large number of cells. However, the high ratio of the breadth to depth dimensions and irregular shape of the reservoir make this process difficult. Appropriate mesh resolution is linked to the hydrodynamic conditions, the flow features, and the discretization schemes.

In finite volume methods for numerical simulations, the hexahedral mesh (hex mesh) is preferred, as compared to the tetrahedral mesh (tet mesh), owing to the reduced error and smaller number of elements [21]. Typically, tet mesh is preferred for filling irregular spaces, since existing algorithms can semiautomatically subdivide most of the spaces [15]. Structured meshes are better suited to shallow reservoirs, while an unstructured mesh matches better to deeper (or with smaller aspect ratio) reservoirs. On the other hand, generating a hex mesh, with desirable qualities, often requires significant geometric decomposition, which makes the meshing process extremely difficult to perform and automate. As a result, it requires considerable user efforts and may take days or even weeks to develop the proper grid in the case of complex shapes.

To fill the shallow space representing the reservoir area, with minimally skewed, hexahedral cells, the typical cell dimension must be small compared to the depth of the water. Thus, the meshing process becomes computationally expensive due to the requirement of a large number of elements.

While developing the CFD model of flow hydrodynamics, preliminary simulations allowed us to select proper density of the numerical grid, at which a convergent and stable solution can be obtained. In order to generate a mesh for the Sulejow reservoir geometry, the capacity of Gambit 2.2.30 commercial software was used. The domain surface was discretized using structural (hexahedral) mesh with 16,787,820 active cells and 17,717,364 nodes, respectively. The total basin volume was 74,721,600 m^{3}, which reflects the real value. Four layers of elements across the wall thickness were added to confirm the numerical stability by specifying the height of the first layer, nearest to the wall (0.01), and the growth factor (0.1).

The process of generating structured grid was much more labor intensive than creating an unstructured mesh; however, elaborated model was more stable and converged quickly.

Exclusion from the model of some near-shore regions of the basin, with the shallow depth (<0.3 m) and inconsiderable slope, was necessary to avoid generation of cells with highly acute angles. Removal of strongly skewed elements from the computational domain would have minimal impact on the overall circulation patterns in the basin due to the high ratio of bottom area to the water volume, resulting with minimal flow and very low velocity in this region. This operation removed approximately 10% of the basin surface area but only a small fraction of the basin volume ∼1%.

A post-processing check on mesh quality, based on assessing the skewness of the generated cells, indicated that the mesh is of high quality and would not compromise solution stability. The use of a coarser grid (<10^{5} and 10^{6} elements) caused rapid solution divergence. However, increasing the mesh number to 18 × 10^{6} had no further influence on the results especially on the location of swirl flows. Figure 9 shows the hexahedral, structured mesh which has been generated for the Sulejow reservoir geometry.

The quality of the numerical grid was determined by the shape and the size of the computing field and the total number of elements used in the generated numerical grid and through the position of the first node relative to the plane of the wall. To assess the quality of the numerical grid elements, three parameters were used (Gambit User Guide):

Parameter

*y+*—determining the quality of the mesh in the boundary layer“Skewness”—determining the quality of the individual grid elements

“Aspect ratio”—defining the degree of deformation of the mesh elements

Ranges of parameters: y^{+}, skewness, and aspect ratio for the numerical grid generated in this work are given in Table 1.

Parameter | Values for numerical grids in this work | The limit values in the numerical simulations |
---|---|---|

y+ | ∼1.8 | ≤2 |

Skewness | ∼0.8 | 0–1 |

Aspect ratio | ∼1.7 | ≥1 |

## 4. Mathematical model

In order to fully present the hydrodynamics in the dam reservoir, volume of fluid (VOF) model should be developed. Owing to large size and complexity of the computational grid, as well as limited computer capacity, such an approach would be very difficult to apply for a big water bodies. In the literature, VOF technique was adopted only for the CFD modeling of reservoir in a downscale model (1:50) [22].

Sulejow reservoir due to the location, shape of the bowl, and uncovered, flat shores is particularly exposed to the effect of wind. Direction and energy of wind determine the waving movement and mixing of the water, so the impact of the factor, on flow distribution in the CFD simulations, should be taken into consideration.

In the model under wind conditions, the speed, at which the plate was moving on the water surface (which reflects the speed and direction of the wind in the area of the Sulejow reservoir), was determined. For this purpose, independent, 2-D, two-phase problem was resolved. Data concerning wind parameters were provided by the Institute of Meteorology and Water Management in Warsaw and come from the meteorological station (Sulejow-Kopalnia) located near the reservoir. Average speed and direction of wind were selected as 3 m/s and southeast, respectively, based on the daily data from year 2007.

### 4.1. Simplified VOF CFD model

The 2-D model of a straight channel with a length of 80 m and a width of 6 m (3 m layer of water, 3 m layer of air) was elaborated. Figure 10 shows the hexahedral, structured mesh which has been generated for the simplified geometry. The grid contains 159,242 elements (159,116 hexahedra, 126 wedges) and 321,872 nodes, respectively. A boundary layer was generated consisting of 10 rows.

Wind accelerates surface fluid particles by imparting momentum to the fluid, through surface stresses. In the analysis, the water flows through the air momentum. In the computational domain, the water phase has two outlets that allow for fluid reversing. Analysis of flow in 2-D model was intended to determine which boundary condition best describes the situation that prevails over the water surface by the effect of wind.

### 4.2. Boundary conditions

For two-dimensional CFD model, the following boundary conditions were imposed:.

Inlet boundary conditions applied for the analyzed domain were two tributaries (Pilica and Luciaza rivers) and outlet (dam) as given in Figure 11.

Simulated inflow boundaries were specified with mass flow rates, normal to the boundary. Velocities at each inlet were calculated from the inlet area measurements of stream flow for the Pilica and Luciaza rivers, made by the Regional Board and Water Management in Warsaw in 2007. The monthly values of mass flow rates in the Pilica and Luciaza rivers are presented in Table 2.

Month | Pilica river (mass flow rate m^{3}/s) | Luciaza river (mass flow rate m^{3}/s) |
---|---|---|

January | 19.86 | 2.80 |

February | 46.31 | 5.50 |

March | 22.33 | 4.92 |

April | 23.04 | 2.19 |

May | 17.36 | 1.39 |

June | 15.65 | 1.64 |

July | 16.45 | 1.91 |

August | 12.83 | 1.39 |

September | 13.21 | 1.59 |

October | 16.37 | 1.69 |

November | 25.50 | 2.21 |

December | 25.65 | 2.55 |

The definition of the inlet requires the values of the velocity vectors and turbulence properties. For the air inlet, the simulations were first conducted at a speed of 3 m/s. Velocity profile obtained at the outlet of the computational domain was loaded as an input file to receive the velocity profile at the inlet. This approach allows to obtain a fully developed velocity profile for a small domain.

Implementation of volume of fluid model requires an additional boundary condition to be specified, namely, the turbulence intensity at the inlet and turbulent viscosity ratio. Introduction of disturbance into the flow reflects the real features of the flow pattern.

The turbulence intensity, *I* [%], is defined as the ratio of the root-mean-square of the velocity fluctuations *u`* to the mean flow velocity *Uavg*. The turbulence intensity at the core of a fully developed duct flow can be estimated from Eq. (1):

Kennedy [15] carried out sensitivity analysis on the turbulence intensity in the simple, channel geometry and found that value 4% was able to match the conditions well. The author concluded that this parameter had a noticeable effect on the solution. Following the above findings, turbulence intensity level of 4% was specified at the inlet.

Another parameter, the turbulent viscosity ratio (the ratio of turbulent to laminar (molecular) viscosity), was used as given in Eq. (2). The default value of 10% was applied in the simulations:

The air outlet of pressure type was defined in the model. Pressure outlet boundary conditions require specification of a static pressure at the outflow. Convergence difficulties were minimized by specified values for the backflow quantities (backflow turbulence intensity and viscosity ratio). A no-slip boundary condition was applied on the wall.

### 4.3. Solution methods

Table 3 summarizes the solution conditions and methods used in the modeling process.

Model volume of fluid | |

Space | Two-dimensional |

Time | Steady |

Discretization method | |

Gradient | Least squares cell based |

Pressure | Body force weighted |

Pressure–velocity coupling scheme | Coupled |

Momentum | Second-order upwind |

Volume fraction | Compressive |

Turbulence energy kinetic | Second-order upwind |

Specific dissipation rate | Second-order upwind |

### 4.4. Results of VOF model

Based on the analysis of two-dimensional two-phase model, an appropriate input data for the wall boundary condition, corresponding to real wind conditions at the surface, could be determined. Figure 12 depicts the velocity profile of water at the outlet of the computational domain. As a result, the replacement velocity of wind, which corresponds to the real, average value in the area of the Sulejow reservoir, was obtained. The mean value of 0.147 m/s was used in the simulations of flow hydrodynamics in the modeled reservoir.

### 4.5. Results of CFD calculations at wind conditions

Due to the fact that the Sulejow reservoir is a shallow, polymictic lake, the wind will be important for the distribution and mixing of the water masses of two tributaries. The results shown in Figure 13 confirm this assumption.

The findings suggest that when steady flow pattern develops in the basin, large regions of recirculation are formed below the outlet of the reservoir. Figure 14 shows contours of velocity field in the Sulejow reservoir in (a) March (b) July, and (c) December.

## 5. Model verification

Numerical simulations have many advantages such as providing results in the entire domain and the ability to make changes to the geometry, boundary, or initial conditions, but numerical models always require validation of the simulations with reliable and appropriate experimental data [22]. The following section describes the measurements for obtaining data necessary for validation of the numerical simulations.

### 5.1. Acoustic measurements

Field measurements of velocity and turbulent quantities have been made much more accessible with the development of acoustic measurement methods. Acoustic Doppler Current Profiler (ADCP), highly efficient and reliable instrument for flow measurements in riverine and open-channel environments, has been used for the first time to determine the velocity profiles in the dam reservoir in Poland.

The measurements are repeated continuously during the movement of the boat. As a result, for a single passage along the cross section, a few hundred to several thousands of partial flow are obtained, which are summed during the measurement process. The measurements are repeated several times. The final result is calculated as the average value of at least four correct runs.

Teledyne RD Instruments’ StreamPro Acoustic Doppler Current Profiler (ADCP) was used to validate the hydrodynamic CFD model. The device was mounted on a boat that moves across a transect of the reservoir channel. This technique can be adopted while complying the assumptions: (1) the water surface is not wavy and (2) velocity of the water is less than 2 m/s.

In order to correctly determine the places selected for the verification, the GPS Garmin eTrex 10 has been used. Flow velocity measurements were made in June 2013 and included four cross sections in the Sulejow reservoir and measurements of the flow rate at the inlet to the reservoir of the Pilica and Luciaza rivers.

The places, which have been selected for the verification of the hydrodynamic model, were arranged along the longitudinal axis of the Sulejow reservoir. The selection of these areas was dictated by expecting a different character of flow in indicated areas (Figure 15):

Barkowice Mokre (1)—the place closest to the backwaters of the reservoir, located in the riverine zone, characterized by small depths (<3 m) and the highest flow rates.

Zarzecin (2)—located in the upper, narrow part of the reservoir, with higher flow velocities, resulting in a half-river character. The depth at this point was about 4 m.

Bronislawow (3)—situated in the central part of the reservoir, near the former water intake for Lodz City. The place is characterized by a low flow, due to greater width (1500 m) and depth of the basin (approximately 6 m).

Tresta (4)—located in the lower part of the artificial lake, where the depth is about 7–8 m. The place is closest to the dam of the reservoir, characterized by, in addition to great depths, the largest cross-sectional width of approximately 2000 m.

### 5.2. Results of model validation

A reasonable agreement between the flow pattern predicted by the model and those deduced from the field data was found (Figure 16). To quantify the comparison, relative error was calculated by taking the difference between numerical and measured values and then dividing the results by the measured values. The measured velocity profiles were provided with a relative accuracy within the range 1–10%.

The possible reasons for the discrepancies are (1) inaccuracies in location of measuring points, (2) point velocity measurement errors, (3) errors in modeling the flow, and (4) errors in modeling the geometry.

The first category is related to the field velocity measurements taken from a boat. Considering the fact that a boat cannot maintain an absolute fixed position due to the waving and wind, errors are introduced in velocity measurements. A deviation of ±20 cm from the fixed position can cause large errors if there is a steep velocity change in the plane of measurements. The magnitude of this error could not be estimated accurately; however, rough estimate of the nearby velocities within a distance of ±20 cm at the measuring point resulted in an error in the range 3–5%.

The second category consists of errors related to the instrument, its volume resolution, the range of operation, and the sampling time. The ADCP device could measure instantaneous 3-D velocity vectors with 1% accuracy. The vertical resolution of the instrument was 0.05 m, which is less than the vertical mesh spacing (∆y) used in the numerical model, that is, 0.08 m. In other words, the instrument resolution error can be ignored.

The third category is related to the numerical methods (discretization and iteration errors), the boundary conditions, and the closure models. For a carefully modeled problem that has well-posed boundary conditions, these errors are relatively low in comparison with other errors.

The fourth category is how the model geometry was built. The modeled geometry was an approximation of the reservoir topography as it was based on measurements of discrete cross sections. The regions between the cross sections were interpolated and may not represent the right topography of the artificial lake. A rapid variation in the topography significantly affects the flow velocity distributions. The spacing used in the present study was selected with special attention to the section properties of the reservoir. However, they may not have captured important changes of the bed.

Based on the above discussion, a total error between computed and measured velocities of about 10% is a reasonable assumption. Proper agreement of theoretical and experimental results shows the correctness of the actions, developed in the frame of this work.

## 6. Conclusions

The objective of this study was to develop and validate a three-dimensional numerical model for simulating flow through the long dam reservoir of a complex bathymetry (17 km length). As a result of the study, a three-dimensional one-phase CFD model of flow hydrodynamics in the large water body on the example of the Sulejow reservoir was developed with an accurate depiction of basin bathymetry and verified on the basis of field measurements.

The results of three-dimensional one-phase CFD model indicate that the flow field in the Sulejow reservoir is transient in nature, with visible swirl flows in the lower part of the lake.

The results of simulations confirm the pronounced effect of wind on the water flow in the reservoir and the accumulation of phytoplankton cells in the epilimnion layer of the lacustrine part of the Sulejow reservoir. Methodology developed in the frame of this work can be applied to all types of storage reservoir configurations, characteristics, and hydrodynamic conditions. Results of the simulation are complementary to the direct measurements of the surface water quality. A well-defined and constructed model can be used while developing a strategy for water environment quality control and can be used as an auxiliary tool for the monitoring and prediction of surface water quality and decision-making in the field of planning.