Selected applications of an open-source three-dimensional computational fluid dynamic code in hydraulic engineering

Exponential growth of computational power in the past decade on one hand and the success of numerical models in studying hydro-environmental flows on the other hand lead to an increasing usage of this technique in hydraulic engineering problems. In particular, three-dimensional numerical models are becoming the industry standard, especially in problems where capturing three-dimensional flow characteristics are important. Additionally, with availability of open-source computational fluid dynamic (CFD) codes, this type of modelling is more accessible. Among these codes, OpenFOAM (Open Source Field Operation and Manipulation) is a state-of-the-art CFD code with an extensive range of features for engineering and scientific applications. In this paper, selected studies are presented which the authors carried out using OpenFOAM. These studies include a wide range of applications from implementation of the new algorithms to modelling a practical application in hydraulic engineering.

Schlüsselwörter OpenFOAM · Turbulente Strömung · Large-Eddy Simulation · Detached-Eddy Simulation · Dünen · Sekundärströmung · Flusskraftwerk 1 Introduction Exponential growth of computational power in the past decade and the success of numerical models in studying hydro-environmental flows lead to an increasing usage of this technique in hydraulic engineering problems. In particular, three-dimensional numerical models are becoming the industry standard where the three-dimensional flow characteristics become important. Additionally, with availability of opensource computational fluid dynamic (CFD) codes, this type of modelling is more accessible. Among these codes, OpenFOAM (Open Source Field Operation and Manipulation) is a stateof-the-art CFD code with an extensive range of features for engineering and scientific applications. In principle, this code is a collection of C++ libraries written with object-orientated approach, hence, it facilitates implementation of new algorithms or modification of the existing codes. Furthermore, the code employs Finite Volume Method (FVM) on unstructured grids to solve the flow equations. In this paper, selected studies that are carried out by the authors are briefly presented. All simulations are carried out using OpenFOAM. The post processing and visualization are carried out using Paraview and Python.
In the first part of this paper an example of implementation of a new solver in OpenFOAM is presented. The aim of the solver is to speed up the computation especially for Large Eddy Simulations (LES). The code is then validated using two classical cases of turbulent flows: turbulent channel flow and flow over a backward-facing step. In the second part, flow over a twodimensional dune is modelled using a novel approach known as Improved Delayed Detached Eddy Simulation (IDDES) technique (Shur et al. 2008). The results are then compared to values from a reference LES study and corresponding experiment.
Furthermore, a brief description of the hydrodynamical processes is presented. Another study, in which the flow in a trapezoidal channel is modelled employing the same technique. The results are compared with theoretical and empirical values. Additionally, the secondary currents in such channels are presented and briefly discussed. Finally, in the last part, an approach flow of a Run-of-River plant (ROR) is modelled using Reynolds-averaged Navier-Stokes (RANS) approach. A series of parameter studies are carried out to identify suitable parameters for such a modelling approach. Additionally, the numerical results are compared to the experimental values. therefore these types of studies are computationally expensive. Furthermore, simulations of high Reynolds number flows in complex geometries require a large number of CPU cores and they typically take days or even months of computation. Therefore, if an algorithm can be employed to reduce the required computational time, it is desirable as it can lower the cost of LES significantly. In OpenFOAM, two solvers are available for transient, incompressible turbulent flow problems. Pressure-Implicit with Splitting of Operators (PISO) and PIMPLE. The latter, is a combination of Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) and PISO. Furthermore, these algorithms are using first-and secondorder implicit time integration schemes for solution of the flow equations.
In LES, on one hand the Courant-Friedrichs-Lewy (CFL) number must be kept well below one. On the other hand, a numerical algorithm with low dissipative properties is desirable. It is shown that a solver based on high-order explicit Runge-Kutta (RK) family of methods in combination with the fractional step can reduce the computational time as well as decrease numerical errors in comparison to PISO algorithm (see e.g. Vuorinen et al. 2014). This algorithm is implemented in OpenFOAM by modifying PISO algorithm. In this study, explicit third-and fourth-order Runge-Kutta fractional step algorithms (here referred to as RK3 and RK4 respectively) are implemented. These algorithms are then tested via several test cases. In the following, two of these cases are presented and briefly discussed. Full description of the method and their implementation are described in Lora (2019).

