Pressure Drop and Energy Recovery with a New Centrifugal Micro-Turbine: Fundamentals and Application in a Real WDN

: Water distribution networks need to improve system efﬁciency. Hydropower is a clean and renewable energy that has been among the key solutions to environmental issues for many decades. As the turbine is the core of hydropower plants, high attention is paid to creating new design solutions and increasing the performance of turbines in order to enhance energy efﬁciency of leakage by pressure control. Hence, design and performance analysis of a new turbine is a crucial aspect for addressing the efﬁciency of its application. In this study, computational ﬂuid dynamics (CFD) modeling is coupled with experimental tests in order to investigate the optimal performance of a new centrifugal turbine. The behavior of the ﬂow through the turbine runner is assessed by means of velocity proﬁles and pressure contours at all components of the machine under different operating conditions. Finally, the turbine geometry is scaled to a real water distribution network and an optimization procedure is performed with the aim of investigating the optimal location of both the designed new centrifugal micro-turbines (CMT) and pressure reducing valves (PRV) in order to control the excess of pressure and produce energy at the same time.


Introduction
Hydropower is a well-advanced technology that convert the potential energy of water into electric energy [1][2][3] and contributes to about 50% of the all existing renewable energy sources [4,5]. It is reported that at least 16% of worldwide energy production results from the hydraulic power generation [6,7].
The turbine represents the core device of any hydropower plant, converting the hydraulic energy in mechanical energy of a rotating shaft [8], which is then converted into an electrical current by means of a generator. Turbines can operate under steady-state or transient operations in hydropower plants [9]: the angular frequency of the runner, as well as the discharge through the turbine, is constant in the former case and varies over time for transient operations [10]. Moreover, steady-states include both best efficiency point (BEP) conditions and off-design conditions, such as part load (PL), high load (HL), and full load (FL), whereas transient operations include start-stops and load variations [11,12]. The operation of a turbine is significantly complex and is characterized by high turbulence,

Experimental Tests
Modeling the geometry of a centrifugal turbine involves defining the impeller and volute components. In this study the geometric modeling of the draft tube has also been considered, as it is one of the crucial components of a turbine. Performance tests have been carried out with reference to a centrifugal turbine prototype ( Figure 1) in order to assess the hydraulic behavior profile under different operating conditions. The experimental setup consists of a free surface reservoir, on the top of which a turbine is installed with a low-return suction tube. The setup is also equipped with an air vessel to control sudden pressure impulses, a flow-control circulation pump, and all required hydraulic and electrical components, such as a permanent magnet generator (or SEIG) and energy dissipation for flow control.
The turbine prototype is characterized by an intake external diameter of 46 mm (40 mm inside) and a shaft diameter of 12 mm (Figure 2). The turbine is installed on a rigid base and connected to a generator. A datum torque meter (model M425-S2-B) is installed between the turbine shaft and the electric brake.

Experimental Tests
Modeling the geometry of a centrifugal turbine involves defining the impeller and volute components. In this study the geometric modeling of the draft tube has also been considered, as it is one of the crucial components of a turbine. Performance tests have been carried out with reference to a centrifugal turbine prototype ( Figure 1) in order to assess the hydraulic behavior profile under different operating conditions.
The experimental setup consists of a free surface reservoir, on the top of which a turbine is installed with a low-return suction tube. The setup is also equipped with an air vessel to control sudden pressure impulses, a flow-control circulation pump, and all required hydraulic and electrical components, such as a permanent magnet generator (or SEIG) and energy dissipation for flow control. The turbine prototype is characterized by an intake external diameter of 46 mm (40 mm inside) and a shaft diameter of 12 mm (Figure 2). The turbine is installed on a rigid base and connected to a generator. A datum torque meter (model M425-S2-B) is installed between the turbine shaft and the electric brake.  In the first stage of the tests, the turbine is experimented with variable pressure, flow and variable speed, keeping the opening guide vane at maximum to establish the performance points. The measured values are flow, pressure, torque and rotational speed. During the tests, the flow was varied within the range of 19-46 L/s, whereas the head drop was kept as 32.8 m (3.17 bar). The measured torque range was 25-50 N·m., whereas the rotational speed varied between 750 and 1850 rpm. Finally, the power on the turbine shaft was between 2 and 10 kW. The experimental points are shown in Table 1.

