*Article* **Simulation of a Hydrostatic Pressure Machine with Caffa3d Solver: Numerical Model Characterization and Evaluation**

#### **Rodolfo Pienika 1,2,\*, Gabriel Usera <sup>1</sup> and Helena M. Ramos <sup>2</sup>**


Received: 21 July 2020; Accepted: 25 August 2020; Published: 28 August 2020

**Abstract:** The Hydrostatic Pressure Machine (HPM) is a novel energy converter for micro and pico hydropower that becomes very suitable for installation in channels with very low head, where conventional hydraulic turbines are inadequate or too expensive. Although this technology has been studied through several experimental tests and also by numerical simulations, open source flow solvers have not been used yet. The research team on Computational Fluid Mechanics of IMFIA-Universidad de la República (Uruguay) has been developing a CFD open source solver named caffa3d, which has obtained great results in a few international challenges, although it has not been used yet for free surface flows or turbomachinery simulations. The present work shows the contributions made within caffa3d in order to enable its use for simulating a HPM. The Large Eddy Simulation (LES) method is used to model the turbulence structures of the flow. Sliding Mesh (SM) and Volume of Fluid (VOF) methods were chosen respectively to resolve the rotation of the wheel and the position of the free surface. The SM module was already validated in the past, but the VOF module needed to be validated in the present work through the simulation of free surface over a semicylindrical dam. Finally, the performance of a small 12-straight-blade HPM was simulated with caffa3d, with quite satisfactory results. Some issues of the solver yet need to be solved before other HPM with more complex designs could be studied.

**Keywords:** hydrostatic pressure machine; micro hydropower; CFD; open source; sliding mesh; volume of fluid; caffa3d

#### **1. Introduction**

Hydropower, being one of the oldest energy sources developed by humankind, continues to also be one of the most prominent ones, in terms of costs, social, and environmental impacts. A renovated interest in micro and pico hydropower has been developing in recent years to accommodate the increasing electricity demand avoiding the environmental impacts of large power stations [1]. Very few non exploited sites along the entire world are suitable for medium and large hydropower, considering the social and environmental negative impacts. In the opposite side, micro to pico hydropower generates minimum impact on the environment in addition to when they are installed in already existing infrastructures such as water distribution systems, irrigation dams or water channels. Conventional hydraulic turbines become very expensive (especially in relation to the total cost of the project) when it comes to very low heads and low power ratings [2], and most of the suitable micro hydropower converters present poor efficiencies or are still in a development stage such as the ones presented in [3–5].

The Hydrostatic Pressure Machine (HPM) is a novel energy converter first introduced by [6] inspired in the ancient water wheels used within mills [7], but with the main difference of extracting energy from the pressure of the water flow instead of the kinetic or potential energy. Among the recent publications, there are several studies on waterwheels' performances and improvements, but a great part of these refer to stream or gravity wheels [8–12] and only a few refer to HPM [13–18] that operate in a different way and should be distinguished from the formers. The HPM becomes very suitable for installation in water channels with very low head, where conventional hydraulic turbines or hydrokinetic turbines are inadequate or just too expensive because of the low head or power ratings. The installation of a HPM requires very little modifications to the existing infrastructure, and has proven to work with acceptable efficiencies (up to 80%), which make them the best choice for very low head waterways like irrigation or waste water channels. The principle of operation of a generic HPM is explained in Figure 1. The presence of the HPM covering almost the entire water passage of a rectangular channel section (leaving only the required gap between it and the channel walls and bed) generates a water level difference between upstream and downstream (just like a sluice gate). Then, like the sluice gate, the blades of the HPM receive a force due to this level or pressure difference, which generates a torque at the rotating centre and consequently a shaft mechanical power.

The rotational speed in this machines is very low, so, on one hand, there is a technological challenge to convert the mechanical to electrical power, but on the other hand and in combination with the narrow bottom opening, enables sediment transport and fish passage through it.

**Figure 1.** Explanation of the operation of a generic HPM, taken from HYLOW Internal Task Report 2.3 [19].

