Variable Speed Control in PATs: Theoretical, Experimental and Numerical Modelling

1. Introduction

The availability of water resources at the global level has significantly decreased. Among the main agents that have caused this situation are climate change, environmental pollution, human activities, and failures in hydraulic structures, among others. It is becoming increasingly complicated to access appropriate sources that meet the quality and quantity of the resource. Despite this, water loss due to leaks in pressurized distribution systems still manages considerable values, with losses of 8% to 24% in developed countries [2]. If it is considered that the need for water is increasing, there is an urgent need to implement sustainable projects that allow the user to be carried out efficiently[3]. This type of project requires the use and development of new technologies that are easy to implement and apply[4].

In the case of sustainable water systems, there are some approaches from which improvements can be proposed. Among these approaches are the determination of water quality parameters [5], optimization of energy efficiency [6,7]reduction of water leaks [8,9], mathematical modelling of management, and optimization of systems [10,11], among others.

Water distribution systems are not energy efficient because they depend on pressure demands that can generate leaks, increasing energy costs [12]. One of the elements that have a negative effect from the point of view of energy efficiency, but is necessary for the hydraulic operation of the systems is the pressure-reducing valve (PRV) [13]. PRVs are used to reduce the pressure at one point by regulating the flow passage. An alternative to the use of these devices, to reduce dependence on non-renewable energy[14]and take advantage of the excess energy of these systems [15], is the use of PATs (Pumps working as a turbine). In addition, PATs have been used as energy-generating devices in micro-hydroelectric power plants as a sustainable solution in the water industry[16]. For this reason, it is a trend to study the use of PATs to optimize different water systems to improve their sustainability [17,18,19].

PATs are pumps that work in reverse mode to generate energy. This machine's cost is cheaper than a conventional turbine of the same size[20], although they have lower hydraulic efficiencies in ranges between 0.6 and 0.7 [21]. When all electromechanical equipment is considered, the overall efficiency decreases to values between 0.5 to 0.6 [1]. The use of pumps operating as turbines (PATs) increases due to their application, availability, and cost advantages [7,22,23,24,25,26,27]. For example, Novara et al. [28]conclude that an installation with PATs could be 5 to 15 times cheaper than a conventional installation with turbines.

The PATs study begins with Thoma and Kittredge [29], who accidentally found that pumps can operate efficiently as turbines when trying to evaluate the complete characteristics of the pumps. In 1957, Stepanoff [30] reported several modes of operation of the pumps on performance curves plotted in quadrants. Once it was discovered that PATs could be applied in the chemical industry and the supply of drinking water, different researchers developed some techniques to predict the operation of this type of machine. In 1962 Childs [31] carried out comparative studies between efficiencies in devices working in both modes (pump-turbine). Subsequently, the first studies were carried out to predict the performance values in turbine mode and get the Best Efficient Point (BEP) through linear equations. The PATs study has been developed with different approaches, such as in water distribution systems, where Jain's [32] research stands out the results of placing PATs in distribution systems. Fecarotta [33] y Morani [12] proposes an analysis regarding the proper location of PATs; the latter focuses its research looking for cost reduction and maximization of production and energy savings. Moazeni [34] investigates finding the optimal number and location of PATs through mixed nonlinear programming models. Macias [19] established a methodology that was applied in an irrigation project in a rural area in the province of Valencia (Spain) that focuses on optimizing the location and selection of PATs based on the influence of leaks. The same author [14] develops a new methodology for self-calibration of leaks to know the injected flow rate and the volume consumed in the water networks. This methodology was applied in the city of Manta, Ecuador.

Since the performance curves are not available in pumps that work in turbine mode [27], different studies and methodologies have been carried out to obtain them, and to select the appropriate machine depending on the type of working condition required. Rossi [35] proposes a general method to predict PATs performance using artificial neural networks (ANN). Based on the datasheets provided by the pump manufacturers, the author obtains the BEP and off-design performance using the ANN methodology. In addition, the resulting predictions were compared with experimental data not used in the training process, which resulted in a high degree of compatibility. The study concluded that the BEP flow rate increases in reverse mode while the specific speed in BEP decreases slightly and also recommends the use of this tool to choose the proper PAT. To estimate the BEP and the characteristic curves of the PATs, Perez-Sánchez [21] proposed new approach equations from an experimental base of 181 machines. In the same sense, Plua [1] presents new empirical expressions to estimate head, efficiency, and power curves for PATs with variable speed. These equations allow the application of various operation strategies in hydraulic simulation tools (e.g., Epanet, WaterGEMS).