Turbulent channel flow
In wall-bounded turbulent flows, fully developed turbulent channel flow is a classical case. The idea behind this numerical test is a flow between two infinitely long plates driven by a constant pressure gradient. This case is extensively studied numerically by LES and Direct Numerical Simulation (DNS) methods for different flow conditions. Here, DNS results of the turbulent chan-nel flow, which is investigated by Moser et al. (1999), are compared to the results obtained with the aforementioned solver (RK3). This is done with two objectives in mind: first, to validate the solver and second-to compute the speed up that might be achieved by utilizing this solver.
The flow simulations are carried out at friction Reynolds numbers Re τ = u τ δ/ν = 180,395 and 590, where u τ is the friction velocity, δ = 1 m is the half distance between the walls and ν is the kinematic viscosity. Here only the results from Re τ = 590 are shown. Results from the other cases can be found in Lora (2019). The computational domain is shown in Fig. 1a. The domain is extended 5δ × 2δ × 2δ in streamwise, spanwise and wall normal directions respectively. The computational grid is generated using blockMesh utility of OpenFOAM. The dimensionless grid spacings in stream and spanwise direction are Δx + ≈ 40 and Δz + ≈ 25 respectively. Along the wall-normal directions, a stretching law is applied in a way that y + ≈ 1.
The computational domain represents a domain with infinite lengths. Hence, cyclic boundary conditions are imposed in streamwise and spanwise directions. Additionally, no-slip condition is assumed for the walls. In the absence of inlet and outlet, the pressure gradient is necessary to drive the flow. This is achieved by adding an external force to the momentum equation and computing the magnitude of this force from the prescribed bulk velocity (U b ). In this case where Re τ = 590, the bulk velocity is set 0.22 m/s. This corresponds to the bulk Reynolds number of Re b = 10935 based on the half of the distance between the walls (δ).
Second-order central differencing scheme is employed for the convective as well as the diffusive terms. Furthermore, all the simulations are performed using LES method with the Wall Adapting Local Eddy viscosity (WALE) model as a subgrid-scale (SGS) model. The description of this SGS model can be found in Ducros et al. (1998). In the PISO case, the time discretization scheme is set to second-order backward.
In the absence of any obstacle, walls are the only source of turbulence in a flow. In reality, imperfection of walls or small perturbations can trigger the transition process to turbulence flow. However, in this numerical case, very little if any imperfection exists. Therefore, an artificial perturbation is necessary to initiate the process. Here, this is done by adding a random perturbation (random noise) to the domain. This together with the numerical noise and round off errors leads to a fully turbulent flow. The simulations are carried out in two phases. In the first or preliminary phase, the domain is initialized and the flow is simulated for 100 flow passes to reach fully turbulent flow. This allows all the transient process related to the initial conditions to be eliminated. Then velocity and pressure are averaged over additional 1000 flow passes. The time step is adjusted at each time step in a way that the maximum Courant number remains under 0.5. Fig. 1b-d show selected results from the simulations. The results are normalized by the friction velocity u τ . Additionally the DNS results of Moser et al. (1999) are also shown for comparison. The normalized mean streamwise velocity profile u + is plotted in Fig. 1b with respect to logarithmic wall coordinate y + . Computational points close to the wall show good agreement with DNS results up to y + ≈ 10. With increasing wall distance, the overestimation increases and stays approximately constant from y + ≈ 50 on. The deviation from the DNS data also can be observed in Fig. 1c when near the wall (y/δ < 0.2) streamwise turbulence intensities are considered. In other directions, however, the computed normalized turbulence intensities are in a good agreement with the DNS data. Finally, the graph of normalized turbulent shear stress depicted in Fig. 1d shows little difference between the simulations. The results of LES obtained by PISO and RK3 solvers are in good agreement with the DNS data across the complete range of y + .
Overall, the differences between RK3 and PISO are negligible. Furthermore, both solvers in comparison with the DNS data achieve excellent agreement. In conclusion, the implemented solver was performed as expected. However, the most important point is that RK3 solver is found to perform approximately 40% faster than PISO in this case.