The operation of some HPM designs has been assessed during the last years, mostly within the HYLOW project (of which CERIS, from IST of Lisbon, was part of), by means of experimental tests in reduced scale models [13] and field tests [14]. In addition, some numerical assessment with commercial software has been done in order to better understand the optimum configuration for the performance of the HPM [15]. The highlighted modifications introduced to the design shown in Figure 1 are the alteration of the bed of the channel below the HPM in order to reduce the leakage flow, and the inclination of the blades remaining straight or radially twisted (see Figure 2), in order to facilitate the filling and emptying of the buckets. This latter improvement would also reduce the impact when the blades enter the water and would make the torque more continuous along the rotation. However, some doubts regarding the shape optimization of the blades to extract the most of the energy still remain. In fact, there are currently several ongoing projects that try to search for the best configuration, not only of the geometry of the wheel [16], but also of the geometry of the channel, water levels and type of generator [17,18].

**Figure 2.** HPM with radially twisted blades and modified bottom shroud, taken from HYLOW Internal Task Report 2.3 [19]. In the present work, only the modification of the bottom shroud was modeled, maintaining the radially straight blades as in Figure 1.

Regarding the use of numerical tools for performance assessment of turbomachinery, the majority of academic researchers use commercial 3D flow solvers, while only few are using open source and freely available solvers. Within the latter, OpenFOAM is by far the most known and used around the world. In addition, a small but growing group of researchers from Uruguay have been developing another CFD solver named caffa3d [20]. It has evolved from a family of 2D flow solvers introduced by Ferziger and Peric [21], it is also open source, freely available, and used by researchers all around the world. The flow solver caffa3d is an implementation of the finite volume method in Fortran 90; it uses block structured body fitted grids, it can be combined with the immerse boundary condition method, it can be parallelized through the domain decomposition under a distributed memory model using the MPI library, and include a few methods for modeling turbulence like RANS (Reynolds Average Navier–Stokes) and LES (Large Eddy Simulation). The research team on Computational Fluid Mechanics of School of Engineering from Universidad de la República (Uruguay) has successfully participated with caffa3d in few CFD international challenges, working with an immerse boundary condition in hemodynamic simulation to predict the rupture of an intracranial aneurysm [22,23] and with the actuator line model in the simulation of wakes behind wind turbines [24]. Other uses of the solver include simulation of wind around buildings, pollutant dispersion, numeric wind tunnel, wind farm evaluation, and more [25]. However, it has not been successfully used yet for the simulation of rotating machinery and free surface flows.

Although the original version of the caffa3d solver from 2008 already included approaches for the solution of this kind of flows, they have not been used since then, and therefore became out of date within the last updates of the solver. Then, a big effort was done in updating the specific part of the code to make it work properly. This paper presents the most important parts of the updating and the results of preliminary simulations of straight-blades HPM in order to validate the solver. The next step will be to modify the geometry of the wheel and assess its optimum shape.

#### **2. Materials and Methods**

The numerical simulation of the operation of a HPM is obviously transient, and the main challenges are: the rotation of the wheel, the presence of a free surface and the generation of an adequate mesh. Further description of the solver caffa3d can be found in [20,25]. Other configuration parameters like boundary conditions and time step are also covered within this section. Finally, the procedure for the calculation of power generated by the HPM is stated.

#### *2.1. Rotation of the Wheel*

The rotation of the wheel can be simulated by means of different approaches, of which the multiple reference frames (MRF) and sliding meshes (SM) are highlighted.

In the MRF approach, the solver runs in a relative reference frame for the rotating part of the domain (working with relative velocities and adding the transport and Coriolis terms to the momentum equations) and in an absolute reference frame for the stationary part of the domain (working with absolute velocities and no addition of terms to the momentum equations). In MRF, all the blocks in the domain remain fixed and there has to be a steady flow condition at the interface.

In the SM approach [21], the entire flow is solved in the absolute reference frame for two separate parts of the domain (the moving and the stationary), with no addition of terms to the momentum equations. Because it is solved in the absolute frame, the mesh of the rotating part of the domain must also rotate, and some fluxes must be corrected because of mesh movement and some connectivity information through the interface must be updated. It is an inherently transient simulation, and captures the actual position of the rotor at any time, making it the best approach to solve the simulation of a HPM.

In the interface between the rotating and the stationary meshes, there is a perfect match between cells only for null relative displacement angle *α* (Figure 3). For other values of *α*, every cell in one side of the interface relates to two cells in the other side.

**Figure 3.** Sketch of the sliding interface, showing the relative displacement angle *α* between stationary blocks (painted) and rotating blocks (white).