Micro hydroelectric power plants (MHP) have become very effective solutions for rural sectors with powers of 5-100 kW. The big problem with these facilities is the high turbine cost concerning the entire project [27]. In the case of MHP, the price of these elements can be higher than 60%–70% [36]. One possibility to reduce this cost is to use PATs instead of a conventional turbine [37], which would favour the expansion of MHP and the reduction of greenhouse gas emissions [28]. In 2012 Pascoa [38]proposed a new approach for a hydroelectric plant with PAT with constant flow. Rossi [39]suggested the economic feasibility of placing PATs in the Merano aqueduct, which resulted in the production of 338 kWh of daily electricity and power of 19.18 kW. Table 1 [27]shows different PATs installations in power generation projects.

CFD techniques have been widely used to predict characteristic curves and performance of pumps in direct and reverse modes and proved to be an effective solution in PATs approaches [40,41,42,43]. Also, an experimental investigation is fundamental for obtaining reliable results in PATs under different optimization stages [1,15,44,45]. Different types of machines such as axial, mixed, and radial PATs; with horizontal and vertical axis; single and multistage [46], have been studied using CFD simulations for fixed and variable speeds [24,25,26]. However, very few studies related to numerical modelling in PATs of variable rotational speed have been executed, so it is imperative to establish equations and laws that predict their behaviour [47]. The numerical simulations were carried out to define the performance of the pump [48], analyze flow in turbine mode [49], predicting and extrapolate characteristic curves [50], among others.

Plua [40] presented research in which the main parameters and techniques that have been simulated for PATs through CFD are shown and whose main simulation ranges are mentioned: specific speed: 0.8 -306; rotational speed: 250-3900 rpm; flow rate: 2.9-300 l/s; mass flow: 13-17.8 kg/s. Concerning numerical simulations, the principal turbulence models used were Reynolds Average Navier Stokes (RANS) and Unsteady Reynolds Average Navier Stokes(URANS). The most used closure model was k-e, followed by k-w and k-w-SST, among others. Regarding packages, ANSYS-CFX was the most used, followed by CFD Code Fluent and OpenFOAM. Respecting mesh generation, the number of cells was 1x10⁶ to 4.2x10⁶, with hexahedral, tetrahedral, mixed blocks, and pyramids. Depending on each situation, boundary conditions such as total pressure, mass flow rate, stagnation pressure, constant total pressure, static pressure, and volumetric flow were placed at the inlet and outlet of the model. In conclusion, it was established that the CFD methodology to predict the performance of a pump working as a turbine presented adequate accuracy based on the comparison of results with the experimental tests. However, numerous errors were also reported in some studies. The authors assumed that the reported errors are due to the geometries between the tests and the simulations not being identical, the loss estimation is not exact, and more experience in computational analysis is required for modelling this type of phenomenon. Finally, the same author [47] evaluated the application of numerical CFD simulation in PATs in comparison with experimental results and obtained conclusions for future numerical analysis. As a result, it was evidenced that there are few simulated cases where flow with variable speed is simulated and that the number of studies with free code computational packages is minimal, and that its use should be promoted due to its outstanding capabilities.

Therefore, the present study is focused on a numerical simulation in the OpenFOAM 3D free code package of PATs that have experimental data to validate the use of the new empirical expressions proposed for machines with different rotational speeds. The particularity in the modelling is that the study of a rotating PAT at different speeds will be carried out, and comparisons will be made with experimental results obtained on a test bench to calibrate the model.

2. Materials and Methods

Figure 1 depicts the main tasks performed to determine the validity of new expressions obtained by Plua [1] to predict the behaviour of PATs with variable speed:

2.1. Preprocess

2.1.1. Computational Domain

The PAT model presented in this study was taken from research conducted by Pérez Sánchez [51] and experimentally tested at the CERIS-Hydraulic Lab of Instituto Superior of Lisbon. Geometry corresponds to an installation of a PAT in a laboratory that allows experiments where flow, pressure conditions, and rotational speed can be varied. The hydraulic facility consists of a 1m³ air-vessel tank, a 50 mm HDPE pipe, a KSB radial impeller centrifugal pump (model Etarnom 232) that operates in turbine mode, regulating tank, pressure transducers, valves, and a flow recirculation pump. The air-vessel tank sends water to reach the PAT, which finally discharges to the open free surface tank and then incorporates it into the system through the recirculation pump. The 3D model was built in SolidWorks CAD system from which the following drawing view is extracted, see Figure 2.