Backward step
In the second validation case, turbulent flow over a backward-facing step is modelled employing RK3 solver. This case can also be considered as sudden expansion of the channel. The numerical model is based on the DNS study by Le et al. (1997) and the experimental study of the same case by Jovic and Driver (1994). The flow simulations are carried out at Reynolds number Re h = u 0 h/ν = 5100 with a step height h and u 0 as the mean inlet free stream velocity. This study is based on an expansion ratio of 1.2 (outlet to inlet heights ratio). The computational domain and the boundary faces are illustrated in Fig. 2. The computational grid is generated using blockMesh utility of OpenFOAM and it consists of approximately 850 thousand cells.
The model boundaries are divided into top and bottom wall, the sides, two inlet boundaries and the outlet. The outlet is modelled as zero-pressure outlet. At the bottom wall the no-slip condition is applied setting all velocity components to zero. The top wall is used as symmetry plane. The cyclic boundary condition is applied to the sides of the domain. One of the chal-lenges in this case was the generation of the appropriate inlet boundary conditions. In the original study by Le et al. (1997), at the inlet, flat-plate turbulent boundary layer velocity profile is prescribed. In order to generate the inflow in this study, the inlet is divided into two parts as shown in Fig. 2. This allows the reproduction of the inlet velocity profiles of the original study accurately via special numerical treatments of each boundary face separately.
Similar to the previous study, the second-order central differencing scheme is employed for the convective as well as the diffusive terms. Furthermore, all the simulations are performed using LES method with WALE SGS model. A precursor simulation is carried out to generate a fully developed turbulent flow. Then the flow is simulated for a total time of approximately 4000h/u 0 at constant time steps of Δt ≈ 0.004h/u 0 . During the second phase of the simulation, first and second order statistics are computed.
The flow is separated at the step and reattaches at x/h = 6.05to6.10. This is in a good agreement with the numerical study of Le et al. (1997) where the reattachment computed to occurs at x/h = 6.28. Additionally Jovic and Driver (1994) report the reattachment length of 6.04 to 6.10. The separated zone is shown in Figs. 3 and 4. The primary and secondary vortex are clearly identifiable in the figure. Additionally a third small recirculation is present at the corner (hardly can be seen in the image). The mean flow field retrieved by DNS shows very similar vortex structures.
The normalized friction coefficient C f = 2τ w /ρu 2 0 is computed and illustrated in Fig. 5. It is shown that the computed values are in good agreement with both DNS and the experimental values. Similar agreements are also observed when mean velocities, diagonal and off-diagonal components of the Reynolds stress tensor are compared. Additionally, the spatial twopoint autocorrelations of the velocity components in spanwise direction are carried out. It is shown that the chosen spanwise length of the domain is sufficient.

Flow over a 2D dune
In rivers, several types of bedforms are present. Among these, two-and three-dimensional dunes have large Comparison of the skin coefficient C f between this study and the references impact on the sediment loads and the discharge capacity of rivers. Several experimental and numerical studies have been carried out to expand our understanding of their effects on the flow. It is known that these types of bedforms in the presence of steady unidirectional flows reach a periodic equilibrium shape (see for example Jackson [1976] and Balachandar et al. [2007]). In this study the dune configuration that is experimentally investigated by Polatel and Muste (2006) and numerically studied by Stoesser et al. (2008) and Omidyeganeh and Piomelli (2011) is modelled.
The computational domain is illustrated in Fig    = 20 mm and the dune wavelength is λ = 20k. The width of the domain in the spanwise direction is 8k. Additionally, the water depth is fixed at 4k. In this case, the bulk Reynolds number Re b is equal to 2.5 × 10 4 based on the mean bulk velocity U b = 0.3 m/s and the maximum water depth. The measurement is carried out by Polatel and Muste (2006) at six vertical lines along the dune using laser Doppler velocimetry (LDV). These lines are illustrated in Fig. 6.
The computational grid is generated using the blockMesh utility of Open-FOAM. It consists of 200 × 80 × 60 cells in streamwise, spanwise and wall-normal directions respectively. In the wallnormal direction, the grid spacings are stretched with the ratio of 10 (largest to smallest element ratio). The number of grid points is approximately 10 times lower than the reference LES study of Stoesser et al. (2008). In the streamand spanwise directions, a periodic boundary condition is applied. The water surface is assumed flat and rigidlid condition is imposed. Finally, a noslip boundary condition is applied to the walls.
The simulation is performed with OpenFOAM using explicit third-order Runge-Kutta with fractional-step method. Moreover, Spalart-Allmaras IDDES (Shur et al. 2008) is used for turbulence modelling. Second-order central differencing is employed for the Selected applications of an open-source three-dimensional computational fluid dynamic code in hydraulic engineering Fig. 9 Instantaneous flow field. The contour plots are streamwise velocities and the isosurfaces are visualized vortices using Q-criterion and coloured with streamwise velocities convective as well as for the diffusive terms. The maximum Courant number is set to 0.5 and the simulation is initially carried out for approximately 20λ/U b time units. Then in the second phase, the flow is simulated for 500λ/U b time units and first-as well as secondorder statistics are computed.
Based on computed time-averaged velocities, turbulence intensities and shear stresses, it is found that the results are in excellent agreement with the experiment and the reference LES. As an example, Fig. 7 shows mean streamwise velocities for six vertical lines where measurements are carried out. Fig. 6 shows the streamlines and the location of the reattachment point. The flow separates at the crest and reattaches downstream at x/k = 5.6. This is in agreement with other experimental and numerical studies on dunes (Kadota and Nezu 1999;Stoesser et al. 2008;Grigoriadis et al. 2009;Piomelli and Omidyeganeh 2013). The values are normalized using the bulk velocity. Three distinct zones can be observed in the figure: the separated shear layer (zone 1), the wake layer (zone 2) and the developing boundary layer (zone 3). The shear layer is generated due to the flow separation and Kelvin-Helmholtz instabilities. Further downstream beyond the reattachment point, a wake layer is produced and extended beyond one dune length. This layer is produced by the vortices in the shear layer advected further downstream. Finally, the development of the boundary layer near the wall generates the upward displacements of the vortices.
One of the advantages of scale-resolving simulation (SRS) over RANS is the possibility of studying the instantaneous flow field and turbulence structures. These provide valuable insight into the hydrodynamic processes and physical mechanism involved in such flows. Fig. 9 shows a snapshot of the instantaneous streamwise flow velocity and the vortices. Studying the flow field reveals that in the separated shear layer spanwise vortices (rollers) are generated because of the Kelvin-Helmholtz instabilities. These vortical structures swept downstream and they undergo secondary instabilities. Consequently, the spanwise vortices are deformed into hairpin-shaped vortices. As these vortices move upward and reach the water surface, they produce phenomena known as boils. Finally, beyond the reattachment point and near the walls, streamwise vortices or streaks are forming.