The momentum balance equation for each moving block must consider the local time derivative of the absolute velocity field (*v*-*A*), and the relative velocities (*v*-*<sup>R</sup>*) for the mass fluxes, in a way as [20]:

$$\int\_{\Omega} \rho \frac{\partial\_R \vec{v\_A}}{\partial t} d\Omega + \int\_{\mathcal{S}} \rho \vec{v\_A} (\vec{v\_R} \times \mathfrak{n}\_{\mathcal{S}}) dS = \int\_{\Omega} \rho \vec{\mathfrak{g}} d\Omega + \int\_{\mathcal{S}} -p\mathfrak{n}\_{\mathcal{S}} dS + \int\_{\mathcal{S}} \mu (\nabla \vec{v\_A} + \nabla \vec{v\_A}^\top) \mathfrak{n}\_{\mathcal{S}} dS \tag{1}$$

where subindex A refers to the inertial reference system and subindex R refers to the reference system that rotates with the blocks, *ρ* is the density, *μ* is the dynamic viscosity, *g* is the gravity force, and *p* is the pressure.

#### *2.2. Free Surface Flow*

For the solution of the free surface flow, there also exists several methods, some of which are explained by Ferzger and Peric [21]. Free surface simulation methods can be classified as Interface-Tracking where each fluid occupies a different part of the grid, and the grid changes each time step to accommodate to the free surface geometry, as Interface-Capturing where the grid is fixed and the free surface geometry is found by different approaches, or as an hybrid method. Among the Interface-Capturing methods, the most popular one is the Volume of Fluid (VOF), developed by Hirt and Nichols [26] with latter updates and modifications [27]. In the VOF method, a scalar (*Vf*) is used to assess the fraction of volume of each fluid that occupies each cell volume at each time. For *Vf* =1, the cell is filled with water, for *Vf* =0, the cell is filled with air, and, in the middle, there is an interface between both fluids. Fluid properties such as density and viscosity are derived from the value of *Vf* for each cell at a given time, and the densities and viscosities of water and air:

$$
\rho = V\_f \rho\_w + (1 - V\_f)\rho\_a \tag{2}
$$

$$
\mu = V\_f \mu\_w + (1 - V\_f)\mu\_a \tag{3}
$$

Then, the corresponding transport Equation (4) is solved together with conservation equations for mass (5) and momentum (6), to find the geometry of the free surface at each time step:

$$\frac{\partial V\_f}{\partial t} + \nabla(V\_f \vec{v}) = 0\tag{4}$$

$$\frac{\partial \rho}{\partial t} + \nabla(\rho \vec{v}) = 0 \tag{5}$$

$$\frac{\partial \rho \vec{v}}{\partial t} + \nabla (\rho \vec{v} \vec{v}) + \nabla p = \rho \vec{g} + \vec{F}\_{\sf s} + \nabla \left[ \mu (\nabla \vec{v}\_{A}^{\*} + \nabla \vec{v}\_{A}^{\*T}) \right] \tag{6}$$

where *ui* is the component of the velocity vector *v* in the direction of *xi*, and - *Fs* is the force due to surface tension.

The solver caffa3d uses the VOF method together with a CICSAM approach (Compressive Interface Capturing Scheme for Arbitrary Meshes) for flux calculation [28,29], and the inclusion of the surface tension force [30,31]. In addition, in order to avoid non physical values of *Vf* , it is implemented a predictor-corrector procedure that resets these values to zero or one [28].

#### *2.3. Model Geometry and Meshing*

The analysed HPM (see Figure 4) is a 12-straight-blade-wheel with 150 mm inner diameter, 450 mm outer diameter, and 235 mm width. The ratio between the channel width and the wheel width is 1.26, meaning that the distance from the wheel to the channel's lateral walls is 30 mm. The thickness of the blades is 2 mm and the tips have a semi-spherical shape to reduce the risk of boundary layer separation. The shroud has a modified arc-shaped bottom that is attached to at least one full bucket every time, looking for a reduction of leakage flow. The gap between the blades and the bottom shroud is 10 mm (smaller gap could not be modeled accurately because it affected the simulation stability), while the gap with the lateral wall is 30 mm. Due to the straight blades, the geometry and thus also the flow, have a symmetry plane in the midsection of the channel. This enables the definition of a symmetry boundary condition at the midsection and thus only half of the entire domain needed to be modeled, reducing the computational time.