The computational domain consists of four parts: Inlet pipe starting in the inlet section and reaching the pump's runner, rotating part that is the impeller of PAT, rotating part of the domain; Casing, stationary part of the pump and; Outlet pipe, which corresponds to the discharge of the pump to the outlet. The original geometry was redefined according to the configuration of the control volumes to obtain optimal meshing. See Figure 3.

Depending on the actual geometry and its characteristics, the .stl files were modified with the Autodesk Inventor Software to get a better-quality mesh. The areas of meshing interest were prioritized: the casing, the impeller, and the blades. Each has different elements and details simultaneously with different levels and definition angles. It allows the surfaces to stick more to the edges, bringing the mesh's geometry closer to the actual configuration. The geometries modelled were the volute, the discharge pipe, the inlet pipe, and the impeller. In the case of the impeller, it was divided into three parts, as seen in Figure 3. The impeller is composed of the lower and upper parts and the blades. These elements are treated independently to improve the mesh quality and then facilitate the visualization of results at the post-processing stage. In addition, six blades were configured inside the impeller, which allows a better study of the phenomenon presented in the PAT.

2.1.2. Mesh

The mesh was created with snappy Hex Mesh, an automatic mesh generator that adjusts to the surface to obtain the required mesh. First, the 3D model was exported to the format.stl through Autodesk® Inventor® software. Later, with the help of HELIX-OS, the BlockMeshDict file was created to generate through BlockMesh utility orthogonal mesh elements for the Casing, Inlet pipe, Impeller, and Outlet, respectively. Once the block meshing was ready, the domain geometries were admitted in the snappyHexMeshDict file. The local refinement was defined through castellatedMesh, and the internal points within the closed domain were entered. Finally, it was necessary to use the topoSet tool to generate zones with movable cells for the runner and merge the meshes with the mergeMeshes utility. Mesh characteristics are presented in Table 2, and generated mesh is in Figure 4.

2.1.3. Approach

MRF Technique (Multiple Reference Frame) is the technique used for modelling rotation in CFD in this case. This methodology establishes a separate reference frame for each region of the domain, for both rotational and static [49]. It is based on the creation of a local region around the rotating object where the relative velocity is determined for each point. First, Navier Stokes equations are built, taking account of centrifugal and Coriolis forces, and then, a set of equations for stationary and rotational regions are created. This technique can accurately capture instantaneous local flow, which depends on the relative position of the rotative element versus static geometry. In the MRF approach, the Navier-Stokes equations are solved in terms of the global/inertial velocity. Since, in this case, there is a separation between the impeller and the scroll, the AMI approach is not applicable. For that reason, a set was used that allows simulation these elements in the MRF approach.

2.1.4. Boundary and Initial Conditions

Initial and boundary conditions should be applied when solving the Navier-Stokes and continuity equations. Table 3 summarizes the initial conditions related to the turbulence models used of this research. For the calibration of the mathematical model, the k-e turbulence model was used (the same one used by [51]). k is turbulent kinetic energy, and e is turbulent dissipation rate. The k-w(specific turbulent dissipation rate)-SST turbulence model was used to analyze the experimental data, the nominal rotational speed curve, and the results of the new expressions contained in [1].

The turbulent kinematic viscosity value “nut” represents the roughness in the walls confirming the domain. Regarding boundary conditions, a constant velocity input condition and a static pressure output condition were used. The boundary conditions of the computational domain are detailed in Table 4.

2.2. Numerical Simulation

2.2.1. CFD

The Navier-Stokes equations were solved using CFD methods based on a continuum mechanics approach for fluid mechanics to define the fluid behaviour in the PATs [50]. For that, two equations were considered that obtain the values of velocity and pressure that allow for defining the average behaviour of the flows. The equations correspond to the conservation of mass and linear momentum and are indicated in a tensor with the following expressions:[47]:

(1)

(2)

where i and j are subscripts for the three axes of space respectively, $\overline{u}$ is the filtered velocity magnitude, $\overline{P}$ is the filtered pressure, the subgrid stress tensor is ${\overline{T}}^{ij}$, and ${\overline{T}}_{ij}^{\prime }$ the filtered viscous stress tensor.

2.2.2. CFD and Solvers