Secondary currents in a trapezoidal channel
Numerous natural and man-made channels are of trapezoidal forms. Prediction of the secondary currents in a straight trapezoidal channel by numerical modelling using RANS family Selected applications of an open-source three-dimensional computational fluid dynamic code in hydraulic engineering of turbulence models is difficult due to the anisotropy of the turbulence. On the other hand, these secondary flows are important because of their effect on the primary mean flow, shear stress and formation of bed forms. In this study, a straight trapezoidal channel is modelled using an advance turbulence model to accurately model primary and secondary flows. This channel represents a straight part of Mur River at laboratory scale (1/40). The results from this simulation are used as an input for a numerical model where approach flow to ROR plants are investigated. The computational domain and the coordinate system of the channel is illustrated in Fig. 10. In this study, x-, y-and z-axis represent streamwise, spanwise and bed-normal directions, respectively. The water depth, H = 0.16 m and width of the channel, B , is approximately 4.6H with the bank angle of 33 degree. In the streamwise direction, the domain is extended to D = 8H. The mean streamwise velocity is 0.075m/s leading to the bulk Reynolds number of Re = 12000 based on the water depth. The computational grid is generated using the blockMesh utility of OpenFOAM consisting of approximately 2 million hexahedra elements. The dimensionless grid spacings based on the bulk wall shear stress in stream-and spanwise directions are Δx + = 50 and Δy + = 30 respectively. In the bed-normal direction, the grid consists of 50 layers where the grid size is coarsened towards the water surface where the ratio of the largest to the smallest element equal to 8. This ratio puts the first grid point at z + ≈ 2 based on the time-averaged values.
The simulation is performed with OpenFOAM using explicit third-order Runge-Kutta with fractional-step method. Moreover, Spalart-Allmaras IDDES is used for turbulence modelling. Second-order central differencing is employed for the convective as well as for the diffusive terms. Rigidlid boundary condition is imposed at the water surface. No-slip boundary condition is applied at the bed and the sidewalls. The simulation is carried out in three steps. First, fully turbulent flow is generated using Divergence Free Synthetic Eddy Method (Poletto et al. 2013). The general idea behind the method is to continuously generate and inject synthetic eddies at the inlet in a way that generated coherent flow structures persist in the domain. In the second part, the boundary conditions in streamwise direction are changed to cyclic and the simulation is performed for 20 flow passes. In the final step, the simulation is continued for approximately 100 flow passes and time averaging is performed. Fig. 11 shows contour plots of the normalized streamwise velocity on a spanwise plane. Only half of the channel is shown due to symmetry of the domain. Moreover, the results obtained by Standard k − ε (SKE) turbulence model (Launder and Spalding 1974) are also shown for comparison. It is shown that the results differ, especially near the sidewalls. In the results obtained from IDDES, the bulging of the velocity is observed toward the corner of the channel. This is also observed in an experimental study by Tominaga et al. (1989) for trapezoidal channel. Similarly, in a numerical study by Wright et al. (2004), where flow in a straight trapezoidal channel is modelled using RANS and LES, only results from LES showed bulging of the velocity toward the corner of the channel.
Selected applications of an open-source three-dimensional computational fluid dynamic code in hydraulic engineering   Fig. 12. These cells are also observed in experimental as well as numerical studies of Tominaga et al. (1989) and Wright et al. (2004), respectively. The bulging of the velocity that was previously observed in Fig. 11 is generated due to the two contra-rotating vortices at the corner. In comparison to the flow in a rectangular channel, in a trapezoidal channel an additional vortex is present near the free surface and the sidewall. Overall, as the sidewall angle decreases, flow characteristics are changing from rectangular to trapezoidal channel.
At the centre of the channel, streamwise velocity profile is following the logarithmic law up to the water surface. Similarly, the turbulence intensities are plotted in Fig. 13 along a vertical line at the centre of the channel. Additionally, values computed by the semi-theoretical expression of Nezu and Nakagawa (1993) for turbulence intensities are added. Overall, there is a good agreement between the numerical model and the semi-theoretical values especially in z/H > 0.2 range with exception of the turbulence intensities in vertical direction near the water surface where the values approach zero due to the boundary condition. Although, in z/H < 0.2, the results differ due to the characteristics of the semi-theoretical expressions, when numerical values are compared to the measured values provided by Nezu and Nakagawa (1993), it is evident that the numerical model has accurately predicted the turbulence intensities. In addition to turbulence intensities, turbulence kinetic energy is also plotted in Fig. 13 along a vertical line at the centre of the channel. It is shown that the computed values are in excellent agreement with the expression proposed by Nezu and Nakagawa (1993). On the other hand, the values from Nikora and Goring (2000) are larger than the other two especially for z/H > 0.2.
In conclusion, the numerical model predicted flow velocities, turbulence intensities and turbulence kinetic energy accurately. It is shown that RANS turbulence models based on eddy-viscosity hypothesis (e.g. SKE) have shortcomings in predicting secondary flow in a straight trapezoidal open channel. Additionally, it is demonstrated that three secondary flow cells exist near the sidewalls. These cells are different from the ones observed in a flow in a rectangular open channel. Finally, a similar study is carried out for a trapezoidalrectangular channel with rough walls based on the experiment by Blanckaert et al. (2011). In this study, the channel is modelled using LES and roughness elements are generated using a method proposed by Stoesser (2010). Details of this study can be found in Steger (2019).