**Figure 4.** Computational model of the analysed HPM.

The full mesh, generated within caffa3d solver, is composed of three parallel bands of blocks (seventy-five in total): two bands for the half width of the HPM and the remaining for the width of the lateral gap. The only blocks that cover the entire width of the domain (half of the channel) are the two ring-shaped blocks corresponding to the sliding interface. Each of the bands including parts of the HPM are composed of sixteen structured grid blocks of H-type, O-type and C-type [21]: twelve body-fitted blocks for the wheel (one for each sector around blades) of the C-type, one H-type block surrounding the ring blocks, and three orthogonal blocks for the inlet, outlet and upper regions of the domain. These latter blocks are long enough to avoid interference between the wheel and the inlet and outlet sections. Moreover, both sections are located 5 times the outer diameter away from the periphery of the wheel. The band of blocks covering the lateral gap is similar to the others, but with more blocks to model the space adjacent to the hub and blades of the HPM, totaling forty-one blocks for this band. Currently, the band of blocks containing the lateral gap is modeled as a narrow passage of constant width along the entire length of the channel. However, in the real HPM, the lateral gap is of course much smaller than the difference between the channel and wheel width, and this is realized by positioning a wall normal to the flow at the plane of the rotation axis. This was not done in the present simulation, just for simplifying the geometry modeling and the simulation itself, but, in the future, this wall will be included in the model with help of the immerse boundary conditions module implemented in caffa3d.

The solver is domain parallelised within regions which are integrated by one or more grid blocks [20,25]. For the implementation of SM approach in caffa3d, the pair of blocks related by the sliding interface must be within the same region. The rest of the blocks are arranged in fifteen regions searching for a rough equality of element count among them all.

The grid is refined near the walls (to correctly model the viscous sublayer) and in a vertical direction around the expected position of the free surface, in order to capture its shape accurately. Because the free surface will occupy variable parts of the interior of the wheel as it rotates, the grid in this region also needs to be refined.

The final version of the mesh has 1.2 million volume cells,strongly concentrated around the HPM, as can be appreciated in Figures 5 and 6. The smallest cells are 0.3 mm in length and correspond to cells close to the tips of the blades.

**Figure 5.** Computational mesh of the entire domain with dimensions in meters.

**Figure 6.** Computational mesh of the region surrounding the HPM (**left**) and around a blade (**right**) with dimensions in meters.

#### *2.4. Boundary Conditions*

Symmetry boundary condition is applied to one of the lateral walls of the domain (that corresponds to the midsection of the actual channel), while the other lateral wall is treated as a non-slip wall, as also are the surfaces of the wheel (blades and hub) and the bottom of the channel.

For the inlet section, it is assigned uniform velocity for the water region and zero velocity for the air region, thus imposing the inlet water flow (discharge).

Zero manometric pressure was assigned to the top of the domain, to simulate the atmosphere, and also to the outlet section of the channel, to simulate a free discharge.

#### *2.5. Other Configuration Parameters*

The rotational speed of the HPM at the start of the simulation is increased smoothly from zero to the nominal value during the first 50 steps, in order to prevent instabilities during the start-up of the simulation. In this simulation, the rotational speed is fixed in its nominal value, but the solver has the possibility to calculate it from the interaction with the fluid if such would be the case.

Regarding the time step, for the SM method, it would be sufficient to be less than the time that takes a cell to rotate one position. However, the VOF method imposes higher restrictions in order to satisfy the Courant–Friedrichs–Lewy condition. This restriction implies that, for any given cell, Equation (7) must verify, where *u* is the velocity of the fluid, Δ*t* is the time step, and Δ*x* is the smallest dimension of the cell. *<sup>u</sup>*Δ*<sup>t</sup>*

$$\frac{u\Delta t}{\Delta x} \le 1\tag{7}$$

The solver is able to adapt the time step to accomplish the latter condition. The simulation is initiated with a time step of 5 × <sup>10</sup>−<sup>5</sup> s, but it is reduced whenever necessary. Actually, the ideal time step would have been close to 1.5 × <sup>10</sup>−<sup>5</sup> s.