The CFD package used is the CFD OpenFOAM 9 which models multiphysics simulations applicable to Computational Fluid Dynamics for incompressible and compressible flows with application in dynamic mesh management to make rotating reference frames with adaptable mesh refinements as required. OpenFOAM uses a directory structure to solve the cases, where the case is the name of the analysis case; the system sets the numerical control to run time and solver; the constant contains the physical properties, modelling and mesh information; 0 has the edge conditions as well as the beginning to modelling and, time-directories that correspond the solutions and derived cases.

Regarding meshing, OpenFOAM has some mesh utilities such as BlockMesh, SnappyHexMesh, foamyHexMesh, and foamyQuadmesh. OpenFOAM also allows the mesh to be generated with other packages since mesh conversion utilities are compatible with popular mesh formats (Gmsh, Fluent, Ideas, and Netgen, among others). As stated above, SnappyHexMesh generated the mesh. SnappyHexMesh utility is an automatic hybrid mesh that divides, refines, and adjusts to the analyzed surface, attaching the mesh with complex details of geometry [52]

For the calibration of the model, the solver simpleFoam was applied to steady-state incompressible flow based on the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm for pressure velocity-coupling [54] with applications in turbulent and transient flows in pipes.

3. Results

3.1. Numerical simulation validation

3.1.1. Mesh quality

The checkMesh tool was used to evaluate the mesh quality, giving the mesh stats, the overall number of cells of each type, topology, geometry, and conclusions about the mesh. Two parameters are used to verify the quality of the mesh; one of them is Ω, which corresponds to the following expression Ω =NE/ND, where ND is the number of nodes and NE is the number of elements.Ω indicating the homogeneity of the mesh, a good mesh quality will present Ω values close to 1, and values close to 2 have very dispersed meshes. For this case study, the calculated Ω value is 0.69, which is acceptable. The other value is the so-called y⁺ which verifies the acceptable range of values for the turbulence model. If this value is less than 1, it is considered that the quality of the mesh is good. In this study, it was found that the average y⁺ values in all the simulations of the mesh are less than 1.

3.1.2. Calibration

For the CFD simulation validation, two calibrations were performed concerning the Pérez-Sánchez study [51]. The first concerns the mathematical model made with SolidWorks-FloEFD and the second concerns experimental research. In the Pérez-Sánchez CFD model, the simulated global variables were: the head (H), the output hydraulic torque (T), the discharge (Q), and the rotational speed (N). Within the simulations, the absolute static pressure contours were obtained for a flow rate of 4.5 l/s and rotation speeds of 810rpm, 930rpm, 1050rpm, 1170rpm, 1275rpm, and 1500rpm. The results show that the pressure decreases from upstream to downstream as the fluid flows within the domains and along the impeller, from the inner to the outer region, as the energy is transmitted to the shaft. On the other hand, it was found that the higher the speed, the lower the pressure value downstream of the impeller.

The results of the simulation performed with OpenFOAM in this study, are shown in Figure 5 and Figure 6. As can be seen, the pressure decreases from upstream to downstream, and the lowest pressure value occurs at point D (before the first elbow of the volute outlet) for maximum speed. The error of this simulation concerning the original work varies in ranges from 0.014% up to 14.297% at points A, B, C, and F of the model see Table 5.

Regarding the calibration of the mathematical simulation with the experimental data, a sensitivity analysis was performed to identify which turbulence model produces the best results. Simulations were executed on the machine's Best Efficient Point (BEP) tested in [51] when operating in turbine mode (Q_BEP=3.6 l/s) for speeds of 200, 600, 880, 1020, 1200 and 1500 rpm, using the k-e k-w-SST models. The results obtained for both simulations are shown in Table 6. As can be seen, the simulations produced errors of similar magnitude. Still, for the nominal rotational speed of 1020 rpm, the k-w-SST model is the one with the lowest error. An error index analysis was performed to define the turbulence model with which the cases of experimental data, nominal rotational speed curve, and the results of the new expressions [1] will be simulated. Figure 9 shows that in all cases, the k-w-SST model has better performance. So k-w-SST turbulence model was adopted for the rest of the cases.

Once the mathematical model has been validated concerning the results obtained in the numerical modelling and experimental works of [51], the curve of the machine working at nominal speed is contrasted, as seen in Figure 10.

Figure 8. The nominal curve obtained with CFD OpenFoam vs Nominal curve in [51].

Figure 8. The nominal curve obtained with CFD OpenFoam vs Nominal curve in [51].

As can be seen, the simulation performed with OpenFOAM presents satisfactory results, and therefore the model is considered validated.

3.2. Analytical Expressions Validation

3.2.1. Analytical Expressions-New expressions to predict PATs behaviour