Turbine Geometry and Boundary Conditions
Providing a qualitative and quantitative prediction of fluid, CFD tools represent an important technological advance toward a detailed analysis of the flow behavior. CFD modeling allows one to study the behavior of both turbulent and laminar flows, as well as exchanges of energy, flow phases, vorticity and turbulence levels. Based on the experimental data provided as boundary conditions, as well as the virtual prototype of the turbine (Figure 1), a turbine modeling was performed with the aid of a CFD software for turbulent 3D flow and a design/drafting software. In particular, the software COMSOL Multiphysics [37] was selected for the CFD modeling in terms of geometry modification, mesh generation, calculation and results analysis, presenting accurate results for several fluid flow problems. COMSOL is a finite element method (FEM) software, relying on the mass conservation and the RANS (Reynolds-averaged Navier-Stokes) equations. FEM is a computational method according to which an object is divided into smaller elements. Each element is assigned to a set of characteristic equations, which are then solved as a set of simultaneous equations to investigate the behavior of the whole object [38].The design of each machine component was, instead, performed by the software AutoCAD [39] and then imported by COMSOL to reproduce the physical problem and post-processing capabilities. The CFD approach was mainly intended to investigate the complex internal flows through all components of the machine at different operating conditions, which may be crucial to improve the performance of the machine at off-design conditions. The simulations were performed by using the k-ε model for turbulent flows (with Reynolds number ranging between 20,850 and 209,000), assuming water as an incompressible fluid and the fluid regime as steady state. This model exhibits a good convergence rate and relatively low memory requirements. A crucial assumption on which the k-ε model is based is that the Reynolds number is high enough. Moreover, such a model assumes that the turbulence is in equilibrium in boundary layers, which means that the production equals the dissipation. Being not always verified, such assumptions limit the accuracy of the model. Furthermore, in case of rotating flows, such a model often shows poor agreement with experimental data [40]. However, the flow close to the surface of the blade is turbulent; thus, the assumptions used to derive the k-ε model are not valid there. The correction of the k-ε model to better describe the flow in regions near the blades is not always desirable due to the significantly high resolution requirements. Analytical expressions, well known as wall functions, are instead preferred to describe the flow in the critical regions [41] and are used in this application. Wall functions ignore the flow field in the buffer region and analytically compute a nonzero fluid velocity [40][41][42]. When wall functions are used, the wall lift-off in viscous units δ + w has to be thoroughly verified to be below 30 in the whole region. In case of δ + w greater than the threshold value, a finer boundary layer mesh is added. All calculations were performed on an Intel 5, CPU 3.90 GHz, RAM 8 GB with four cores and threads running in parallel. The model uses a physics-controlled mesh, according to which the size attributes and operation sequences are automatically defined; thus, the mesh is automatically generated and adapted for the model's physics settings. Hence, once the geometric model was defined, the automatic mesh generation was performed by a careful selection of the related parameters, in order to achieve good precision results with an acceptable computational effort. The mesh near the surface of the blade was refined by an additional fine mesh to solve turbulence equations. Based on the assessment of the accuracy of the calculated maximum pressure, a mesh convergence analysis was tested by the use of a finer mesh and, thus, a greater number of degrees of freedom (DOF). In particular, local sizing and curvature size functions were used to obtain a sufficiently fine mesh to properly capture regions characterized by a rapid change in key variables, such as pressure and flow velocity. According to the mesh convergence analysis, for a mesh of approximately 150,000 DOF, the computed pressure value increased from 0.1 to 2.0 bar for a mesh of approximately 1,000,000 DOF, showing that 500,000 DOF provides the same precision as 1,000,000 DOF. Finally, the fluid domain was constructed with 280,056 tetrahedral components, with 85% of quality element. This number of elements is a compromise between the mesh quality as computed by the automatic mesh generator and the computational capacity. Nevertheless, as demonstrated in [43], the number of elements and their size ensure for limited error in the calculation of global parameters. The good agreement between calculation and experimental results in terms of global parameters, as shown below in Section 3.2, demonstrates the accuracy of the numerical solution.
The analyzed domain consists of an inlet, an outlet and two periodic interfaces. The boundary conditions at the inlet and outlet parts of the turbine were set according to the experimental measurements. In particular, the pressure value recorded by the pressure transducer located downstream of the draft tube was assumed as boundary condition at the outlet, whereas the flow discharge was considered as the the boundary condition at the inlet.
Hence, for each simulation, the rotation of the impeller was assumed as the rotational speeds achieved experimentally, and such a rotating domain was considered as a moving mesh. This definition regulates the spatial frame and applies to all interfaces of physics where domains rotate. Thus, in a rotating coordinate system in the inner domain, and in the outer one in fixed coordinates, these interfaces formulate the Navier-Stokes equations for incompressible fluids which are expressed in Equation (1) [44]: where µ is the dynamic viscosity, ρ the fluid density, p the pressure and F is the external forces applied to the fluid. Non-rotating parts can be considered as fixed parts, which were coupled with the rotating parts by means of an identity pair set, that is, the fluid continuity boundary condition.