The large eddy simulation (LES) technique, with the Smagorinsky model [32], is attached to the solver to model turbulence. The size of the cells adjacent to wall surfaces in normal direction (*y*) needs to be controlled in order to be included inside the viscous sublayer. According to Kirkgoz and Ardiclioglu (1997) cited in [33], the dimensionless distance must satisfy *<sup>y</sup>*<sup>+</sup> <sup>=</sup> *<sup>ρ</sup>u*<sup>∗</sup> *y <sup>μ</sup>* ≤ 10, where

*<sup>u</sup>*<sup>∗</sup> <sup>=</sup> *τ<sup>o</sup> <sup>ρ</sup>* and *τ<sup>o</sup>* is the boundary shear stress.

In order to increase stability of the simulation, strong under relaxation in the time scheme, and an advance for momentum equations is introduced, although with the risk of obtaining results with too many errors.

Diffusive fluxes are discretized in space using CDS (Central Differencing Scheme), while the contribution of CDS vs. UDS (Upwind Differencing Scheme) in convective fluxes is blended with a factor of 0.5.

The convergence of the iterations at each time step was checked when the residuals of each of the equations solved attained a value lower than 1 × <sup>10</sup>−6, while the maximum number of iterations for each time step was set to 5.

#### *2.6. Power Calculation*

The available (maximum) power is the theoretical hydraulic power of the flow through the HPM:

$$P\_h = \rho \text{gQ} \Delta H = \rho \text{gQ} (h\_1 - h\_2 + \frac{v\_1^2 - v\_2^2}{2g}) \tag{8}$$

where *Q* is the discharge, Δ*H* the head difference upstream and downstream of the HPM calculated from the water levels (*h*<sup>1</sup> and *h*2) and the free surface velocities (*v*<sup>1</sup> and *v*2). The dynamic head (last term in Equation (8)) can be neglected as the velocities are very low. Due to intrinsic losses during the energy transfer, the hydraulic power exerted on the HPM is lower than the available hydraulic power.

The total hydraulic moment (*M <sup>h</sup>*) acting over the HPM should be calculated from the simulation results, through the contribution of pressure and viscous tensions acting on every cell face of the rotating grid blocks corresponding to the wall type boundary condition (those adjacent to the wheel), multiplied by the distance between the center of the HPM and the center of the cell face (*d*-):

$$
\vec{M}\_h = \int\_S \vec{d} \times [-p n \mathbf{\hat{s}} + \mu (\nabla \vec{v}\_A + \nabla \vec{v}\_A^\dagger)] dS \tag{9}
$$

Because of the particular shape of the analysed HPM, the contribution of viscous tensions could be neglected, as its direction is parallel to the distance vector for almost every cell face. This was confirmed by the results of the simulations made by [15], even for the HPM with radially twisted blades.

Then, the hydraulic power (*Ph*) exerted on the HPM could be obtained by multiplying the total moment and the rotating speed (*ω*-):

$$P\_h = \vec{M}\_h \times \vec{\omega} \tag{10}$$

#### **3. Results and Discussion**

As the solver has not been used extensively for free surface flows or turbomachinery simulations, it was advisable to assess each approach individually. The SM method has been previously used and verified within caffa3d [20]. First, a simulation of a free surface flow over a rigid and stationary object using the VOF module was performed. Then, both modules were incorporated into the solver, and a simulation of the HPM performance was conducted with the configuration mentioned in Section 2. Different meshes and time steps were used for each part of the assessment process, since each approach has distinct simulation demands. The post-processing was made in Paraview.

#### *3.1. VOF Assessment*

For the assessment of the VOF module, it was performed a simulation of a free surface flow over a semicylindrical object of 0.196 m diameter lying down roughly on the middle of the bed of a 2.4 m long and 0.2 m width channel. The case was chosen equal to the one analysed in [33] by means of numerical simulations in 2D (in ANSYS-Fluent v.12.1 with RANS method and several turbulence closure models) and experimental tests, in order to validate the simulations run with caffa3d. The mesh used here (see Figure 7) was different than the one in Section 2.3 (mainly because the HPM is not part of the simulation). The maximum value of the normal dimensionless distance to the wall is *max*(*y*+) = 7.0.

**Figure 7.** Computational mesh around the semicylindrical dam.