Considering that in the case of PATs, the information to select the suitable machines is not known because it is not provided by the manufacturers [55], polynomial expressions have been proposed as a function of semi-empirical methods to estimate the characteristic curves in PATs when the rotational speed is constant [21,55,56,57]. However, considering that flow rates in water systems are variable due to user demand, an optimal energy analysis for PATs cannot be performed if the rotational speed is considered constant. Therefore, strategies have been proposed to maximize energy when the machine works at different rotational speeds, called variable operation strategy (VOS) (19). Plua et al. [1]propose new empirical expressions applying the VOS strategy in water systems for different rotational speeds of 15 different machines and analyzing 87 different curves with 56450 operating points.

Through a mathematical analysis of 10 general expressions (polynomial 6 and 4 potential) considering specific variables as the ratio of rotational speed α, and the ratio Q/QBEP, it was possible to adjust a polynomial function for experimental values of head and efficiency, and a potential function for power. These expressions are observed in equations (3) to (7) and presented the lowest error (30 to 50% compared to other models) in the respective analyses performed where the RMSE, MAD, MRD, and BIAS indices were calculated. Figure 6 of [1] shows a head and efficiency curve comparison between the proposed model, experimental data, and other models.

$$q=-0,1525\left(\alpha \frac{Q}{Q_{BEP}}\right)+0,1958{\left(\frac{Q}{Q_{BEP}}\right)}^2-0,0118\left(\frac{Q}{Q_{BEP}}\right)-0,6429{\alpha }^2+1,8489\alpha -0,2241$$

(3)

$$h=-0,31070\left(\alpha \frac{Q}{Q_{BEP}}\right)+0,3172{\left(\frac{Q}{Q_{BEP}}\right)}^2-0,0546\left(\frac{Q}{Q_{BEP}}\right)+0,242{\alpha }^2+1,1708\alpha -0,3426$$

(4)

$$e=0,8271\left(\alpha \frac{Q}{Q_{BEP}}\right)-0,3187{\left(\frac{Q}{Q_{BEP}}\right)}^2-0,1758\left(\frac{Q}{Q_{BEP}}\right)-1,035{\alpha }^2+1,1815\alpha +0,5019$$

(5)

$$p={\alpha }^{2,4762}$$

(6)

$$q={\alpha }^{0,7439}$$

(7)

where

$$q=\frac{Q_i}{Q_{BEP}};h=\frac{H_i}{H_{BEP}};e=\frac{{\eta }_i}{{\eta }_{BEP}};p=\frac{P_i}{P_{BEP}}=qhe$$

3.2.2. Analytical Expressions Validation

The Head value H for different rotational speeds at points close to the BEP was compared with the experimental Head obtained in [51], the expressions proposed in [1], and the mathematical model. See Figure 11. As can be seen, the relationship between Q and H is increasing in all cases. The predictions made in the numerical simulation with OpenFOAM and with the new expressions present values close to the experimental ones when the operation of the machine approaches the BEP.

3.2.3. Error Analysis

The error rates obtained in the predictions made as a function of the rotational speed are presented in Figure 12. As can be seen, as the conditions approach those of the BEP, the predictions reflect values closer to reality. The range of calculated absolute error for rotational speeds 880 rpm, 1020 rpm,1200 rpm, and 1500 rpm is between 5% to 24%; 2% to 17%, 0% to 12%, and 1% to 24%, respectively.

4. Conclusions

This research proposes to validate as a prediction methodology of Flow (Q) vs Head (H) curves of variable speed PATs, the numerical simulation with OpenFOAM 3D Free Code Package depending on its configuration and working conditions. It also proposes to validate the new expressions submitted by Plua [1]. It was demonstrated that the simulation presents adequate results since the mathematical model and the nominal curve of the Pérez-Sánchez [51] research were calibrated. Furthermore, based on experimental data from a PAT, the Q vs. H curves were calculated through the new expressions [1] as well as with the numerical simulation performed in Openfoam, presenting satisfactory results as the operation point of the work approaches the BEP, since the trend of the generated curves, the slope thereof and the error indices demonstrate acceptable values. However, when moving away from the BEP conditions, the error increases.

The computational domain was divided into four parts: Inlet pipe, rotating part, casing, and outlet pipe. The mesh was created with snappyHexMesh con elements type Hexahedra, Polyhedra, Prism in several 827578. The control conditions imposed on the model were flow and static pressure, and the k-e turbulence model and the k-w-SST model were used, showing that the latter has better results.