Numerical Results
The design of a turbine impeller is generally based on given values of differential head, flow rate and rotating speed. This means that the specific speed is practically imposed. Hence, the novel geometry was designed to target the specific speed resulting from the experimental results. Thus, the turbine was simulated for different rotational speeds (varying between 700 and 1800 rpm), different flow values (ranging between 19 and 50 L/s), and keeping the head drop value at 32.8 m. Figure 3 shows the total pressure contours across the runner of the turbine for different rotational speeds. head, flow rate and rotating speed. This means that the specific speed is practically imposed. Hence, the novel geometry was designed to target the specific speed resulting from the experimental results. Thus, the turbine was simulated for different rotational speeds (varying between 700 and 1800 rpm), different flow values (ranging between 19 and 50 L/s), and keeping the head drop value at 32.8 m. Figure 3 shows the total pressure contours across the runner of the turbine for different rotational speeds.  Figure 4 shows the total pressure drop at the inlet and outlet rotor for different rotational speeds. According to this figure, such a difference in terms of pressure drop is quite small for most rotational speeds. Instead, a more significant difference is detected for a rotational speed equal to 1850 rpm.   Figure 4 shows the total pressure drop at the inlet and outlet rotor for different rotational speeds. According to this figure, such a difference in terms of pressure drop is quite small for most rotational speeds. Instead, a more significant difference is detected for a rotational speed equal to 1850 rpm.
posed. Hence, the novel geometry was designed to target the specific speed resulting from the experimental results. Thus, the turbine was simulated for different rotational speeds (varying between 700 and 1800 rpm), different flow values (ranging between 19 and 50 L/s), and keeping the head drop value at 32.8 m. Figure 3 shows the total pressure contours across the runner of the turbine for different rotational speeds.  Figure 4 shows the total pressure drop at the inlet and outlet rotor for different rotational speeds. According to this figure, such a difference in terms of pressure drop is quite small for most rotational speeds. Instead, a more significant difference is detected for a rotational speed equal to 1850 rpm.  A comparison of the flow velocity for different rotational speeds is presented in Figure 5. It is clearly shown that the geometry is more resistive to the flowrate for N < 1500 rpm for a given range of the total pressure drop. This results in a higher change in absolute kinetic energy over the range of pressure differentials, causing a higher change in dynamic pressure. rpm for a given range of the total pressure drop. This results in a higher change in absolute kinetic energy over the range of pressure differentials, causing a higher change in dynamic pressure.  Figure 6 represents for each rotational speed the velocity isosurface at the first instant the turbine is crossed by flow and starts to rotate. According to Figure 6, the turbine response becomes faster by increasing the rotational speed. The fastest response is attained at a rotational speed at 1850 rpm.  Figure 7 shows the velocity field at six operating points in a section plan of the turbine runner. The characteristics of the flow were analyzed in critical regions, such as near the hub and tip zones. With reference to Figure 7a, high flow acceleration is present from the leading to the trailing edge, but at 68% along the blade chord, a sudden increase in static pressure is observed, causing flow deceleration. Figure 7b,c shows almost the same conditions as in Figure 7a, with a sharper velocity magnitude near the trailing edges. Figure  7c-f shows some differences in terms of flow velocities, with higher values in the regions near the hub. The degree of flow deceleration is lower when compared with N < 1500 rpm, inducing lower internal losses as well.  Figure 6 represents for each rotational speed the velocity isosurface at the first instant the turbine is crossed by flow and starts to rotate. According to Figure 6, the turbine response becomes faster by increasing the rotational speed. The fastest response is attained at a rotational speed at 1850 rpm.
rpm for a given range of the total pressure drop. This results in a higher change in absolute kinetic energy over the range of pressure differentials, causing a higher change in dynamic pressure.  Figure 6 represents for each rotational speed the velocity isosurface at the first instant the turbine is crossed by flow and starts to rotate. According to Figure 6, the turbine response becomes faster by increasing the rotational speed. The fastest response is attained at a rotational speed at 1850 rpm.  Figure 7 shows the velocity field at six operating points in a section plan of the turbine runner. The characteristics of the flow were analyzed in critical regions, such as near the hub and tip zones. With reference to Figure 7a, high flow acceleration is present from the leading to the trailing edge, but at 68% along the blade chord, a sudden increase in static pressure is observed, causing flow deceleration. Figure 7b,c shows almost the same conditions as in Figure 7a, with a sharper velocity magnitude near the trailing edges. Figure  7c-f shows some differences in terms of flow velocities, with higher values in the regions near the hub. The degree of flow deceleration is lower when compared with N < 1500 rpm, inducing lower internal losses as well.  Figure 7 shows the velocity field at six operating points in a section plan of the turbine runner. The characteristics of the flow were analyzed in critical regions, such as near the hub and tip zones. With reference to Figure 7a, high flow acceleration is present from the leading to the trailing edge, but at 68% along the blade chord, a sudden increase in static pressure is observed, causing flow deceleration. Figure 7b,c shows almost the same conditions as in Figure 7a, with a sharper velocity magnitude near the trailing edges. Figure 7c-f shows some differences in terms of flow velocities, with higher values in the regions near the hub. The degree of flow deceleration is lower when compared with N < 1500 rpm, inducing lower internal losses as well.   Figure 8 presents the velocity streamlines between the runner and the guide vane. According to the streamline trend, it is clear that the fluid moves radially outward in the rotor boundary layer. With reference to Figure 9, the velocity contour shows that the maximum velocity reaches about 12 m/s in the region close to the leading edge of the rotating blade, whereas the minimum velocity occurs close to the right side of the vane.   Figure 8 presents the velocity streamlines between the runner and the guide vane. According to the streamline trend, it is clear that the fluid moves radially outward in the rotor boundary layer. With reference to Figure 9, the velocity contour shows that the maximum velocity reaches about 12 m/s in the region close to the leading edge of the rotating blade, whereas the minimum velocity occurs close to the right side of the vane.   Figure 8 presents the velocity streamlines between the runner and the guide vane. According to the streamline trend, it is clear that the fluid moves radially outward in the rotor boundary layer. With reference to Figure 9, the velocity contour shows that the maximum velocity reaches about 12 m/s in the region close to the leading edge of the rotating blade, whereas the minimum velocity occurs close to the right side of the vane.  Dimensional analysis is a useful and powerful tool to analyze complex problems and solve different technical/physical problems. Thus, non-dimensional hydraulic feature terms of the machine were described in accordance with the geometric and kinetic resemblance principles: Non-dimensional quantities defined by Equations (2)-(4) represent the head coefficient, the flow coefficient and the torque coefficient of the centrifugal turbine, respectively. Figure 10 shows the comparison between experimental measurements and numerical results. Scale effects, due to design parameters not perfectly adherent to real parameters, lead to a deviation between the prototype and the numerical model. In addition, since the Reynolds number of the model is typically higher than the prototype, losses of friction resulting from the CFD analysis are larger than the prototype. In addition, the model efficiency is generally greater because of the viscous effects associated with vortex or eddy effects. The maximum efficiency of the tested is attained as 65% for a value of the flow coefficient as about 0.45. Dimensional analysis is a useful and powerful tool to analyze complex problems and solve different technical/physical problems. Thus, non-dimensional hydraulic feature terms of the machine were described in accordance with the geometric and kinetic resemblance principles: Non-dimensional quantities defined by Equations (2)-(4) represent the head coefficient, the flow coefficient and the torque coefficient of the centrifugal turbine, respectively. Figure 10 shows the comparison between experimental measurements and numerical results. Scale effects, due to design parameters not perfectly adherent to real parameters, lead to a deviation between the prototype and the numerical model. In addition, since the Reynolds number of the model is typically higher than the prototype, losses of friction resulting from the CFD analysis are larger than the prototype. In addition, the model efficiency is generally greater because of the viscous effects associated with vortex or eddy effects. The maximum efficiency of the tested is attained as 65% for a value of the flow coefficient as about 0.45.
According to Figure 10, the comparison between CFD and experimental tests predicted that CFD results are slightly higher than the experimental results.
Based on experimental results and affinity laws, the following relationships were obtained: where q T and n T are experimental coefficients equal to 0.141 and 22.76, respectively, whereas α is the runner thickness. According to Figure 10, the comparison between CFD and experimental tests predicted that CFD results are slightly higher than the experimental results.
Based on experimental results and affinity laws, the following relationships were obtained: where q T and n T are experimental coefficients equal to 0.141 and 22.76, respectively, whereas α is the runner thickness.