Numerical modelling of ROR approach flow
Design and operation of ROR plants is far from straightforward. The design must be flexible and robust to maximize the electricity production in a highly dynamic environment. This requires detailed and comprehensive study of the flow. Among many aspects, an optimum design of the approach channel and intake structure are especially important. The design procedure must consider elimination of unfavourable flow at the intakes and reduction of energy losses to a minimum level. The investigations are usually carried out using physical model test and more recently by numerical models. Physical model tests are well established and refined over decades. However, numerical models are relatively new in modelling a complex approach flow. Hence, it is necessary to evaluate its performance and identify suitable key parameters involved in modelling. The present study is focusing on these aspects by means of a case study. In principal, turbines are designed to work at their full capacity for the uniform flow or ideal intake condition. However, in reality an ideal flow condition is hardly achieved. In high-head hydropower plants (HPP), the intake structure is far away from its turbines and usually flow upstream the turbines is well developed. On the other hand, in the low-head HPPs, the intake structure is short and flow condition becomes important.
Unfavourable flow con-ditions can lead to vibration-induced damages and reduction of energy production. Therefore, early detection of possible flow problems and optimization of the intake is crucial in the design phase.
The case study resembles a typical ROR block layout with bed enlargement where the powerhouse is located on the right side at the enlarged side of the river with two Kaplan (pit) turbines. Additionally, three weirs structures are located on the axis of the river. The general layout of the physical model is illustrated in Fig. 14. This project is chosen for this study because of its similarity to other ROR projects in several aspects. Furthermore, availability of a physical model provides the possibility to compare measured values to those computed by the numerical model. Successful numerical modelling of this case facilitates the future investigation for other similar projects and possibly optimization of the designs via numerical models.
The physical model study is carried out at the Hydraulic Laboratory of the Institute of Hydraulic Engineering and Water Resources Management of Graz University of Technology. After completion, wide range of studies are carried out including evaluation and optimization of the approach flow. The physical model was constructed based on Froude similarity law with 1:40 scale. The total length of the model is approximately 18.5 m with 2.4 m width at the inlet of the model. The powerhouse and the weirs are located approximately 6.5 m from the model's inlet. In order to improve inflow condition, a flow straightener was placed at the inlet. The intakes are made from a transparent acrylic glass. Openings are placed at the top of the intakes for measurement device (Fig. 15). In this study, the load case where both turbines are in operation is studied with the design discharge of 20.25 l/s (205 m 3 /s in prototype).
The flow velocities are measured with Acoustic Doppler Velocimetry (ADV) device. A side-looking Vecterino+ probe from Nortek AS is used for the measurement.
The device has one transmitter and four receivers capable of measuring the velocity components with high frequency. The probe was orientated in a way that the large component of the velocity is approximately perpendicular to the transmitting direction. The measurement was carried out for 100 points for each intake, 200 points overall (Fig. 15). At each point flow velocities were recorded for 90 s. It was verified that this period is adequate for computation of mean values by performing a series of tests. The raw data is then processed using WINADV (Wahl 2000) software. The correlation (COR) and Signal to Noise Ratio (SNR) thresholds are set to 70 and 15 respectively. Additionally, de-spiking filter based on phase-space threshold is applied. The method is proposed by Goring and Nikora (2002) and further modified by Wahl (2003). Finally, the results are visualized using Python. The aim of the numerical model was to evaluate the accuracy of the numerical model in predicting flow inside the intakes. This is done by comparing computed results to the measured values. Furthermore, series of sensitivity studies are carried out to identify important numerical parameters that significantly influence the results. In order to realize the objectives of the study, a numerical domain is generated, representing the physical model test. The domain is extended from the runner of the turbines to the location of the inlet of the physical model. Three sets of computational grids with different levels of refinements are created. They consist of all hexahedra elements along with refinement at walls to place the first point within the applicability of the wall functions. These grids are used to carry out a grid convergence study.
The simulations are carried out using OpenFOAM. The numerical domain and its corresponding boundary faces are shown in Fig. 16. At the walls, noslip condition is assumed for the velocity followed by zero gradient for pressure. Appropriate wall function is then applied for turbulence parameters. The wall functions vary depending on the turbulence model used. At the outlets, pressure is set to zero along with zero gradient condition for other parameters. Assuming that the variation of the water level in the domain is negligible, rigid-lid condition is employed for the water surface. Finally, velocity and turbulence parameters are imposed at the inlet of the model. These values are computed via precursor simulations of a periodic trapezoidal channel.
The incompressible RANS equations are solved using Semi-Implicit Method Selected applications of an open-source three-dimensional computational fluid dynamic code in hydraulic engineering  (Patankar 1980). Three RANS turbulence models, SKE, Spalart-Allmaras or SA (Spalart and Allmaras 1992) and Shear Stress Transport or SST (Menter et al. 2003), are employed to investigate the role of turbulence models. Spatial discretization scheme for convection term is set to secondorder upwind for momentum as well as turbulence equations. In order to ensure convergence, besides monitoring equation residuals, several probes are defined and their values monitored during the computation.
A grid convergence study is carried out to determine the effects of grid refinement on the results. Three sets of grids, Coarse, Medium and Fine, are generated with approximate global size of 2, 1 and 0.5 cm respectively. Moreover, five points for each intake, one for each corner (B2, B9, I2 and I9) and one approximately at the centre of the sections (point F6). The Grid Convergence Index (GCI) is then computed for the aforementioned points by Approximate Error Spline (AES) method proposed by Celik et al. (2005). It is shown that in general GCI values are low for all the grid sets. Furthermore, GCI values are smaller in the right intake (intake 1) in comparison with the left intake (intake 2). Finally, it was concluded that a coarse computational grid is sufficient for preliminary and optimization studies. However, a finer mesh is necessary to obtain the final solution and to reduce errors related to the grid size. Therefore, a medium size grid with global size of 1 cm is chosen for the rest of the studies as the difference between the medium and fine grid sizes are negligible.
In hydraulic engineering the physical model tests play an important role in studying hydrodynamics and sediment transport processes. Often a physical model in the laboratory is a scaled down representation of a large prototype using similarity laws. Considerable differences can be observed between the up-scaled results from the model test and the measured values from the prototype due to scale effects, if an appropriate value for the scaling factor has not been chosen. A series of numerical tests are carried out to investigate the scale effects. In these numerical tests, results from the laboratory scale model are compared to the prototype and it is found that the differences between the computed velocities Selected applications of an open-source three-dimensional computational fluid dynamic code in hydraulic engineering are negligible. Furthermore, in these tests, a wide range of roughness values, corresponding to concrete material, are considered for the walls. These studies are carried out for both laboratory scale and prototype models. It is found that although considering roughness in the computation has a great impact on the magnitude of the shear stresses, it does not influence the velocity distribution inside the intakes.
In a single-phase internal flow, it is common that the boundary conditions are velocity inlet and pressure outlet. In this set up, the pressure at the outlet is defined (usually 0 in case of one outlet) and velocity along with turbulence parameters are imposed on the inlet. Usually the values are defined as constant using empirical formulas or they are obtained via precursor analy-sis. The velocities can be computed for the entire section with assumption that it is constant over the entire section. For turbulence values (e.g. kinetic energy and dissipation rate), these values can be computed using mathematical expressions which are available in the literature. Although imposing values based on the description above is sufficient due to prominent geometrical features downstream, in all the studies here the inlet values are generated by precursor simulations using periodic channel. Besides imposing accurate values for the inlet, the location of the inlet with respect to the zone of interest is highly important. In order to investigate the effect of inlet's location on the results, three numerical tests have been carried out. In the first model, the inlet is placed immediately before the forebay. In the second and third model, the trapezoidal approach channel is extended upstream by eight and sixteen times the channel's water depth respectively. It is found that as far as the velocity distribution inside the intakes is of concern, the inlet of the model can be placed at the beginning of the forebay. This can reduce the size of the computational domain and consequently decrease the computational resources required to model this type of HPPs. In addition, it indicated that the geometrical features that are present upstream of the intakes have strong influence on the flow. Furthermore, if the results beyond the intake are required, the approach channel must be included and extended at least to eight times the water depth in the channel.
As discussed above, as long as the distribution of the velocity in the intakes are a concern, the results are not significantly affected by grid refinement, scale and roughness effects nor the location of the inlet. One of the important parameters in modelling turbulent flows is the turbulence modelling. In order to investigate the role of the turbulence models, three commonly used RANS models SKE, SA and SST models are employed. The results obtained by the SST model are presented and discussed here because of the better performance of this model in predicting the flow. Fig. 17 shows the contour plots of the out-of-plane velocities at the measured section. The figure compares the computed values by the SST model to those that are measured. In intake 1, computed velocity distribution and their magnitudes are in good agreements with the experiment. This is also shown in Fig. 18 where out-of-plane velocities are plotted along selected vertical lines. In addition, the agreement between the computed and the measured values are found to be good for the other components of the velocity. Finally, all the turbulence models predicted similar flow distribution at this section.
In intake 1, the most notable feature is the sharp gradient of the velocity in vertical direction where high out-ofplane velocities are located at the upper regions reaching up to 1.4 times the average velocity. The values then reduce to 0.8 times the average out-ofplane velocity at the lower zones. It is shown in an additional study that reducing the slope before the intakes can significantly reduce the velocity gradi-Selected applications of an open-source three-dimensional computational fluid dynamic code in hydraulic engineering ent. Finally, on the right side of the intake 1 near the wall (z/H ≈ 0.4), a low velocity region is predicted by the SST model. This zone corresponds to the vortical structure that begins from the water surface and ends approximately at the study section. This vertical structure is visualized (No. 5) in Fig. 20 via isosurfaces of the Q-criterion.
In intake 2 (left intake), overall the computed flow velocities are in an agreement with the measured values, however, there are differences especially on the left side of the intake, where predicted values are lower than those which are measured. Furthermore, it is at this zone where different turbulence models predicted different flow distributions. On the other hand, the flow on the other side of the intake is found to be similar to intake 1. Fig. 19 shows the computed and measured out-of-plane velocities along the selected vertical measurement lines. It is shown that on the right side, the results from the numerical model have predicted the out-of-plane velocities accurately. However, on the left side of the wall up to line L3, the numerical model underestimated the out-ofplane velocities. This is partly due to prediction of higher swirling flows and higher in-plane velocities near the left wall.
In order to explain the unique characteristic and the origin of the flow in intake 2, vortical flow structures are visualized in Fig. 20. It can be seen that at least four distinct vertical structures are present. Vortical structure 1 is a clockwise swirling flow that is generated due to the flow from the left side passing over the divider wall and flowing into the left intake. In fact, it is found that almost the entire flow in intake 2 is originated from the left side of the intakes. Moreover, due to the high velocity flow from the left side an additional vortical structure is generated above the divider wall. This is indicated by number 2 in Fig. 20. The submerged pier on the left side of the intakes also generates a vortical structure (No. 3). It is shown that at this place, flow must change direction completely and as it passes this structure, it separates and generates a longitudinal vortex toward the intake. Finally, a vertical vortex is predicted on the left side (No. 4), which originated from the water surface and ended midway through the entrance of the intake. This free surface vortex is also observed in the experiment, however, the location and presence of this vortex varied significantly through the experimental tests.
Typically, each turbine manufacturer has its own set of criteria for acceptable flow condition. Nevertheless, Fisher and Franke (1987) proposed a set of guidelines based on their own experiences and three turbine manufacturers known as Fisher-Franke (F&F) guidelines. Later, Godde (1994) carried out experimental investigations on the bulb turbines intake. The author applied and compared existing criteria to the physical models. According to these guidelines, flow must be evaluated at a control section. The control section is usually located after the trash-rack and upstream of the turbine due to its accessibility for measurement in the physical model tests. These criteria are applied to the velocity values that are obtained from numerical as well as experimental study. Fig. 21 shows flow evaluation based on F&F criterion for the study section. Overall, in intake 1, the experimental values are within the range proposed by Fisher and Franke (1987). Moreover, the results from the numerical model show slightly higher values in the upper range. On the other hand, in the lower range, the values are in a good agreement with the experiment and therefore within the range of the criterion. In contrast to intake 1, in intake 2, although experimental values are within the range, the numerical values are significantly lower for a large area. This is aligned with the contour plots of out-of-plane velocities in Fig. 17. Finally, in the upper range both experimental and numerical values are slightly deviating from the criterion. A series of tests have shown that at this section, the flow is not fully developed. Hence, a new location for the reference section is proposed. It is shown that a section further downstream is better suitable for flow evaluation.
Another commonly used criterion is deviation of the quadrant's discharges from the ideal conditions. Typically, the reference section is divided into quadrants (Fig. 15) and flow rate is computed for each quadrant. Then the deviation of the flow rate from the ideal condition is computed and plotted. Fig. 22 shows this condition for both intakes. In intake 1 not only the predicted values are in good agreement with the experiment but also the deviations are relatively low (< 4%). On the other hand, in intake 2, deviations that are computed Selected applications of an open-source three-dimensional computational fluid dynamic code in hydraulic engineering