Lateral sides of the channel were treated as non-slip walls, while the boundary conditions for the inlet, top and outlet were the same as in Section 2.4. Upstream of the object the flow is subcritical with a depth *ho* = 0.147 m (roughly 1.5 times the radius of the cylinder), and an inlet uniform velocity *vo* = 0.187 m/s (corresponding to a discharge *Q* = 0.055 m3/s as in [33]), resulting in a Froude number *Fr*<sup>1</sup> = *v*1/ *g* × *h*<sup>1</sup> = 0.156. At *t* = 0 the water was at rest occupying the entire channel up to a depth *h*<sup>1</sup> = 0.147 m.

Downstream of the object the flow becomes supercritical with *Fr*<sup>2</sup> = 2.26. This change from subcritical to supercritical also occurs during the operation of the HPM, thus this validation becomes more attractive.

After a few seconds, the flow becomes stable, so only the results corresponding to the last time (*t* = 20 s) are shown. In addition, although the simulation was done in 3D, the flow is quite bi-dimensional, and thus the results are shown only for the midsection of the channel. The comparison of the free surface profile over the semicylindrical dam is shown in Figure 8. The profile in Figure 8b was obtained through experiments and numerical simulations performed with the Reynolds Stress Model (RSM) as turbulence closure by [33]. Similar comparison for the computed streamlines over the dam are shown in Figure 9.

**Figure 8.** Comparison of free surface profile in the channel over the semicylindrical dam: (**a**) obtained with caffa3d and (**b**) presented in [33] through exprimental (o) and numerical (-).

These results show a quite well performance of the simulation, and are very similar to the results presented for similar studies [33,34].

The free surface level upstream of the dam obtained with caffa3d is a little higher than the one obtained by [33], where it remained in the initial value of *h*<sup>1</sup> = 0.147 m. This might be related with the fact that the streamlines in Figure 9a show a separation region bigger than in Figure 9b, thus raising the level of water. However, the free surface level and streamlines downstream of the dam are very similar in both figures, even with the separation region obtained with caffa3d being a little bigger. In [33], the authors also found differences in the size and shapes of the separation regions among the several turbulence models analysed, so the difference encountered here may be due to an error in the implementation of LES and/or the Smagorinsky model.

There is a small recirculation of air from the top surface to the outlet section above the free surface, but as this occurs just near the outlet section far away from the dam, it does not appear any effect of it in Figure 9a.

**Figure 9.** Comparison of streamlines in the channel over the semicylindrical dam: (**a**) obtained with caffa3d and (**b**) presented in [33].

The presence of spurious velocities in the air region has been previously reported by several authors to be caused by errors related to the calculation of the surface tension force [27,35], and several methods to solve it are also reported in [31,36–38]. However, when performing simulations without surface tension force, the recirculation persisted.

Another cause of spurious velocities could be the use of a predictor–corrector procedure [28] to reset values of *Vf* greater than 1 or lower than 0, which could also affect stability of the solution. Because non physical values of *Vf* are rare, the simulation without the corrector step was performed. Nevertheless, the recirculation did not disappear.

Afterwards, the recirculation was ignored, as it appears not to affect the flow near the object (this is where the HPM would be located).

Three horizontal velocity profiles at different sections of the channel, upstream of the dam (*x* = 0.8 m), over the dam (*x* = 1.0 m) and downstream of the dam (*x* = 1.3 m), are shown in Figure 10, enabling a direct comparison between the results with those presented in [33] for all the analysed turbulence models and the experimental test.

**Figure 10.** Horizontal velocity profiles obtained with caffa3d for sections: (**a**) *x* = 0.8 m; (**b**) *x* = 1.0 m; (**c**) *x* = 1.3 m, and by experiments and simulations [33] for sections: (**d**) *x* = 0.8 m; (**e**) *x* = 1.0 m; (**f**) *x* = 1.3 m.

One should compare the results obtained in caffa3d with the experimental and the numerical results obtained with the RSM, as it is the one that showed the best performance among the six analysed models in [33]. In that case, the shape of the profiles are similar, but the maximum velocity reached with caffa3d simulations in the three sections observed in Figure 10 are lower than those presented in [33].

#### *3.2. Simulation of a HPM*