Application to the Water Distribution Network (WDN) of Funchal City-Madeira, Portugal
Once the design and performance of the novel turbine were investigated by experiments and CFD analyses, an application to a real case study network was then considered. The case study concerns the water distribution network supplying Funchal city in the Madeira region (PT) (Figure 11). Such a network consists of a number of links (l) and nodes (n) as 2844 and 2801, respectively, and it is supplied by five reservoirs, whose water levels were assumed as constant. The total flowrate supplied by the reservoirs is 300 L/s. Out of five reservoirs, two are located at 360 and 400 m.a.s.l., whereas the elevation of the remaining three reservoirs is approximately 250 m.a.s.l. The network layout is presented in

Application to the Water Distribution Network (WDN) of Funchal City-Madeira, Portugal
Once the design and performance of the novel turbine were investigated by experiments and CFD analyses, an application to a real case study network was then considered. The case study concerns the water distribution network supplying Funchal city in the Madeira region (PT) (Figure 11). Such a network consists of a number of links (l) and nodes (n) as 2844 and 2801, respectively, and it is supplied by five reservoirs, whose water levels were assumed as constant. The total flowrate supplied by the reservoirs is 300 L/s. Out of five reservoirs, two are located at 360 and 400 m.a.s.l., whereas the elevation of the remaining three reservoirs is approximately 250 m.a.s.l. The network layout is presented in Figure 11. According to the current configuration of the network, several pressure reducing valves (PRVs) and throttle control valves (TCVs) are installed with the aim of controlling the pressure in the whole network. Moreover, further valves are present in order to close some links and determine a kind of district area of the network.
Due to the large topographical variability, a significant energy potential is available. Hence, the installation of one or more turbines could allow for a large production of energy, as well as a large water savings due to the reduction of the overall pressure. Nevertheless, in the links where the flow is significantly low, the energy production is low as well, and installing a valve instead of a turbine could be a more viable solution. Hence, in this study, the optimal location of turbines and valves was simultaneously investigated. An optimization procedure is implemented with reference to the daily average flow condition, with the aim of searching for the best scale design for turbines and, simultaneously, the best number and location of turbines and valves in order to maximize both water and energy savings, and minimize the investment cost.  Figure 11. According to the current configuration of the network, several pressure reducing valves (PRVs) and throttle control valves (TCVs) are installed with the aim of controlling the pressure in the whole network. Moreover, further valves are present in order to close some links and determine a kind of district area of the network. Due to the large topographical variability, a significant energy potential is available. Hence, the installation of one or more turbines could allow for a large production of energy, as well as a large water savings due to the reduction of the overall pressure. Nevertheless, in the links where the flow is significantly low, the energy production is low as well, and installing a valve instead of a turbine could be a more viable solution. Hence, in this study, the optimal location of turbines and valves was simultaneously investigated. An optimization procedure is implemented with reference to the daily average flow condition, with the aim of searching for the best scale design for turbines and, simultaneously, the best number and location of turbines and valves in order to maximize both water and energy savings, and minimize the investment cost.