Fig. 22
Deviation of the quadrant's discharges from the ideal condition by the numerical model are reaching up to 5%, significantly higher than the measured. It must be noted that part of these differences is due to the underlying assumption in calculation of the discharges in the experiment where the velocities outside of the measured sections are assumed to have the values of the outermost measurement points. However, in the numerical study it is shown that a significant velocity gradient exists near the walls.
In conclusion, it is found that if the modelling of the flow and distribution of the velocity inside intakes are considered, the effect of grid refinement, roughness and inflow condition are negligible. On the other hand, a suitable turbulence model is crucial as it is found that changing the models significantly affects the prediction especially on the right intake where a complex flow was expected. Furthermore, there is a good agreement between the numerical and experimental results. Finally, it is found that it is better to evaluate the flow at a section as far as possible from the intake's entrance as the flow near the entrance is not fully developed.

Summary and conclusion
In this paper, a series of studies are presented where a three-dimensional open-source computational fluid dynamic code OpenFOAM is employed to carry out the numerical simulations. These studies cover a wide range of applications: from implementation of a new algorithm to modelling of a practical hydraulic engineering problems. Via access to the source code, a new solver is implemented. This solver is then validated via two classical turbulent flow cases. It is shown that the solver significantly reduced the computational time in comparison with PISO solver. In another case, flow over a twodimensional dune is modelled using a novel hybrid approach. The results are compared to experimental and LES of the same case. Furthermore, a brief description of the hydro-dynamical processes is presented. This study is followed by the numerical study, where the same hybrid approach is used to model secondary currents in a straight trapezoidal channel. The secondary flow cells are identified and the results are compared to the ones observed in the experimental studies. Additionally, it is found that the novel hybrid approach is capable to accurately model the flow with relatively lower computational costs in comparison with LES. Finally, in the last part of this paper, approach flow of a ROR plant is studied numerically. A series of parameter studies are carried out to identify important numerical parameters which affect the predictions with respect to the flow inside the intakes. It is found that turbulence models play an important role in this case and must be chosen carefully. Moreover, a numerical model can be used for preliminary and optimization design phases. However, it is recommended to validate the results from the numerical model with the measured values from the physical model.
In conclusion, although the learning curve of such tools are somewhat different from their commercial counterparts, these tools have the advantage that the source code is freely available, hence, implementation of new algorithm is possible. In addition, these codes do not have licensing costs that are usually associated with the commercial codes. Finally, CFD codes like OpenFOAM typically have large user base across most areas of engineering and science and, therefore, it is constantly being updated and additional capabilities are added.
Funding Open access funding provided by Graz University of Technology.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://