The mesh, boundary conditions and general configuration of this simulation correspond to the ones explained in Section 2. For the initial condition of this simulation, the water level for the upstream zone of the channel corresponded to the height of the top of the HPM's hub, and to the height of the bottom of the HPM's hub for the downstream zone of the channel. All the fluid is at rest at the initial time. From the laboratory tests performed with the full width HPM, the actual flow rate corresponding to the best efficiency point at a rotational speed of 15 rpm was found to be *Qbep* = 11.5 × <sup>10</sup>−<sup>3</sup> <sup>m</sup>3/s. The wheel rotates at a constant speed of 15 rpm and there is an inlet flowrate *<sup>Q</sup>* = 5.75 × <sup>10</sup>−<sup>3</sup> m3/s (resulting from a uniform velocity *v*<sup>1</sup> = 0.13 m/s and a water depth *h*<sup>1</sup> = 0.31 m).

The predictor–corrector procedure was affecting the stability of the simulation, especially in the zone close to the inlet of the channel. This was suspected to happen according to [28], so only the predictor step was used. Despite disabling the corrector step, non-physical values of fraction of volume were not observed.

The maximum value of the normal dimensionless distance to the wall is *max*(*y*+) = 40 in the vertical plane, but reaches much bigger values for the *z* direction in band of blocks containing the lateral gap (up to 200). Even though this will probably affect the numerical stability, it was a compromise between this and the computational extra time of dealing with a finer mesh.

3D images of the water flow and the streamlines at time *t* = 4.0 s (after one exact full revolution of the wheel has been performed) in the surroundings of the HPM are shown in Figures 11 and 12, respectively, where the water flows from left to right of the image. In addition, Figure 13 shows the magnitude of the velocity in a plane normal to the flow at the section of the rotation axis (*x* = 0.0 m), in order to better assessed the quantity of flow through the lateral gap.

**Figure 11.** 3D image of water flow near the HPM, at time *t* = 4.0 s.

**Figure 12.** Streamlines of water flow near the HPM, at time *t* = 4.0 s.

**Figure 13.** Magnitude of the velocity at the section of the rotation axis (*x* = 0.0 m), at time *t* = 4.0 s.

The importance of reducing the width of this gap at the section of the rotation axis becomes clear observing the high velocities of water passing through the lateral gaps in Figures 12 and 13. In the present simulation, the channel was modeled with constant width along its entire length, meaning that the lateral gap was equal to the distance between the wheel and the channel's wall (see Section 2.3). This distance needed to be big enough to enable the water entering laterally (besides longitudinally) facilitating the filling process, as well as to enable the water leaving the buckets laterally (besides longitudinally) to facilitate the emptying process. It was observed during the study that, without modeling this distance (as in a 2D simulation), the filling and emptying process were greatly hindered and, as a result, a poor performance of the HPM was obtained, as also observed by [15]. In fact, it has been pointed out by [18] that an important gain in power and efficiency is obtained by increasing the ratio of the channel width to the wheel width up to a value of 1.3 (for bigger ratios, both power and efficiency would remain almost constant). This was the reason for choosing the width of the gap as 0.3 m. In summary, the filling and emptying of the buckets of the wheel were well performed (as can be seen in Figure 11), but a great amount of water passed through the lateral gap reducing the power absorbed by the HPM and affecting its performance.

The pressure acting over the surfaces of the blades and hub of the HPM can be analysed in Figure 14, where the water flows from the left to right of the image.

**Figure 14.** Manometric Pressure acting on the surface of the HPM, at *t* = 4.0 s.

Negative values of pressure can be observed at the edges of some blades, the edge of the hub and over the surface of the shaft, due to flow separation in those zones.

Slice views at the midsection of the channel at *t* = 4.0 s of fraction of volume, manometric pressure and 2D velocity vectors are shown in Figures 15–17, respectively. The good performance of the filling and emptying of the wheel is also observed in Figure 15. The reduction of the pressure of the water while passing through the wheel also has a good and expected behaviour. A small zone of water with negative pressure is observed in the bucket leaving the bottom shroud at the downstream side, but the value does not represent any risk of cavitation. Lastly, the velocity vectors look quite well in almost every zone, noting a reverse flow near the free surface just upstream the wheel that could be due to the impact of the blade entering the water.

**Figure 15.** Slice view at the midsection of the channel of fraction of volume, at *t* = 4.0 s.

**Figure 16.** Slice view at the midsection of the channel of manometric Pressure, at *t* = 4.0 s.

**Figure 17.** Velocity vectors of the slice view at the midsection of the channel, at *t* = 4.0 s.

In the following, a sequence of images for three instants of the rotation of the wheel, showing the fraction of volume and pressure fields around the HPM in a slice view at the midsection of the channel, is presented in Figure 18.