The Optimization Procedure
In this study, the optimization problem was coupled with the hydraulic resolution of the network in one single mathematical model. Instead of using a hydraulic external solver (e.g., EPANET) to compute flow through links and pressure at nodes, the hydraulic behavior of the network is modeled within the optimization procedure by means of the mass continuity equations at the nodes and momentum balance equations along the pipes. The mathematical model will therefore consist of both linear and non-linear inequalities constraining the optimization problem, as well as non-linear equations modeling the hydraulic resolution of the network. With regard to the variables of the optimization, these are both binary (i.e., the installation of a device within a branch of the network) and continuous (i.e., pressure in the nodes, the discharges within the pipes, head loss within the devices). The resulting mathematical model is therefore a mixed integer non-linear programming (i.e., MINLP). Due to the non-convexity of some constraints, among the existing MINLP solvers, the Solving Constraint Integer Programs (SCIP) [45] was chosen to perform the optimization, since it is a robust global optimization solver, suitable for general MINLP problems-no matter if the optimization problem is convex or non-convex.

The Optimization Procedure
In this study, the optimization problem was coupled with the hydraulic resolution of the network in one single mathematical model. Instead of using a hydraulic external solver (e.g., EPANET) to compute flow through links and pressure at nodes, the hydraulic behavior of the network is modeled within the optimization procedure by means of the mass continuity equations at the nodes and momentum balance equations along the pipes. The mathematical model will therefore consist of both linear and non-linear inequalities constraining the optimization problem, as well as non-linear equations modeling the hydraulic resolution of the network. With regard to the variables of the optimization, these are both binary (i.e., the installation of a device within a branch of the network) and continuous (i.e., pressure in the nodes, the discharges within the pipes, head loss within the devices). The resulting mathematical model is therefore a mixed integer non-linear programming (i.e., MINLP). Due to the non-convexity of some constraints, among the existing MINLP solvers, the Solving Constraint Integer Programs (SCIP) [45] was chosen to perform the optimization, since it is a robust global optimization solver, suitable for general MINLP problems-no matter if the optimization problem is convex or non-convex. SCIP relies on spatial branch-and-bound algorithm with linear relaxation, various heuristics, and bound-tightening procedures. This solver can handle only differentiable functions; hence, it cannot perform the optimization in the presence of absolute values, if/sign function, etc. Morani et al. (2021) [46,47] found proper formulations of both variables and constraints in a general optimal location problem, avoiding the use of any non-differentiable functions within the mathematical model. Most of these formulations will be used in this study and presented later on.

The Variables
The location of a general device within a water network consisting of l links and n nodes is represented by a binary variable (i.e., I T k and I V k for turbine and valve, respectively) assuming a value of 1 if the device is installed and 0 otherwise. As previously highlighted, since the hydraulic behavior of the network is modeled within the optimization procedure, flow through links (Q k ) and head at nodes (H i ) are continuous variables of the optimization problem. With regard to H i , it is constrained between a lower and an upper bound, as follows: where p min and p max represent the minimum and maximum pressure values allowable at each node, z i the elevation of the i-th (i = 1 . . . n) node, and γ the specific weight of water. When the node is a reservoir (r), the head H i (t) is a boundary condition, and the variable is given by the discharge q r , which flows into or out of the reservoir itself. Regarding the discharge (Q k ) flowing through the k-th (k = 1 . . . l) link, it is expressed according to Morani et al. (2021) [46], as follows: where q + k and q − k are, respectively, the positive and negative parts of the discharge, and Q max is the upper bound. With regard to ζ k , it is a binary variable modeling the direction of the flow within a branch, and it is equal to 1 if the discharge flows according to the direction of the pipe k (i.e., the flow is given by the only positive component q + k ) and 0 if otherwise (i.e., the flow is given by the only negative component q − k ). The head loss within the device (H T k and H V k for turbine and valve, respectively) is a further continuous variable of the optimization problem. As presented in Morani et al. (2021) [46], these variables can also be split into their positive and negative components, as follows: where H k max is the maximum value of head loss achievable within the devices. Further variables of the problem are related to the design of the turbines. The machine presented in the previous sections is a small-scale prototype; hence, its diameter (D) should be properly scaled to be suitable for such a real application. Such a diameter represents a further variable of the optimization problem and has is bounded as: (16) where D min and D max represent the lower and upper bound, respectively, of the variable D in the case of turbine installation (i.e., I T k = 1). A further design variable of the optimization problem is represented by the impeller thickness (α). As the diameter D, α is constrained as: (18) where α min and α max represent the lower and upper bound, respectively, of the variable α in the presence of a turbine (i.e., I T k = 1).

Non-Linear Constraints
Most of the non-linearities in the optimization problem are given by the equations modeling the turbine design. In addition to the diameter (D), the installed turbine should also be characterized in terms of flow (Q T ), efficiency (η T ) and producible power (P T ). With regard to the flow through the turbine, it can be obtained from Equation (5), as follows: In order to set Q T as the total discharge flowing through the k-th pipe if the turbine is installed, the following constraint was defined: (20) According to Equation (20), when the turbine is not installed in the k-th pipe (i.e., I T k = 0), the Q T is forced to be equal to 0.
The turbine efficiency can be expressed according to Equation (7) as a function of (H T+ k + H T− k ): The producible power is hence given by the following expression: where η G is the generator efficiency, assumed equal to 90%. As constraints (19)- (22) are written with reference to each link of the network, the total number of such constraints is equal to 4l.
Further non-linearities are represented by the hydraulic equations modeling the resolution of the network, such as the mass continuity equations within the nodes, and momentum balance equations along the links.
Given a node i (i = 1 . . . n), the mass continuity equation can be expressed as follows: where q + k + q − k is the total discharge in the k-th link belonging to set K i of links approaching the node i, and the superscripts in and out indicate whether the discharge flows into or out of the i-th node, respectively. The term f i (H i − z i ) β represents the amount of leaked water at the i-th node, which will be better defined in Section 4.1.4. As the constraint (23) was written for each node i (i = 1 . . . n), the number of mass continuity equations can be accounted for as n.
Given a link k (k = 1 . . . l), the momentum balance equation is given by the following expression: where H i and H j are the head at the initial (i-th) and final (j-th) node of the k-th link, and r k L k is the head loss due to the resistance along the k-th pipe presenting a length L k . The unit head loss r k was calculated by the Hazen Williams formula: where C k and D k are the roughness coefficient and the diameter of the k-th pipe, respectively. As constraint (24) is expressed for each link k (k = 1 . . . l), the total number of such constraints can be accounted for as l.

Linear Constraints
Some linear constraints were used and were already presented in Section 4.1.1 (see Equations (9)- (14)) to define the variables of the problem. An additional set of l linear constraints (i.e., one for each link k) was defined in order to avoid the installation of both turbine and valve within a same link: Finally, further linear constraints were introduced in order to push the head loss within the devices to be equal to 0 if the device is not installed, and to be bounded otherwise.
As constraints (27) and (28) defined for each link k of the network, the total number of such devices is equal to 4l.