The three instants correspond to times *t* = 2 s (half revolution of the wheel), *t* = 3 s (three quarters of revolution of the wheel) and *t* = 4 s (full revolution of the wheel). The wheel in the images appears not to be moving, but it is just a visual effect because between each of the three instants the wheel rotates exactly a quarter of revolution.

**Figure 18.** Sequence of fraction of volume (**left**) and pressure (**right**) around the HPM at times *t* = 2 s (**up**), 3 s (**middle**) and 4 s (**down**).

From the evolution of the free surface until the first full revolution of the wheel, it seems that the channel is emptying downstream of the wheel. This could be due to the zero manometric pressure acting on the outlet section of the channel, and perhaps it would be better to try a hydrostatic pressure profile or install a little dam to generate the water level. The evolution of the pressure field follows the behaviour of the evolution of the fraction of volume field.

Regarding the hydraulic power exerted on the HPM, at time *t* = 5.0 s, it reached a value of *Ph* = 1.0 W, and before that time it presented negative values. It is thought that this low value could be associated with the high leakage flow through the lateral gap. From previously conducted experiments ([19]), the shaft mechanical power generated by the HPM for the analysed discharge was *Pexp* = 5.4 W (corresponding to the best efficiency point of this machine). Then, if we would consider the mechanical losses, it could be concluded that the power exerted on the HPM obtained through the numerical simulations are of the same order of magnitude of the experimental value (if such would have been calculated).

Finally, it is worth noting that a mesh independence analysis was not carried out during the present study, but it would be interesting to do it, as it was found that a steeper refinement near the walls would be needed to improve the numerical accuracy.

#### **4. Conclusions**

The recent academic research on HPM has been reviewed and, despite all the advances in its knowledge, a lack of open source solvers usage in the numerical simulations was recognised. Then, the solver caffa3d, developed within an academic environment in Uruguay, was chosen for performing the simulations. This is a well proven solver; however, it had not been used until now for free surface flows or for turbomachinery simulations.

The Sliding Mesh approach for solving the rotation of the HPM was chosen among other methods available in the solver. It has been previously validated, but yet it required some adjustments because it was outdated.

For the solution of the free surface flow, the VOF method is implemented in the solver. Some issues regarding the use of VOF still exist, such as spurious velocities in the air zone, recirculation of air from the top surface to the outlet section and instabilities near the inlet section of the channel caused by the predictor–corrector procedure. It was tested through a simulation of a free surface flow over a semicylinder lying on the bottom of the channel, and the results were similar to the results obtained by other authors.

Preliminary quite acceptable results from the merging of both modules (SM and VOF) were obtained after a simulation of the performance of a small scale HPM with straight blades. Thanks to the significant space between the wheel and the lateral walls of the channel, the filling and emptying process of the buckets of the HPM were adequately performed. However, in contrast, a great amount of leakage flow passed through the lateral gap, hindering the performance of the HPM. The boundary condition applied to the outlet section of the channel should be changed to a hydrostatic profile (of the desired depth), or a small dam should be positioned, in order to assure the downstream level of the free surface, avoiding the emptying of the channel. The value of the power exerted on the HPM assessed by the pressure field obtained through the numerical simulations resulted of the same order of magnitude (and very close) of the shaft mechanical power obtained through previous experimental tests.

The next steps in this line of research will be:


Once these aspects are covered, modifications of the geometry of the HPM pursuing better performance could be assessed with this open source flow solver.

**Author Contributions:** Conceptualization, H.M.R.; methodology, R.P.; software, R.P. and G.U.; validation, R.P.; formal analysis, R.P., G.U. and H.M.R.; writing–original draft preparation, R.P.; writing–review and editing, R.P., G.U. and H.M.R.; supervision, G.U. and H.M.R.; project administration, H.M.R.; funding acquisition, R.P. and H.M.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors wish to thank to the project REDAWN (Reducing Energy Dependency in Atlantic Area Water Networks) EAPA 198/2016 from INTERREG ATLANTIC AREA PROGRAMME 2014–2020 and the programme Young Professors and Researchers of Santander Universidades, for the corresponding grants.

**Acknowledgments:** The simulations reported in this article were performed in ClusterUY, a newly installed platform for high performance scientific computing at the National Supercomputing Center, Uruguay.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