The Objective Function
As already highlighted, the aim of the optimization is to search for the best design and location of both turbines and valves, in order to maximize both the energy and water savings and minimize the investment cost. As in Morani et al. (2021) [46,48], the objective function that can be suitable for such an optimization problem is the net present value (NPV) of the investment, whose expression is given in Equation (29).
In Equation (29) the terms c T k I T k and c V k I V k represent, respectively, the total cost of turbines and valves occurring only at the beginning of the investment, whereas E p y and W s y are the annual inflow cash due, respectively, to the energy production and reduction of leaked water. Moreover, Y is the total number of y years, assumed as 10 years, and λ is the discount rate (set as 5%). In order to evaluate the annual income resulting from the energy production, an energy unit selling price equal to 0.1 €/kWh was assumed. With regard to the cost of the valve, the model cost suggested by Garcia et al. (2019) [49] was adopted, according to which the total cost of PRVs is dependent on pipe diameter ( Figure 12 With reference to the total cost of the turbines, since the analyzed turbine is a prototype, a function to properly model its cost is not actually available. Hence, in this study, the cost models of both a traditional turbine and a pump as turbine (PAT) (Figure 13) were alternatively used to carry out the optimization, since the effective cost model of the novel turbine is likely to comprise these two cost models. Among all the models presented in Figure 13, in this study the costs of the devices were evaluated according to the interpolations of experimental points presented by the authors in [50] (i.e., black dots for turbine unit costs and grey dots for PAT unit costs). With reference to the total cost of the turbines, since the analyzed turbine is a prototype, a function to properly model its cost is not actually available. Hence, in this study, the cost models of both a traditional turbine and a pump as turbine (PAT) (Figure 13) were alternatively used to carry out the optimization, since the effective cost model of the novel turbine is likely to comprise these two cost models. Among all the models presented in Figure 13, in this study the costs of the devices were evaluated according to the interpolations of experimental points presented by the authors in [50] (i.e., black dots for turbine unit costs and grey dots for PAT unit costs). type, a function to properly model its cost is not actually available. Hence, in this study, the cost models of both a traditional turbine and a pump as turbine (PAT) (Figure 13) were alternatively used to carry out the optimization, since the effective cost model of the novel turbine is likely to comprise these two cost models. Among all the models presented in Figure 13, in this study the costs of the devices were evaluated according to the interpolations of experimental points presented by the authors in [50] (i.e., black dots for turbine unit costs and grey dots for PAT unit costs). Figure 13. Unit cost per rated power of PATs and conventional turbines coupled to generators [50], compared to other literature cost models [51][52][53][54] .
With regard to the annual income due to water saving (W y s ) in Equation (29), it was evaluated as the difference between the total volume of water leaked without and in presence of devices within the network [46,48]: The water unit cost c w was fixed as 0.3 €/m 3 , as representative of the analyzed case study area. According to the formulation proposed by Araujo et al. (2006) [55], the total Figure 13. Unit cost per rated power of PATs and conventional turbines coupled to generators [50], compared to other literature cost models [51][52][53][54].
With regard to the annual income due to water saving (W s y ) in Equation (29), it was evaluated as the difference between the total volume of water leaked without and in presence of devices within the network [46,48]: The water unit cost c w was fixed as 0.3 €/m 3 , as representative of the analyzed case study area. According to the formulation proposed by Araujo et al. (2006) [55], the total leaked discharge in the network Q l can be evaluated as the sum of the discharge leaked in all single nodes: where β is an exponent depending on both the material and the shape of the orifice [56], which was set equal to 1.18, according to Araujo et al. (2006) [48,55]. With regard to f i , it represents a leakage coefficient whose expression is following presented: where c is a coefficient equal to 0.00001 L/s m 1+β [48,55] and K i is the set of nodes linked to the i-th node by a pipe with length as L i,j .

The Mathematical Model
The combination of all variables, objective function and linear/non-linear constraints presented in the previous sections results in the mathematical model in Equation (33).
Regarding the model in Equation (33), H k max , that is the upper bound of the variables head losses within the devices, was set as 200 m; p min and p max were set as 25 m and 150 m, respectively. Regarding the bounds of the impeller diameter of the turbines, the lower bound D min was set as 50 mm, whereas the upper bound D max varies for each pipe k and is equal to four times the diameter of the pipe itself. Finally, the upper bound of the variable discharges, which is Q max , was set as 500 L/s, whereas α min and α max were assumed equal at 10% and 200%, respectively.  [46], two small feasibility tolerances (tol H k , tol Q i ) were introduced within the momentum balance equation and mass continuity equation, respectively, in order to push the solver to achieve a solution. Such tolerances are small quantities (tol H k = 0.01 m and tol Q i = 0.005 L/s) and can be considered as acceptable as all other uncertainties affecting the accuracy of the problem (e.g., evaluation of end-used demands, use of Hazen Williams formula).
The total number of variables amounts to 42,617, whereas the sum of both linear and non-linear constraints is accounted for as 51,149. Table 2 shows the results of the optimization obtained by using the cost models presented in Section 4.1.4. In Table 2, (I) and (II) refer to results obtained when a traditional turbine and a PAT cost model, respectively, are employed to evaluate the cost of the installed turbines. When the turbine cost is evaluated according to the traditional turbine cost model (see (I)), the optimization solver installs a number of turbines and valves as 5 and 4, respectively, with an investment cost of EUR 482,298.

Optimization Results
The resulting average produced power was accounted for as 68 kW, and the amount of water saved per day is equal to 8507 m 3 . When the literature PAT cost model is, instead, employed to evaluate the cost of the installed turbines (see (II)), the solver selects more turbines than valves (i.e., 6 against 1- Figure 14), since the installation of valves could not be convenient compared to turbines, especially where the pipe diameters are large. However, compared to the results in case (I), the increased value of NPV in case (II) (i.e., EUR 11,942,920) is due to the larger water savings, which is equal to 13,673 m 3 per day, and smaller investment costs (EUR 55,786). Regarding the power production, it is quite similar between the two cases, obtaining around 65 kW. When the turbine cost is evaluated according to the traditional turbine cost model (see (I)), the optimization solver installs a number of turbines and valves as 5 and 4, respectively, with an investment cost of EUR 482,298.
The resulting average produced power was accounted for as 68 kW, and the amount of water saved per day is equal to 8507 m 3 . When the literature PAT cost model is, instead, employed to evaluate the cost of the installed turbines (see (II)), the solver selects more turbines than valves (i.e., 6 against 1- Figure 14), since the installation of valves could not be convenient compared to turbines, especially where the pipe diameters are large. However, compared to the results in case (I), the increased value of NPV in case (II) (i.e., EUR 11,942,920) is due to the larger water savings, which is equal to 13,673 m 3 per day, and smaller investment costs (EUR 55,786). Regarding the power production, it is quite similar between the two cases, obtaining around 65 kW. The turbines are installed in both the cases in the middle region of the network, whereas the valves are installed in the peripheral area, where the water saving potential is large but the recoverable energy is contained. Tables 3 and 4 show the results in terms of turbine geometry and operating conditions at BEP for case (I) and case (II), respectively. With reference to N, it is the rotational speed resulting from the optimal diameter (D) and head drop (H T ) according to Equation (6).  The turbines are installed in both the cases in the middle region of the network, whereas the valves are installed in the peripheral area, where the water saving potential is large but the recoverable energy is contained. Tables 3 and 4 show the results in terms of turbine geometry and operating conditions at BEP for case (I) and case (II), respectively. With reference to N, it is the rotational speed resulting from the optimal diameter (D) and head drop (H T ) according to Equation (6).  Figure 15a shows the pressure contour plots according to case (I) when devices are installed within the network: the installation results in a significant reduction in the overall excess pressure when compared with the scenario without installed devices (Figure 15b). The same is valid for case (II), whose pressure contour plot is presented in Appendix A ( Figure A1).   Figure 15a shows the pressure contour plots according to case (I) when devices are installed within the network: the installation results in a significant reduction in the overall excess pressure when compared with the scenario without installed devices ( Figure  15b). The same is valid for case (II), whose pressure contour plot is presented in Appendix A ( Figure A1).
The produced power, as well as the head drop within the installed devices, for case (I) is shown in Figure 16a,b, respectively. According to Figure 16b, the reduction in head resulting from the device installation is significant, even reaching a value of 160 m. The same figure for case (II) is presented in Appendix A ( Figure A2).  The produced power, as well as the head drop within the installed devices, for case (I) is shown in Figure 16a,b, respectively. According to Figure 16b, the reduction in head resulting from the device installation is significant, even reaching a value of 160 m. The same figure for case (II) is presented in Appendix A ( Figure A2).

Conclusions
In this paper a new centrifugal micro-turbine (CMT) was investigated in terms of performance and design. Based on a CFD modeling coupled with experimental tests on a turbine prototype, the unstable phenomena affecting the turbine components during its operation have been investigated in detail. In particular, the velocity and pressure profiles within turbine components were presented for different operating conditions. According to the results, the difference in terms of pressure drop at inlet and outlet of the rotor is significantly small, except for rotational speeds greater than 1850 rpm. According to the velocity profiles at the runners, high flow acceleration was detected from the leading edge to the trailing edge, followed by a significant flow deceleration due to an increase in static pressure. Moreover, the degree of flow deceleration was lower for high rotational speeds, at which the internal losses were therefore lower. A maximum velocity value of 12 m/s was detected in the region close to the leading edge of the rotating blade, whereas the minimum velocity of around 1 m/s is presented at the right side of the vane. However, the validation of the CFD numerical simulation was proven by a good agreement with experimental data in terms of global parameters. Based on affinity laws and the experimental tests, some relationships were obtained to design the centrifugal turbine at its best efficiency point. Finally, based on this new centrifugal micro-turbine design, a real water distribution network (WDN), in the Funchal system, was assumed as the case study and an optimization procedure was carried out in order to search for an optimal scale design of turbines, as well as an optimal location of such designed turbines together with pressure reducing valves (PRV), in order to maximize the energy production and minimize the leakage water volume. Due to the lack of information about the purchase cost of such a new centrifugal turbine, two different literature cost models were assumed to evaluate the investment cost: the effective solution will definitely comprise the two obtained solutions. According to the results, the use of a PAT cost model pushes the solver to install more turbines than valves. Moreover, compared to the solution resulting from the use of a traditional turbine cost model, when the turbine cost is evaluated according to a literature PAT cost model, the resulting water savings is large (i.e., 13,673 m 3 per day against 8507 m 3 per day) and the investment cost is significantly more confined (EUR 55,786

Conclusions
In this paper a new centrifugal micro-turbine (CMT) was investigated in terms of performance and design. Based on a CFD modeling coupled with experimental tests on a turbine prototype, the unstable phenomena affecting the turbine components during its operation have been investigated in detail. In particular, the velocity and pressure profiles within turbine components were presented for different operating conditions. According to the results, the difference in terms of pressure drop at inlet and outlet of the rotor is significantly small, except for rotational speeds greater than 1850 rpm. According to the velocity profiles at the runners, high flow acceleration was detected from the leading edge to the trailing edge, followed by a significant flow deceleration due to an increase in static pressure. Moreover, the degree of flow deceleration was lower for high rotational speeds, at which the internal losses were therefore lower. A maximum velocity value of 12 m/s was detected in the region close to the leading edge of the rotating blade, whereas the minimum velocity of around 1 m/s is presented at the right side of the vane. However, the validation of the CFD numerical simulation was proven by a good agreement with experimental data in terms of global parameters. Based on affinity laws and the experimental tests, some relationships were obtained to design the centrifugal turbine at its best efficiency point. Finally, based on this new centrifugal micro-turbine design, a real water distribution network (WDN), in the Funchal system, was assumed as the case study and an optimization procedure was carried out in order to search for an optimal scale design of turbines, as well as an optimal location of such designed turbines together with pressure reducing valves (PRV), in order to maximize the energy production and minimize the leakage water volume. Due to the lack of information about the purchase cost of such a new centrifugal turbine, two different literature cost models were assumed to evaluate the investment cost: the effective solution will definitely comprise the two obtained solutions. According to the results, the use of a PAT cost model pushes the solver to install more turbines than valves. Moreover, compared to the solution resulting from the use of a traditional turbine cost model, when the turbine cost is evaluated according to a literature PAT cost model, the resulting water savings is large (i.e., 13,673 m 3 per day against 8507 m 3 per day) and the investment cost is significantly more confined (EUR 55