Water: Slip Factor Correction in 1-D Performance Prediction Model For Pats

Download as pdf or txt
Download as pdf or txt
You are on page 1of 23

water

Article
Slip Factor Correction in 1-D Performance Prediction
Model for PaTs
Tommaso Capurso * , Michele Stefanizzi , Giuseppe Pascazio , Sergio Ranaldo,
Sergio M. Camporeale , Bernardo Fortunato and Marco Torresi
Department of Mechanics, Mathematics and Management (DMMM), Polytechnic University of Bari,
70125 Bari, Italy; [email protected] (M.S.); [email protected] (G.P.);
[email protected] (S.R.); [email protected] (S.M.C.); [email protected] (B.F.);
[email protected] (M.T.)
* Correspondence: [email protected]; Tel.: +39-327-915-0352

Received: 29 December 2018; Accepted: 13 March 2019; Published: 19 March 2019 

Abstract: In recent years, pumps operated as turbines (PaTs) have been gaining the interest of
industry and academia. For instance, PaTs can be effectively used in micro hydropower plants (MHP)
and water distribution systems (WDS). Therefore, further efforts are necessary to investigate their
fluid dynamic behavior. Compared to conventional turbines, a lower number of blades is employed
in PaTs, lowering their capability to correctly guide the flow, hence reducing the Euler’s work; thus,
the slip phenomenon cannot be neglected at the outlet section of the runner. In the first part of the
paper, the slip phenomenon is numerically investigated on a simplified geometry, evidencing the
dependency of the lack in guiding the flow on the number of blades. Then, a commercial double
suction centrifugal pump, characterized by the same specific speed, is considered, evaluating the
dependency of the slip on the flow rate. In the last part, a slip factor correlation is introduced based
on those CFD simulations. It is shown how the inclusion of this parameter in a 1-D performance
prediction model allows us to reduce the performance prediction errors with respect to experiments
on a pump with a similar specific speed by 5.5% at design point, compared to no slip model, and by
8% at part-loads, rather than using Busemann and Stodola formulas.

Keywords: PaT; slip phenomenon; performance prediction; OpenFOAM; CFD

1. Introduction
The world electricity demand is increasing year by year, and the shift to renewable energy from
fossil fuels is attracting many investments; at the same time, many financial resources are being spent
to improve the efficiency of existing power plants and energy distribution systems. Hydropower,
which comprises hydroelectric power plants and water distribution systems, is the main field where,
in recent decades, pumps as turbines (PaTs) have been installed en masse. Many researchers strongly
suggest the use of PaTs, i.e., conventional pumps used in reverse mode to recover energy from fluids
subjected to a considerable pressure drop, as an alternative to throttling devices. PaTs are becoming
more and more tempting because pumps are mass-produced, cheaper than conventional turbines and
they cover a wide range of specific speeds and sizes.
Some technical reports by international pump manufacturers, like Sulzer and KSB, started to be
published in the early 80s with Laux [1], who studied reverse-running multistage pumps as energy
recovery turbines in oil supply systems; and Apfelbacher et al. [2] also studied the application of a PaT
in a reduction pressure station in Aachen (Germany).
Large hydropower plants feed the national grid, whereas typical off-grid micro hydropower
plants (MHP), in the range of 5–100 kW, are the most popular solutions for electrification among

Water 2019, 11, 565; doi:10.3390/w11030565 www.mdpi.com/journal/water


Water 2019, 11, 565 2 of 23

rural communities [3]. In MHP, hydraulic turbines represent the critical technological component
for energy production. Their cost significantly impacts the project budget. PaTs fulfill an important
role in hydraulic pumped-storage plants to supply power peak demands; indeed, the oscillations of
the power produced by renewable sources are going to promote the need for energy storage such as
micro hydro, where the same machine needs to be operated both in pump mode and in turbine mode.
Furthermore, PaTs can represent stand-alone small hydro systems in all those remote communities,
where the grid connection is not economical and/or practically possible. In this case, PaTs can satisfy
the energy demand by exploiting natural water sources.
Many works assess the installation of PaTs into Water Distribution Networks (WDNs) in order to
save energy dissipated by pressure reducing valves (PRVs). PaTs found applications not only in the
hydraulic sector, but also in other fields of the process engineering.
Nowadays, the prediction of the performance of a PaT is crucial in order to select the most
suitable machine, starting from its geometric parameters and its specific speed, since pumps are
not usually tested in turbine mode by manufacturers. Because of the lack of experimental data,
considerable attempts have been made by several authors in the prediction of PaT performance,
not only from a theoretical point of view, but also with experimental campaigns of pumps tested in
reverse mode, as shown by Barbarelli et al. [4], Derakhshan and Nourbakhsh [5], Nautiyal et al. [6],
Pugliese et al. [7], Rossi et al. [8], Singh and Nestmann [9], Shi et al. [10], Yang et al. [11], Tan and
Engeda [12], Carravetta et al. [13].
Several models have been proposed for the prediction of the characteristic curve of a PaT. Some
models, based on statistical and normalization approaches or artificial neural networks [14,15], are
adopted when the geometric parameters of the runner are unknown [16], and they are useful to
evaluate the flow rate and the head exploited by the PaTs running at their best efficiency point.
At the same time, other models are based on theoretical approach as for the evaluation of the flow
characteristics through the machine, as proposed by Barbarelli et al. [17], Gülich [18]. These are usually
used by manufacturers, which have detailed information about their machines. The latter models
currently show a lack of accuracy in the prediction of the characteristic curve, especially at part loads
and over loads because of the simplification of the impeller channel and the volute geometry. For this
reason, manufacturers need a tool that could support them not only to predict the turbine mode
operation, but also to design ad hoc PaTs for specific applications by investigating more deeply the
hydraulic behavior of their machines running in reverse mode.
Similarly to in centrifugal pumps, a slip phenomenon can also occur in turbines with a finite
number of vanes, at the outlet section of the runner. As shown by Ventrone [19], the relative velocity
vectors at the outlet of a centripetal turbine (e.g., Francis turbines and PaTs) are subject to a deflection
with respect to the direction defined by the blade congruent angle (slip). Ventrone [19] and Shi et al. [20]
justify the slip phenomenon by demonstrating the presence of fluid flow vorticity between the blades
in planes orthogonal to the rotational axis. Thus, counter rotating vortices appear inside the impeller
channels. This phenomenon shall be more evident in radial PaTs, where the number of blades is lower
than in conventional turbines. In mixed axial-radial runners, additional dynamic effect should be
taken into account.
The resulting flow deflection leads to an increase in the absolute tangential velocity and
consequently a reduction of the extracted hydraulic energy with respect to the available one. Thus,
in this work, the slip phenomenon in PaTs was studied in detail with the aim at focusing on the
influence of the slip phenomenon on the performance prediction of PaTs. Regarding this aspect, the
slip phenomenon is often neglected in prediction models and there are few works in the literature,
which have investigated in detail the effect of the slip phenomenon in PaTs. Indeed, the evaluation of
the slip effect is difficult and there are not many experimental correlations.
Differently from a previous work [21], herein the authors extended their investigations on the
slip phenomenon at the outlet of PaTs. Arranging their work in three parts, the first, (Sections 2
and 3) and the last (Sections 5 and 6) parts being new. Moreover, additional details with respect to
Water 2019, 11, 565 3 of 23

the previous work are added. Initially, the reader is guided through a description of the problem
(from a phenomenological point of view). Indeed, numerical investigations have been carried out
on a simplified PaT runner geometry, designed according to potential flow solution, in order to
study the slip phenomenon and be able to more easily describe the sources of the flow deflection
(see Sections 2 and 3). In the second part, the work deals with the prediction of the performance of
a real commercial PaT, numerically simulated in both direct and inverse mode, with an emphasis
on the flow field through the impeller [21] (see Section 4). The proposed numerical investigation
aims to put in evidence the slip phenomenon on a real configuration. The computational domain
reproduces the entire geometry of the commercial pump, including the suction pipe, the impeller
and the discharge pipe The introduction of a slip factor, σturbine , whose definition is equivalent to the
one used for centrifugal pumps [22], accounts for the difference between the actual work and the
theoretical one. This analysis is based on the evaluation of the axial and tangential velocities at the
mean radius of the trailing edge of the blades at the outlet of the centripetal impeller. Moreover, this
analysis shows the presence of different phenomena, e.g., pressure gradient and secondary flows,
which may occur at the outlet of a centripetal turbine and affect the performance of the machine.
In the last part, a slip factor (σturbine ) correlation has been derived from these numerical results.
The definition of σturbine provides a great support in developing and tuning 1-D models for the
prediction of the turbine performance and it could provide the guidelines to successfully improve
PaT geometries. A relationship between the new parameter (σturbine ) and the flow rate is presented
and discussed. This new correlation has been applied to the 1-D prediction model proposed by
Stefanizzi et al. [22] in order to predict the entire characteristic of the same PaT, starting from the
knowledge of its geometrical parameters. Finally, the same model has been applied to another real
PaT with a similar specific speed (nq ) in order to compare the curve predicted by the model with the
experimental one. In conclusion, the results of the new model show a more accurate prediction of the
PaT performance under design and part-load conditions compared to the results of the same prediction
models, employing slip factor correlations already available in the literature [23,24] (see Sections 5
and 6). The work focuses on a specific speed hydraulic turbomachinery; thus, it does not aim to be
exhaustive—to cover the entire range of specific speeds, further analyses are required.

2. The Slip Phenomenon in Turbomachinery


As stated in the introduction, in turbines with a limited number of vanes, as in the case of PaT,
a slip phenomenon occurs at the outlet section of the runner. The main effect of this phenomenon is
related to the reduction of the Euler’s work, Y, due to the flow deflection increasing the kinematic
relative flow angle, β 1 , with respect to the outlet congruent blade angle, β 1B (Figure 1). Unfortunately,
there are few works in the literature that have investigated the effect of the slip phenomenon in power
generating turbomachinery in detail—particularly in PaTs. Indeed, the evaluation of the slip effect is
difficult, and there are not many experimental correlations.
Among the sources of the flow deflection, the vortex in the blade-to-blade plane (in particular in
purely radial runners, as shown by Ventrone [19] and Shi et al. [20]) and the pressure gradient due to
the blade loading together with the Coriolis force should be mentioned.
The flow inside the mixed axial-radial impeller is subjected to those complex phenomena; however,
it is evident that the lower the number of blades, higher the blade loads, and hence the greater the
flow deflection.
Figure 1 represents this phenomenon by showing both the actual (solid line) and the blade
congruent velocity triangle (dashed line) at the runner outlet section. It is important to highlight that
in this work, the subscripts 1 and 2 refer to the runner outlet and inlet, respectively. Moreover, all the
quantities (pressure and velocities) will be considered as values averaged on the inlet and outlet areas.
The flow deflection at the turbine outlet section 1 involves the increase of the tangential component
Water 2019, 11, 565 4 of 23

of the absolute velocity, cu1 , causing a decrease in terms of specific work, Y, obtained by the turbine,
according to Euler’s equation
Water 2019, 11, x FOR PEER REVIEW 4 of 24
Y = u2 cu2 − u1 cu1 (1)
Water 2019, 11, x FOR PEER REVIEW 4 of 24

Figure 1. Velocity triangles at the inlet (2) and at the outlet of the runner (1). Both actual (solid line)
and blade
Figure
Figure 1. congruent
Velocity
1. Velocity (red dashed
triangles
triangles at the line)
at the
inlet velocity
inlet
(2)(2)
and attriangles
and at the
the areofrepresented
outlet
outlet of
thethe at (1).
runner
runner the
(1). outlet.
Both
Both actual
actual (solid
(solid line)
line)
and blade
and congruent
blade (red
congruent dashed
(red line)
dashed velocity
line) triangles
velocity areare
triangles represented at the
represented outlet.
at the outlet.
Neglecting secondary flows effects, the slip phenomenon at the outlet of mixed axial-radial
runners is mainly
Neglecting
Neglecting due to the
secondary
secondary flowsblade
flows loading.
effects,
effects, Considering
thethe
slip phenomenon
slip the flow
phenomenon at atatthe
the the PaTofoutlet
outlet
outlet mixed
of (Figure
mixed 2) with
axial-radial
axial-radial
tangential
runners
runners is is and
mainly
mainly meridional
due to to
due thethe velocity
blade
blade components
loading.
loading. and no
Considering
Considering the radial
theflow
flowatcomponent,
thethe
at PaTPaT the(Figure
outlet
outlet fluid
(Figureelement,
2) 2)
with
with
displayed
tangential in
and Figure 2,
meridional is subjected
velocity to a
components radialandCoriolis
no acceleration,
radial component,
tangential and meridional velocity components and no radial component, the fluid element, which
the does
fluid not
element,contribute
displayed to
modifying
in displayed
Figure 2, is the pressure
insubjected
Figure 2, togradient along
is asubjected
radial Coriolisthe tangential
acceleration,
to a radial direction,
Coriolis which but rather
does not
acceleration, along
which the
contribute radial
does not direction.
to contribute
modifying to
Figure
the 2 schematically represents the acceleration terms in a qualitative
modifying the pressure gradient along the tangential direction, but rather along the radialFigure
pressure gradient along the tangential direction, but rather along the radial way. Based
direction. on this2
direction.
representation,
schematically the Coriolis
represents the acceleration
acceleration may
terms in introduce
a qualitativesecondary
way.
Figure 2 schematically represents the acceleration terms in a qualitative way. Based on this flows,
Based on which
this will not
representation, be
investigated
therepresentation, in this
Coriolis accelerationthework.
mayUnder
Coriolis thesesecondary
introduce
accelerationsimplifications,
may flows, the Coriolis
which
introduce will not
secondaryacceleration
be flows, does in
investigated
which not influence
this
will work.
not be
the
Under pressure
these gradient along
simplifications, the
the tangential
Coriolis direction;
acceleration conversely,
does not it
influenceis the
the
investigated in this work. Under these simplifications, the Coriolis acceleration does not influence main source
pressure of fluid
gradient flow
along
deflection
thethe
tangentialin centrifugal
pressuredirection; pumps.
gradient conversely, it is the main
along the tangential source conversely,
direction; of fluid flowitdeflection
is the main in centrifugal pumps.
source of fluid flow
deflection in centrifugal pumps.

Figure2.2.Representation
Figure Representationofofthe
themain
mainvelocity
velocityvectors
vectorsand
andaccelerations
accelerationsreferred
referredtotothe
theparticle
particleofoffluid
fluid
atatFigure
the
theoutlet
outletof
ofthe
therunner.
runner.The
Thevector
vectormagnitudes
magnitudes are
aredrawn
drawn qualitatively.
qualitatively.
2. Representation of the main velocity vectors and accelerations referred to the particle of fluid
at the outlet of the runner. The vector magnitudes are drawn qualitatively.
Although
Althoughthethe above-mentioned
above-mentioned works [19,20]
works have introduced
[19,20] the slip phenomenon
have introduced in hydraulicin
the slip phenomenon
turbines,
hydraulic it has often it
turbines, been
hasneglected
often beenin the development
neglected of
in thehave 1-D theoretical
development models, specifically in
Although the above-mentioned works [19,20] introducedof the1-D theoretical
slip models,
phenomenon in
predicting
specificallyPaT
in performance
predicting in
PaT turbine mode
performance operation
in turbine by taking
mode into
operationaccount
by the
taking geometry
into of
account the
the
hydraulic turbines, it has often been neglected in the development of 1-D theoretical models,
geometry of in
specifically thepredicting
machine andPaTsome complexin
performance phenomena
turbine modelike operation
hydraulic bylosses.
takingObviously, a deep
into account the
insight
geometryintoofthis
the phenomenon in realcomplex
machine and some impellers could improve
phenomena the accuracy
like hydraulic losses.ofObviously,
1-D theoretical
a deep
approaches
insight intoin this
the performance
phenomenonprediction of pumpscould
in real impellers operating
improveas turbines.
the accuracy of 1-D theoretical
approaches in the performance prediction of pumps operating as turbines.
Water 2019, 11, 565 5 of 23

machine and some complex phenomena like hydraulic losses. Obviously, a deep insight into this
Water 2019, 11, xin
phenomenon FOR PEER
real REVIEW could improve the accuracy of 1-D theoretical approaches in5 the
impellers of 24

performance prediction of pumps operating as turbines.


3. Influence of the Number of Blades on the Slip Phenomenon
3. Influence of the Number of Blades on the Slip Phenomenon
In this section, an early-stage evaluation of the slip phenomenon is carried out by simulating a
simplified archetype
In this section, anof a centrifugal
early-stage pump vane
evaluation byslip
of the onlyphenomenon
changing theisnumber
carried of outvanes (namely 6,a7
by simulating
and 8) with
simplified the purpose
archetype of quantifying
of a centrifugal pump thevane
effect
byof thechanging
only pressure thegradient
number acting in the(namely
of vanes tangential
6,
7direction.
and 8) with Thisthegeometry
purposedoes not consider
of quantifying thethe rounded
effect of thetrailing
pressure edge and itsacting
gradient effects.
in the tangential
The This
direction. simulation
geometry campaign
does notwas limited
consider the to single trailing
rounded vanes (under
edge and theitssimplified
effects. hypothesis of
periodic behavior).campaign
The simulation These geometries
was limitedwithto different
single vanes numbers
(underof blades
the werehypothesis
simplified designed by means of
of periodic
the modelThese
behavior). proposed by d’Agostino
geometries [25] and
with different investigated
numbers of blades bywere
means of numerical
designed by means simulations.
of the modelThe
d’Agostino model expresses the 3D incompressible, inviscid, irrotational
proposed by d’Agostino [25] and investigated by means of numerical simulations. The d’Agostino flow through helical blades
with slow
model axial the
expresses variations of the pitchinviscid,
3D incompressible, and backsweep by superposing
irrotational flow throughahelical2D cross-sectional
blades with slowaxial
vorticity
axial correction
variations of theto apitch
fullyand
guided flow with
backsweep byaxisymmetric
superposing stagnation velocity in axial
a 2D cross-sectional the meridional
vorticity
plane. It isto
correction important to highlight
a fully guided flow that
withthe three geometries
axisymmetric have novelocity
stagnation blade thickness, since in this
in the meridional early
plane.
Itstage of the sliptophenomenon
is important highlight that analysis, the geometries
the three rounded leading have edges of the
no blade runner could
thickness, sincecomplicate the
in this early
investigation.
stage of the slip phenomenon analysis, the rounded leading edges of the runner could complicate
All the runners have the same specific speed, nq = 20, and the value for the absolute tangential
the investigation.
velocity calculated
All the runnersby havemeans of thespecific
the same blade congruent
speed, nq theory at the
= 20, and thecenter
valueline (𝑐 absolute
for the ) at the outlet of the
tangential
impeller vane is 𝑐 = 4.78 m/s.
velocity calculated by means of the blade congruent theory at the center line (cline ) at the outlet of the
impellerThevane
specific
is cu1speed
= 4.78 is defined
m/s. as
The specific speed is defined as p𝑛 𝑄 /𝑓
𝑛 n= Q BEP / fq (2)
nq = 𝐻 .
0.75
(2)
H
where 𝐻 is the hydraulic head (m), 𝑄 is the flow rate (m 3/s), 𝑛 the rotational speed (rpm) and 𝑓
where H is the hydraulic head (m), Q is the flow rate (m3 /s), n the rotational speed (rpm) and f q
is the number of impeller entries (𝑓 = 2 for a double suction impeller). The previous expression
is the number of impeller entries ( f q = 2 for a double suction impeller). The previous expression
(Equation (2)) has been adopted to calculate both centrifugal pump and PaT specific speeds, since it
(Equation (2)) has been adopted to calculate both centrifugal pump and PaT specific speeds, since it
is commonly used in empirical prediction model to forecast the BEP under turbine mode operation
is commonly used in empirical prediction model to forecast the BEP under turbine mode operation
starting from the pump specific speed and vice versa.
starting from the pump specific speed and vice versa.
Further geometric details are given in Table 1. The three single-vaned computational domains,
Further geometric details are given in Table 1. The three single-vaned computational domains,
designed according to d’Agostino, with a straight pipe at the inlet, are displayed in Figure 3.
designed according to d’Agostino, with a straight pipe at the inlet, are displayed in Figure 3.
Table 1. Geometric data about the centrifugal pump archetypes studied in this paragraph.
Table 1. Geometric data about the centrifugal pump archetypes studied in this paragraph.
D2 (mm) l2 (mm) Zblades β2B (°) β1B (°) ◦n (rpm)
D2 (mm) l2 (mm) Zblades β2B (◦ ) β1B ( ) n (rpm)
252 12 6-7-8 17 25 3900
252 12 6-7-8 17 25 3900

Figure 3. Three single vane geometries of the selected centrifugal pump with different number
of blades.
Figure 3. Three single vane geometries of the selected centrifugal pump with different number of
The meshes were generated by means of Pointwise Gridgen® using unstructured tetrahedral
blades.
elements. To guarantee the same grid density, a fixed node spacing has been applied to the connectors.
The meshes were generated by means of Pointwise Gridgen® using unstructured tetrahedral
elements. To guarantee the same grid density, a fixed node spacing has been applied to the
connectors. This means that the wider the vane (i.e., lower the number of blades), the larger the
number of cells. The total number of cells for each geometry is reported in Table 2. The three
Water 2019, 11, 565 6 of 23

Water 2019,
Water 11, 11,
2019, x FOR PEER
x FOR REVIEW
PEER REVIEW 6 of6 24
of 24
This means that the wider the vane (i.e., lower the number of blades), the larger the number of cells.
computational
The
computational domains
total numberdomains areare
of cells displayed
for in in
Figures
each geometry
displayed 4 and
is reported
Figures 5,in
4 and where
5,Tabletheir
where details
2.their
The three at thethe
inlet
computational
details at and
inlet outlet
domains
and outlet
sections
are areare
shown.
displayed
sections in Figures 4 and 5, where their details at the inlet and outlet sections are shown.
shown.

Table 2. Number
Table
Table 2. of cells
Number
2. Number of of the
of cells
cells of centrifugal
of the
the pump
centrifugal
centrifugal vane
pump
pump archetypes
vane
vane studied
archetypes
archetypes in this
studied
studied in paragraph.
in this
this paragraph.
paragraph.
6 6Blades
Blades
6 Blades7 Blades 8 Blades
7 Blades
7 Blades 8 Blades
8 Blades
Discharge pipe
DischargeDischarge 6
pipe (×10 pipe
) (×10 6) 0.46
(×106) 0.46 0.46 0.40
0.40
0.40 0.34
0.34 0.34
Channel ×106 ) (×10
Channel
(Channel 6)
(×10 )6 0.80
0.80 0.80 0.70 0.70
0.70 0.60
0.60 0.60

Figure
Figure
Figure 4. Grids
4. Grids
4. Grids of the
of the
of the single
single vane
vane
single geometries
geometries
vane with
with
geometries different
different
with number
number
different of
of blades.
of blades.
number blades.

Figure 5. Grid details of the single vane geometry with 7 blades where the outlet and inlet section of
the pump under turbine operating conditions are specified.

Figure
For 5. Grid
this details of case,
the single vane geometry with 7 blades
the where thethe
outlet andand
inlet section ofFluent ®,
Figure 5. simplified
Grid details of thesimulations
single vane were run
geometry with
with commercial
7 blades where CFD
outlet code ANSYS
inlet section of
thethe
pump
applying pumpunder
a single turbine
reference
under turbineoperating
frame conditions
technique.
operating are
conditions specified.
Inare
the next paragraph, when bigger meshes have to be
specified.
analyzed, an open-source CFD code with no restriction on parallelization has been employed in order
to For this simplified case, simulations were run with thethe
commercial CFD code ANSYS Fluent ®,
reduce
For the
this computational
simplified case, time. With respect
simulations were torun
the centrifugal
with pump the
commercial inlet
CFD and
code the outlet
ANSYS as well
Fluent ®,

applying
as a single
the angular
applying a single reference
velocity frame
haveframe
reference technique.
been technique.
reversed. AtIn the next paragraph,
inlet paragraph,
In the next when
the absolute bigger
radial
when meshes
and tangential
bigger have
meshes have to be
velocity
to be
analyzed,
components
analyzed, an anopen-source
have beenCFD
open-source CFDcode
imposed codewith
with withnono
flowrestriction
turbulent
restriction onintensity
parallelization
on equal to
parallelization has
3%,
hasbeen employed
while
been uniform
employed inpressure
order
in order
to was
reduce
to set the
reduce upthe computational
at the outlet. Due
computational time.
to With
the
time. With respect
y+ value to to
thethe
calculated
respect centrifugal
on pump
the stator
centrifugal andthe
pump rotorinlet
the and
surfaces
inlet andthe +outlet
(ythe =150),
outletasthe
as
well
k-ωas
wellSSTthe angular
model
as the angular velocity
by Menter have
velocity[26]
havebeen
withbeen reversed.
thereversed. At
application the inlet
of the
At the the
wall
inlet absolute
thefunctions radial and
has been
absolute radial tangential
andused velocity
for turbulence
tangential velocity
components
closure.
components have
This have been
turbulence imposed
been modelwith
imposed theflow
is with turbulent
standard
flow intensity
for performing
turbulent intensityequal to to
numerical
equal 3%,3%,while
analyses
while uniform
in pressure
hydraulic
uniform turbo
pressure
was set up
was set up at
machinery. at the outlet. Due
the outlet. Due
It automatically to the
uses y + value calculated on the stator and rotor surfaces (y+ =150), the
to the yk-ωvalue
+ model calculated on the stator
in the near-wall and
regions, the k-ε model
rotor surfaces
whereas (y =150),
+ the
is used
k-ωaway
SST
k-ω model
SST model
from bywalls.
the Menter
by Menter
The[26]
k-ωwith
[26] with
SST the application
the
model can giveof
application anthe wall
ofaccurate
the functions
wall functions
predictionhas been
has
of been
flow used
usedforfor
separation,turbulence
turbulence
explaining
closure.
closure.
its This
common Thisturbulence
turbulence
use for themodel
model is the
numerical is thestandard
standard
investigationsforforperforming
performing
of numerical
flow inside numerical analyses
analyses
the centrifugal in hydraulic
in
pumps [27].turbo
hydraulic Aturbo
zero
machinery.
machinery.
shear It
stress It automatically
has automatically uses
been applieduses the
to the thek-ω model
k-ω belonging
walls in
model in thethe near-wall
to thenear-wallregions, whereas the k-ε
regions, whereas the k-ε model
stator parts. model is is
used away from the walls. The k-ω SST model can give an accurate
used away from the walls. The k-ω SST model can give an accurate prediction of flow separation, prediction of flow separation,
explaining
explaining its its
common
common useuse
forfor
thethe
numerical
numerical investigations
investigations of flow inside
of flow insidethethecentrifugal
centrifugal pumps
pumps [27].
[27].
A zero shear stress has been applied to the walls belonging
A zero shear stress has been applied to the walls belonging to the stator parts. to the stator parts.
Water 2019, 11, 565 7 of 23
Water 2019, 11, x FOR PEER REVIEW 7 of 24
Water 2019, 11, x FOR PEER REVIEW 7 of 24

Evaluation of theofPressure
Evaluation andand
the Pressure thethe
Velocity
VelocityDeflection
Deflection
Evaluation of the Pressure and the Velocity Deflection
At the Atend
the end of the
of the simulation,
simulation, theabsolute
the absolute tangential velocity
velocity(𝑐 ) and thethe
static pressure werewere
At the end of the simulation, the absolute tangential
tangential velocity (𝑐(c)u1and
) and static
the static pressure
pressure were
calculated
calculated at
at the the outlet of the channel for the three geometries. The results were calculated along an an
calculated at outlet
the outletof the channel
of the channelforforthe
thethree
three geometries.
geometries. The Theresults
resultswerewere calculated
calculated along along
an
arc atarc
arc
at mean
mean radius
at mean
radius
radius
(0.0435
(0.0435 m), m),
(0.0435 see see
m),
Figure
Figure
see 6a 6a
Figure and
6a
and
andonon iso-clips,see
iso-clips,
on
see Figure6b,
iso-clips, seeFigure
6b, in correspondence
Figure 6b, in
correspondence with
in correspondence with with the
the channel outlet section. As explained in Section 2, a flow deflection occurs at the outlet of the
channel
the outlet
channelsection. As explained
outlet section. in Section
As explained 2, a flow
in Section 2, a deflection occurs
flow deflection at the
occurs outlet
at the of the
outlet runner
of the
runner with respect to the blade congruent angle. Furthermore, considering that the higher the
with runner
respectwith
to therespect
bladetocongruent
the blade angle.
congruent angle. Furthermore,
Furthermore, consideringconsidering
that thethat the higher
higher the of
the number
number of blades, the lower the blade loading, the pressure gradient at the outlet between the
number
blades, of blades,
the lower the lower
thesuction
blade the the
loading, blade loading,gradient
the pressure gradient at the outlet between the
pressure and the side should bepressure
lower in the at the
case with outlet
8 blades between
than thewith
in the case pressure and the
6 blades.
pressure
suction side and
should thebesuction
lower side
in should
the case be lower
with 8 in the
blades case
thanwith
in 8
theblades
case than
with in
6 the case
blades. with
The 6 blades. ratio
p/p
The 𝑝⁄𝑝 ratio was plotted along an arc at the mean radius in correspondence of the vane outlet mean
The 𝑝⁄𝑝 ratio was plotted along an arc at the mean radius in correspondence of the vane outlet
was plotted along an arc at the mean radius in correspondence of the vane outlet (z = 0 mm).
(z = 0 mm).
(z = 0 mm).

(a) (b)
(a) (b)
Figure 6.
6. View View of
of the the arc
arcarc (a)
(a)(a)
andand the
the iso-clips (b) where the pressure and the velocity were evaluated.
Figure
Figure 6. View of the and theiso-clips
iso-clips(b)
(b) where thepressure
where the pressureand
and
thethe velocity
velocity were
were evaluated.
evaluated.

Figure
Figure 7 points
7 points out the
out that that the higher the number of blades, the lower the pressure difference
Figure 7 points out thathigher the number
the higher of blades,
the number the lower
of blades, the pressure
the lower difference
the pressure between
difference
between
the two the
sidesthe two
of two sides
the sides
blade.of the blade.
Moreover, Moreover, the
the maximum maximum
p/pmean𝑝 ⁄𝑝 value occurs for the impeller
between of the blade. Moreover, the maximum 𝑝⁄value
𝑝 occurs
value for the
occurs impeller
for the with 6
impeller
with 6 blades and all the maxima are placed roughly between 15 and 35% of each arc from the
blades and6all
with the maxima
blades are maxima
and all the placed roughly between
are placed roughly15between
and 35%15ofand
each35%
arcof
from
eachthe
arcpressure
from theto the
pressure to the suction side.
pressure
suction side. to the suction side.

Figure 7. Pressure trends vs. the fraction of each arc from pressure side to suction side along the mean
Figure 7. Pressure
Figure trends
7. Pressure vs.vs.
trends thethe
fraction
fractionofofeach
each arc from
frompressure
pressureside
side
to to suction
suction sideside along
along the mean
the mean
radius at the vane outlet of the three geometries.
radius
radius at theatvane
the vane outlet
outlet of the
of the three
three geometries.
geometries.

Table 3 highlights the increase of the pressure difference between pressure and suction side with
the decrease in the number of blades. These values were calculated from the previous Figure 7, as the
difference between the extreme ends of the arcs.
Water 2019, 11, x FOR PEER REVIEW 8 of 24
Water 2019, 11, 565 8 of 23
Table 3 highlights the increase of the pressure difference between pressure and suction side with
the decrease in the number of blades. These values were calculated from the previous Figure 7, as the
difference
Table 3. Pressure the extreme∆p,
betweendifference, ends of the arcs.
between pressure and suction side of the blades for the three
geometries with different number of blades and nq equal to 20.
Table 3. Pressure difference, ∆𝑝, between pressure and suction side of the blades for the three
geometries with different number of blades and nq equal
6 Blades to 20.
7 Blades 8 Blades
∆p (Pa) 10,500
6 Blades 9000 8 Blades 6050
7 Blades
σturbine Δp (Pa) 0.947 10500 0.951
9000 6050 0.953
σturbine 0.947 0.951 0.953
After that, the velocity profiles were plotted on a dimensionless axis indicating the fraction of each
After that, the velocity profiles were plotted on a dimensionless axis indicating the fraction of
arc from the pressure side to the suction side, where ∆θ = 360◦ /Zblades ; specifically, the velocities refer
each arc from the pressure side to the suction side, where ∆𝜃 = 360°/𝑍 ; specifically, the
to the blade congruent absolute tangential velocity, cu1,th . By comparing their trends, the maximum
velocities refer to the blade congruent absolute tangential velocity, 𝑐 , . By comparing their trends,
deflection occurs towards
the maximum deflection30% of the
occurs arc, and
towards 30% this increases
of the by decreasing
arc, and this increases by the number
decreasing theofnumber
blades, see
Figureof8.blades, see Figure 8.

FigureFigure 8. Absolute
8. Absolute tangential
tangential velocity
velocity calculatedvia
calculated via CFD
CFD (𝑐
(cu1))over
overthe blade
the congruent
blade value
congruent for the
value for the
three three geometries
geometries (𝑐 )., ).
(cu1,th

TheseThese results
results are coherent
are coherent with
with thethe trend
trend ofofthe
thestatic
staticpressure
pressure at
at the
the outlet
outletofofthe
thevane
vane(Figure
(Figure 7).
The 6 blades geometry shows steep variations close to the blades which leads to deflections upup
7). The 6 blades geometry shows steep variations close to the blades which leads to deflections to to
+100%
+100% compared to the blade congruent value; on the other hand, the 8 blades geometry is subjected
compared to the blade congruent value; on the other hand, the 8 blades geometry is subjected to
to smooth excursions and it displays a flat profile. To describe the flow deflection at the outlet of a
smooth excursions and it displays a flat profile. To describe the flow deflection at the outlet of a
centripetal hydraulic turbo machine, the slip factor, 𝜎 , has here been defined as the specific
centripetal hydraulicvia
work calculated turbo
CFD,machine,
𝑌 , over the
theslip factor, σspecific
theoretical turbine , has here
work, 𝑌 been defined as the specific work
,
calculated via CFD, YCFD , over the theoretical specific work, Y1D,th
𝑌
𝜎 = (3)
𝑌 ,
YCFD
σturbine = (3)
The results obtained by applying Equation (3) are Y1D,th
summed up in Table 3, where it is evident that
the slip factor decreases with the reduction of the number of blades. Hence, it can be stated that the
The results
higher obtained
the number by applying
of the blades, theEquation
better the (3)
floware summed
guidance at up
the in
exitTable
of the3,vanes
where it is evident
(Figure 8); this that
the slip factor decreases
guarantees with
lower losses andthe reduction
greater of the number
work extraction by the of blades. Hence, it can be stated that the
runner.
higher the number of the blades, the better the flow guidance at the exit of the vanes (Figure 8); this
guarantees lower losses and greater work extraction by the runner.
Moreover, the value of the absolute tangential velocity was evaluated on iso-clips (Figure 6b) by
means of the following expression
s
(cu u)ρc∆dA
cu1 = (4)
Gu1
Water
Water 2019, 11, x565
FOR PEER REVIEW 9 9of
of 24
23

where G is the ideal mass flow rate (kg/s) that flows through each iso-clip in Figure 6b. This calculation
allows the evaluation of the variation of the absolute tangential velocity moving tangentially from
pressure side (θ/∆θ = 0) to suction side (θ/∆θ = 1), see Figure 9. In this case, the velocity component
is weighted
Water 2019, 11,by thePEER
x FOR massREVIEW
flow rate, which actually flows through each iso-clip. 9 of 24

Figure 9. Trend of the absolute tangential velocity calculated on the iso-clips with Equation (4).
𝑐 , , = 6.03 m/s, 𝑐 , , = 4.77 m/s.

Moreover, the value of the absolute tangential velocity was evaluated on iso-clips (Figure 6b) by
means of the following expression

Trend of ∬ 𝑐 𝑢)𝜌𝑐̅ ∙ 𝑑𝐴 on the iso-clips with Equation (4).


Figure
Figure 9. 9. Trend of the
the absolute
absolute tangential velocity calculated
𝑐 =velocity calculated on the iso-clips with Equation (4). (4)
tangential
c𝑐u1,th,mean ==6.03 m/s,cu1,th,
6.03m/s, 𝑐 c line ==4.77
4.77m/s.
m/s.
𝐺𝑢
, , , ,
where 𝐺 is the ideal mass flow rate (kg/s) that flows through each iso-clip in Figure 6b. This
In this calculation, the results are influenced by the distribution of the flow rate, which is not
Moreover,
calculation the the
allows valueevaluation
of the absolute tangential
of the variationvelocity
of thewas evaluated
absolute on iso-clips
tangential (Figure
velocity 6b) by
moving
homogeneous. Indeed, due to the Coriolis force, F coriolis = 2mw × ω, the mass flow rate is thickened
means of thefrom
tangentially following side (𝜃⁄∆𝜃 = 0) to suction side (𝜃⁄∆𝜃 = 1), see Figure 9. In this case, the
expression
pressure
close to the hub.
velocity component is weighted by the mass flow ∬rate, which
𝑐 𝑢)𝜌𝑐̅ 𝑑𝐴actually flows through each iso-clip.
The meridional velocity over the theoretical value at the∙ outlet of the vane is displayed in Figure 10.
𝑐 = by the distribution of the flow rate, which is not
In this calculation, the results are influenced (4)
In conclusion, pressure gradients along the tangential𝐺𝑢 and radial direction and secondary flows due
homogeneous. Indeed, due to the Coriolis force, 𝐹 = 2𝑚𝑤 𝜔, the mass flow rate is thickened
to the Coriolis
where 𝐺 is effect
the occur
ideal at the
mass vane
flow exit.
rate These,
(kg/s) thattogether
flows with the each
through bladeiso-clip
loading,inpromote
Figure the
6b. flow
This
close to the hub.
deflection at the outlet of the vane.
calculation allows the evaluation of the variation of the absolute tangential velocity moving
tangentially from pressure side (𝜃⁄∆𝜃 = 0) to suction side (𝜃⁄∆𝜃 = 1), see Figure 9. In this case, the
velocity component is weighted by the mass flow rate, which actually flows through each iso-clip.
In this calculation, the results are influenced by the distribution of the flow rate, which is not
homogeneous. Indeed, due to the Coriolis force, 𝐹 = 2𝑚𝑤 𝜔, the mass flow rate is thickened
close to the hub.

Figure 10. Contours of the meridional velocity over the theoretical value and representation of the
Figure 10. Contours of the meridional velocity over the theoretical value and representation of the
Coriolis force and its effect at the outlet of the vane.
Coriolis force and its effect at the outlet of the vane.

The meridional velocity over the theoretical value at the outlet of the vane is displayed in Figure
10. In conclusion, pressure gradients along the tangential and radial direction and secondary flows

Figure 10. Contours of the meridional velocity over the theoretical value and representation of the
Coriolis force and its effect at the outlet of the vane.
Water 2019, 11, x FOR PEER REVIEW 10 of 24
Water 2019, 11, 565 10 of 23
due to the Coriolis effect occur at the vane exit. These, together with the blade loading, promote the
flow deflection at the outlet of the vane.
A local
A local reversed
reversed flowflow (negative
(negative velocity)
velocity) cancan be evidenced
be evidenced across
across thethe outlet
outlet of the
of the vane
vane andand
thethe
inlet of the draft tube, as depicted by the velocity contours on the meridional plane (see Figure 11).11).
inlet of the draft tube, as depicted by the velocity contours on the meridional plane (see Figure

Figure 11. Contours of the meridional velocity over the theoretical value on a meridional plane, across
the vane
Figure outlet and
11. Contours the meridional
of the initial part velocity
of the draft
overtube.
the theoretical value on a meridional plane, across
the vane outlet and the initial part of the draft tube.
4. Numerical Study of a Real PaT
After the
4. Numerical evaluation
Study ofPaT
of a Real the slip phenomenon at the outlet of the runners with no blade thickness,
numerical investigations of a real double suction centrifugal pump have been carried out with the
After the evaluation of the slip phenomenon at the outlet of the runners with no blade thickness,
open-source CFD code OpenFOAM by solving the 3D U-RANS equations. The unsteady RANS
numerical investigations of a real double suction centrifugal pump have been carried out with the
equations were considered adequate to model the flow through the pump, where quantities have to be
open-source CFD code OpenFOAM by solving the 3D U-RANS equations. The unsteady RANS
considered averaged over a time period short enough with respect to global unsteady phenomena
equations were considered adequate to model the flow through the pump, where quantities have to
but long enough for statistical significance. The transient simulations were run with an OpenFOAM
be considered averaged over a time period short enough with respect to global unsteady phenomena
application, pimpleDyMFoam, which is able to run transient simulations, including moving meshes,
but long enough for statistical significance. The transient simulations were run with an OpenFOAM
with incompressible flow. This application is based on the PIMPLE algorithm: a combination of
application, pimpleDyMFoam, which is able to run transient simulations, including moving meshes,
the PISO (Pressure Implicit with Splitting of Operator) and SIMPLE (Semi-Implicit Method for
with incompressible flow. This application is based on the PIMPLE algorithm: a combination of the
Pressure-Linked Equations). The turbulence model applied for the system closure is the k-ω SST
PISO (Pressure Implicit with Splitting of Operator) and SIMPLE (Semi-Implicit Method for Pressure-
by Menter et al. [26].
Linked Equations). The turbulence model applied for the system closure is the k-ω SST by Menter et
This procedure has been validated in the past against experimental data of a double suction
al. [26].
centrifugal pump [28].
This procedure has been validated in the past against experimental data of a double suction
centrifugal pumpDomain
4.1. Numerical [28]. and Boundary Conditions

The geometry
4.1. Numerical Domain of anda Boundary
commercial pump with a specific speed nq ∼
Conditions = 21 was simulated in both direct
(pump) and inverse mode (turbine). The numerical domain studied in this work is represented in
The geometry of a commercial pump with a specific speed 𝑛 ≅ 21 was simulated in both direct
Figure 12. It includes the suction pipe, the impeller and the discharge pipe. Herein, the impeller
(pump) and inverse
investigated mode (turbine).
is a double The numerical
suction impeller, actuallydomain studied
two single in this
suction work is represented
centrifugal pump impellers in
Figure 12. It includes the suction pipe, the impeller and the discharge pipe.
with seven blades in a back-to-back configuration. This kind of centrifugal pump is employed to Herein, the impeller
investigated
minimize is thea double suction
net positive impeller,
suction headactually
requiredtwo single
when suctionincentrifugal
operating pump mode. pump impellers
Moreover, thewith
volute
seven blades in
is doubled toabalance
back-to-back configuration.
the radial loads on theThisrotorkind
and of centrifugal
allow pump
high-speed is employed
operation, also attopart
minimize
load [18].
theThe
netcharacteristic
positive suction curve,head required
see Figure 13, when operating
of the pump was in pump mode.
calculated by meansMoreover, the volute
of unsteady is
simulations.
doubled to balance theand
The methodologies radial loads on
the results thebeen
have rotordeeply
and allow high-speed
described operation,
by Capurso et. al.also
[21].at part load
[18]. The characteristic curve, see Figure 13, of the pump was calculated
For all the analyzed cases, the mass flow rate has been imposed at the inlet considering by means of unsteady
a uniform
simulations. The methodologies and the results have been deeply described
inlet velocity distribution and a constant turbulent intensity equal to 3%. Moreover, by Capurso et.the
al. [21].
value of
For all the analyzed cases, the mass flow rate has been imposed at the inlet
the mass flow leakage, which flows through the annular seal, has been modeled as exiting from theconsidering a uniform
inlet velocitycase
impeller distribution and aaxially
and incoming constant turbulent
upstream theintensity
impellerequal to 3%.
eye with Moreover,
a 45 the value
◦ of swirl with of the
respect to the
mass flow leakage, which flows through the annular seal, has been modeled
tangential direction. The leakage was calculated a priori according to a one-dimensional empirical as exiting from the
impeller
modelcasebecauseand incoming
the geometryaxially upstream
of the seal hasthe notimpeller eye with
been included in athe
45° of swirl with domain
computational respect to the
[21].
Water 2019, 11, x FOR PEER REVIEW 11 of 24
Water 2019, 11, x FOR PEER REVIEW 11 of 24

tangential direction. The leakage was calculated a priori according to a one-dimensional empirical
tangential direction. The leakage was calculated a priori according to a one-dimensional empirical
model because the geometry of the seal has not been included in the computational domain [21].
Water 2019,
model 11, 565the geometry of the seal has not been included in the computational domain [21].
because 11 of 23

Suction pipe
Suction pipe

PUMP
PUMP
OPERATION
OPERATION

Discharge pipe
Discharge pipe
Double volute
Double volute

Figure 12.
Figure 12. Front
Front and
and side
side view
view of
of the
the CAD
CAD representing
representing the the centrifugal
centrifugal pump;
pump; the
the areas
areas of
of Flange
Flange 11
Figure 12. Front and side view 22of the CAD representing the centrifugal pump; the areas of Flange 1
and 22 were
and were equal
equal to
to 0.02840 m and
0.02840 m and 0.01941
0.01941 m
m 22, respectively.
, respectively.
and 2 were equal to 0.02840 m2 and 0.01941 m2, respectively.

Figure 13.
Figure 13. Characteristic
Characteristic curve
curve ofof the
the commercial
commercial pump
pump calculated
calculated by
by means
means of
of numerical
numerical simulations
simulations
Figure 13. Characteristic curve of the commercial pump calculated by means of numerical simulations
with CFX and 4040 million
million of
of cells.
cells.
with CFX and 40 million of cells.
A uniform pressure distribution was imposed at the outlet of the the domain.
domain. Straight pipes were
A uniform pressure distribution was imposed at the outlet of the domain. Straight pipes were
added atat the
the inlet
inletand
andthetheoutlet
outletto
toextend
extendthe
thenumerical
numericaldomain
domainininorder
ordertoto
reduce
reduceuncertainties due
uncertainties to
due
added at the inlet and the outlet to extend the numerical domain in order to reduce uncertainties due
boundary
to boundary conditions.
conditions. TheThe
three parts
three of the
parts geometry,
of the which
geometry, werewere
which merged together,
merged communicate
together, communicatewith
to boundary conditions. The three parts of the geometry, which were merged together, communicate
each
with other
each by means
other by of interfaces.
means In OpenFOAM,
of interfaces. the boundary
In OpenFOAM, condition used
the boundary on these
condition interfaces
used is
on these
with each other by means of interfaces. In OpenFOAM, the boundary condition used on these
named cyclicAMI
interfaces is named(i.e., Cyclic Arbitrary
cyclicAMI (i.e., CyclicMesh Interface).
Arbitrary Mesh Interface).
interfaces is named cyclicAMI (i.e., Cyclic Arbitrary Mesh Interface).
An important feature is the wall roughness; indeed, different values of the equivalent sand grain
An important feature is the wall roughness; indeed, −−5 different values of the equivalent−sand
5 m) grain
roughness
roughness have
havebeenbeenimposed
imposedtotothe thestator
stator(5.6
(5.6××1010 m) and
andthe
therotor
rotorsurfaces
surfaces(5.6 10−66 m). These
(5.6××10
roughness have been imposed to the stator (5.6 × 10−5 m) and the rotor surfaces (5.6 × 10−6 m). These
values are typical values measured on real impellers manufactured manufactured for experimental
experimental test [29]. Further
values are typical values measured on real impellers manufactured for experimental test [29]. Further
aboutgrid
details about grid refinement
refinement study study and boundary
and boundary conditions
conditions employedemployed for the simulations
for the numerical numerical
details about grid refinement study and boundary conditions employed for the numerical
simulations
are describedare in adescribed
previous in a previous
work [21]. Theworkresults[21]. The in
shown results shown inparagraphs
the following the following
wereparagraphs
calculated
simulations are described in a previous work [21]. The results shown in the following paragraphs
were
on calculated
a grid made of on11a million
grid made of represented
cells, 11 million cells, represented
in Figure 14. in Figure 14.
were calculated on a grid made of 11 million cells, represented in Figure 14.
Water 2019, 11, 565 12 of 23
Water 2019, 11, x FOR PEER REVIEW 12 of 24

Figure 14. Front and side view of the CAD representing the centrifugal pump; the areas of Flange 1
Figure 14. Front and side view of the CAD representing the centrifugal pump; the areas of Flange 1
and 2 were equal to 0.02840 m22 and 0.01941 m22 , respectively.
and 2 were equal to 0.02840 m and 0.01941 m , respectively.
4.2. Pump Mode Operation
4.2. Pump Mode Operation
The numerical set up has been assessed in a previous work [21] by comparing the results obtained
The numerical
with OpenFOAM andset
theup has been assessed
consolidated commercial in aCFD
previous
code CFX.workThese
[21] results
by comparing
show a fairlythe results
good
obtained with OpenFOAM and the consolidated commercial
agreement with each other, as already experienced by Nilsson [30]. CFD code CFX. These results show a
fairlyAll
good agreement with each other, as already experienced by Nilsson [30].
the transient flow simulations were carried out with a time ∆t = T/(256 Zblades ) = 4.5 × 10 s, − 6

where AllT the transient


is the flow of
time period simulations
a complete were carried
rotation andout with=a7time
Zblades is the∆tnumber
= T/(256ofZblades,
blades) = 4.5 × 10−6 s,
for a time
where T is the time period of a complete rotation and Z blades= 7 is the number of blades, for a time lapse
lapse equal to 0.085 s, which corresponds to five complete impeller revolutions. The results in terms of
equaland
head to 0.085 s, which
efficiency corresponds
were averaged overto five
thecomplete impeller revolutions.
last 3 revolutions. The mean and Thetheresults
maximumin terms of head
values of
and efficiency were averaged over the last 3 revolutions. The mean and the maximum
the Courant number have been approximately equal to 0.03 and 6, respectively, in all the simulations. values of the
Courant
The number have been
numerically approximately
obtained equalcurves
characteristic to 0.03 ofandthe6, respectively,
impeller working in all the simulations.
either in pump
The numerically obtained characteristic curves of the impeller
(see Figure 12) or in turbine mode (see Figure 15) are represented on the same diagram working either in pumpwith (see
Figure 12) or inaxes turbine mode (see Figure 15) are represented on the same diagram with dimensionless
dimensionless (H/H BEP vs Q/Q BEP ) in Figure 16, where HBEP and Q BEP are the head and the
axes (𝐻 ⁄ 𝐻 vs 𝑄 ⁄𝑄 ) in Figure 16, whereFor
flow rate at the BEP of the pump, respectively. 𝐻 completeness,
and 𝑄 are the head
HBEP,P /HBEP,Tandandthe Q flow rate at the
BEP,P /Q BEP,T
BEP of the pump, respectively. For
are equal to 0.958 and 0.832, respectively.completeness, 𝐻 , ⁄ 𝐻 , and 𝑄 , ⁄ 𝑄 , are equal to 0.958
and 0.832,
The head,respectively.
H (m), is evaluated as follows
The head, 𝐻 (m), is evaluated as follows
Ptot − Ptot
H = 𝑃 2 − 𝑃1 (5)
𝐻= ρg (5)
𝜌𝑔
2 ), the
whereP𝑃
where tot is is
thethe
total pressure
total (N/m
pressure (N/m subscripts
2), the 1 and1 2and
subscripts correspond to the inlet
2 correspond to theand outlet
inlet andflanges,
outlet
see Figure
flanges, see1,Figure
whereas the efficiency,
1, whereas η, is defined
the efficiency, 𝜂, isasdefined as

∆P ∆𝑃 Q𝑄
ηP𝜂= = 𝐶𝜔 (6)
tot
(6)

where the total pressure variation is 𝑃 −𝑃 , 𝐶 is the torque (Nm) resulting from the forces
where the total pressure variation is Ptot,2, − Ptot,1, , C is the torque (Nm) resulting from the forces
acting on the blades, hub and shroud of the impeller, comprising the disk friction losses, and 𝜔 is
acting on the blades, hub and shroud of the impeller, comprising the disk friction losses, and ω is the
the rotational speed (rad/s).
rotational speed (rad/s).
4.3. Turbine
4.3. Turbine Mode
Mode Operation
Operation
Eventually,aacampaign
Eventually, campaign ofof simulations
simulations waswas carried
carried out out
withwith the machine
the machine operated
operated in reverse
in reverse mode
mode applying the same setup previously described for the pump mode. However, the
applying the same setup previously described for the pump mode. However, the flow direction andflow direction
androtation
the the rotation
of theofimpeller
the impeller
werewere inverted,
inverted, see Figure
see Figure 15. 15.
In this case, the efficiency, 𝜂 , is defined
In this case, the efficiency, ηT , is defined as as
𝐶𝜔
𝜂 =Cω (7)
ηT = ∆𝑃 𝑄 (7)
∆Ptot Q
Water 2019, 11, x FOR PEER REVIEW 13 of 24
Water
Water 11,11,
2019,
2019, 565 PEER REVIEW
x FOR 13 13
of of
2423

Discharge pipe

Discharge pipe

TURBINE
OPERATION
TURBINE
OPERATION
Suction pipe

Double volute Suction pipe

Double volute

Figure 15. Numerical domain with the turbine operating conditions.


Figure
Figure 15.15. Numerical
Numerical domain
domain with
with thethe turbine
turbine operating
operating conditions.
conditions.
It is worth highlighting that no geometry modifications were applied to the impeller working in
It
reverse ismode,
worthashighlighting
usually happensthat no ingeometry modifications
the reality, where the were were applied
trailing edge of to the impeller working in
It is worth highlighting that no geometry modifications applied to the
the blades
impeller areworking
rounded into
reverse
reduce mode,
the flowas usually
incidence happens
losses. in the reality,
Generally, where
pumps the trailing
operated as edge
turbines of the
workblades
with are
a rounded
higher massto
reverse mode, as usually happens in the reality, where the trailing edge of the blades are rounded to
reduce
flow rate the flow incidence
thanincidence
when they losses. Generally,
are Generally, pumps
operated inpumps operated
directoperated as turbines
mode. Therefore, work with a higher mass flow
reduce the flow losses. as turbinesnumerical
work withsimulations
a higher mass were
rate than when they are operated in direct mode. Therefore, numerical simulations were performed
performed
flow rate than starting they 𝑄are ⁄
when from 𝑄 , equal
operated to 70%,
in direct mode.which is a value
Therefore, close to the
numerical runawaywere
simulations curve
starting
calculated from Qmin /Q
according to the
BEP,P equal to model
empirical 70%, which is a by
proposed value close
Gülich [18]to(see
theFigure
runaway 16a),curve
to 𝑄 calculated
⁄𝑄
performed starting from 𝑄 ⁄𝑄 , equal to 70%, which is a value close to the runaway curve ,
according
equal to 140%. to theThe
empirical
maximum modelflow proposed
rate allowsby Gülich
us toby [18] (see
identify Figure 16a),
the to Qmax /Q BEP,P equal to
calculated according to the empirical model proposed Gülichthe BEP
[18] (seeofFigure turbine (see
16a), to 𝑄 Figure
⁄𝑄 16b).
,
140%.
The to
massThe maximum
flow leakagesflow rate
wereflow allows
calculated us to identify
by means the BEP of the turbine (see Figure 16b). The mass
equal 140%. The maximum rate allows us to of the same
identify theequations
BEP of theprovided by Gülich
turbine (see Figure[18]
16b).for
flow leakages were calculated by means of the
the volumetric
same equations provided by Gülich [18] for constant
centrifugal
centrifugal pumps. Applying this model, efficiency (𝜂
The mass flow leakages were calculated by means of the same equations provided by Gülich [18] forat
) for a PaT results
pumps.
various flow Applying this model, the volumetric efficiency (ηv ) for a PaT results constant at various
rate. Applying
centrifugal pumps. this model, the volumetric efficiency (𝜂 ) for a PaT results constant at
flow rate.
various flow rate.

(b)
(a)
(b)
(a) hydraulic turbo machine performance in both pump and turbine operating
Plot of the
Figure 16. Plot
modes: (a) head curves calculated with Equation (5)(5)
with runaway curve; (b) Efficiency calculated with
Figure 16.(a) head
Plot curves
of the calculated
hydraulic with
turbo Equation
machine with
performance runaway curve;
in both pump (b) Efficiency
and calculated
turbine operating
Equations
with (6) and
Equations (6)(7).
and (7).
modes: (a) head curves calculated with Equation (5) with runaway curve; (b) Efficiency calculated
4.4.with Equations
Evaluation (6) and
of the (7). at the Outlet of the Runner
Velocity
4.4. Evaluation of the Velocity at the Outlet of the Runner
Actually,ofthe
4.4. Evaluation the main flow
Velocity is subject
at the Outlet ofto a deflection in the same direction of the machine rotation,
Actually, the main flow is subject tothe Runner in
a deflection the same direction of the machine rotation,
which could be justified by different sources as explained in Sections 2 and 3, e.g., the vorticity inside
which could the be justified by different
subjectsources as explainedtheinto Sections 2 andof 3, the
e.g.,machine
the vorticity inside
theActually,
channel and the main flow
pressure isgradient. to a deflection
This deflectionin
leads same direction
a reduction of the rotation,
work extracted by the
the
which channel
could and the
be justifiedpressure gradient.
by different This
sources as deflection
explained leads
in Sections 2 and 3, e.g., the vorticity insideby
to a reduction of the work extracted
runner according to specific work equation (Equation (1)).
the
the runnerand
channel according
the of to specific
pressure work This
gradient. equation (Equation
deflection leads(1)).to
With the aim quantifying this phenomenon, using thea results
reduction of the
of the work extracted
numerical by
simulations,
the With the aim
runner according of quantifying this phenomenon, using the results of the numerical simulations,
absolute tangential to andspecific
axial work equation
velocity (Equation
components were (1)).evaluated in Section 1, see Figure 1, by
absolute
With tangential
the aim of and axial velocity
quantifying this components were
phenomenon, using evaluated
the results in Section 1, see Figure
of the exit.
numerical 1, by means
means of area-weighted averages on surfaces downstream the impeller Thesesimulations,
values were
of area-weighted
absolute tangential averages
and axial on surfaces downstream the impeller exit. These values1, were time-
time-averaged over the lastvelocity components
3 rotations were evaluated
of the machine. The mean in Section
value 1,ofsee
theFigure
absolute by means
tangential
averaged
ofvelocity over
area-weighted the last
averages3 rotations of
on surfaces the machine.
downstream The mean value
the impeller of the absolute
exit. These tangential
values velocity
were time-
is 4.40 m/s. This value has been compared with the theoretical value [21] computed at the
averaged over the last 3 rotations of the machine. The mean value of the absolute tangential velocity
Water 2019, 11, x FOR PEER REVIEW 14 of 24
Water 2019, 11, x FOR PEER REVIEW 14 of 24

is 4.40 m/s. This value has been compared with the theoretical value [21] computed at the same flow
is 4.40
Water
rate,
2019,m/s. This with
11, 565
congruently valuethehasblade
beenangle
comparedat thewith theoftheoretical
outlet value
the impeller (𝛽 [21]
). Forcomputed
the machine at the same
14 offlow
studied
23
in
rate,
this congruently
work, with the
the theoretical blade angle
tangential at the
velocity, outlet of by
calculated themeans
impeller (𝛽 blade
of the ). Forcongruent
the machine flowstudied
theoryin
atthis work,
the BEP the theoretical tangential velocity, calculated by means of the blade congruent flow theory
same flowof thecongruently
rate, turbine, is equal
withtothe
2.08 m/s.angle
blade Furthermore, the presence
at the outlet of the blades
of the impeller (β 1B ). induces a strong
For the machine
at
and the BEP
localized of the turbine, is equal to 2.08 m/s. Furthermore, the presence of the blades induces a strong
studied in this deflection of the flowtangential
work, the theoretical field. It was noted
velocity, that this by
calculated might
meansbe ofjustified
the bladebycongruent
the tangential
flow
and localized
pressure gradient deflection
that of the
occurs flow field.
between It was
pressure andnoted that
suction thisofmight
side the be justified
blades (see by the
Figure tangential
17). This is
theory at the BEP of the turbine, is equal to 2.08 m/s. Furthermore, the presence of the blades induces
pressure
also shown gradient
in Figure that
18, occurs between
displaying the pressurealong
pressure and suction
an arc side
at theofmean
the blades
radius.(see Figure 17). This is
a strong and localized deflection of the flow field. It was noted that this might be justified by the
also shown in Figure 18, displaying the pressure along an arc at the mean radius.
tangential pressure gradient that occurs between pressure and suction side of the blades (see Figure 17).
This is also shown in Figure 18, displaying the pressure along an arc at the mean radius.

Figure
Figure 17.
17. Contours
Contours of the static
of the static pressure
pressure (N/m
(N/m22))evidencing the
the pressure
pressure gradient
gradient across
across the
the blade
Figure 17. Contours of the static pressure (N/m 2evidencing
) evidencing the pressure gradient across
blade
the blade
trailing edge and at the mean radius.
trailing edge and at the mean radius.
trailing edge and at the mean radius.

Figure 18. Pressure over the mean pressure (p/p ) along the mean radius at the trailing edge of the
Figure 18. Pressure over the mean pressure (p/pmean
mean) along the mean radius at the trailing edge of the
blade at the
Figure 18. BEP.
Pressure over the mean pressure (p/p mean) along the mean radius at the trailing edge of the
blade at the BEP.
blade at the BEP.
To point out this contribution, axial and tangential velocity components were evaluated along the
To point out this contribution, axial and tangential velocity components were evaluated along
mean radius of the
To radius
point out blades close to theiraxial
thisblades
contribution, trailing edge. The curve trends show a periodic behavior along
due
the mean of the close to their and tangential
trailing velocity
edge. The curvecomponents
trends showwere evaluated
a periodic behavior
tothe
themean
presence of the blades
radius ofofthe (see Figure 19).
due to the presence theblades
bladesclose
(see to their 19).
Figure trailing edge. The curve trends show a periodic behavior
This
dueThis means
to the that the
presence of main
the flow (see
blades at the outlet19).
Figure of the turbine is subject to a deflection not only due to
means that the main flow at the outlet of the turbine is subject to a deflection not only due
the vortex inside the channel but also to the pressureof gradient acting in the same direction.not
To only
further
to the This means
vortex insidethat
thethe main flow
channel at thetooutlet
but also the turbine
the pressure is subject
gradient acting to
in athe
deflection
same direction. due
To
analyze
to the vortex inside the channel but also to the pressure gradient acting in the same direction.toTo
the velocity vector components let’s have a look at a short arc of the curve close to a blade
evaluate the tangential velocity component downstream of the runner.
Water 2019, 11, x FOR PEER REVIEW 15 of 24

Water 2019, 11, x FOR PEER REVIEW 15 of 24


further analyze the velocity vector components let’s have a look at a short arc of the curve close to a
blade to evaluate the tangential velocity component downstream of the runner.
further analyze the velocity vector components let’s have a look at a short arc of the curve close to a
Water 2019, 11, 565 15 of 23
blade to evaluate the tangential velocity component downstream of the runner.

Figure 19. Representation of the axial and tangential components of the absolute velocity vector along
Figure Representation
19. radius
the mean of theedge
at the trailing axialofand
thetangential
blade (BEPcomponents of the absolute velocity vector along
operating point).
Figure 19. Representation of the axial and tangential components of the absolute velocity vector along
the mean radius at the trailing edge of the blade (BEP operating point).
the mean radius at the trailing edge of the blade (BEP operating point).
The tangential velocity was plotted on the same curve at different axial distance from the trailing
edgeThe tangential
(closer velocity
the curve, higherwastheplotted onproximity
peak). In the same to
curve at different
the blade, there isaxial distance
a great from
increase of the
the
The
trailing tangential
edge velocity
(closer the was
curve, plotted
higher on the same
the gradient curve
peak). Incomparedat different
proximitytotothe axial
themain
blade,distance
there from
is a greatthe trailing
increase
tangential component due to the pressure flow, whereas towards the
edge
of the(closer the curve,
tangential higher due
component the peak).
to the In proximity
pressure to thecompared
gradient blade, thereto is a great
the main increase
flow, of the
whereas
discharge, the velocity becomes more uniform and the peak gradually becomes smoother (see Figure
tangential
towards thecomponent
discharge, due to the pressure
the velocity becomes gradient compared
more uniform andto the
the main
peak flow, whereas
gradually becomes towards the
smoother
20).
discharge,
(see Figurethe20).velocity becomes more uniform and the peak gradually becomes smoother (see Figure
20).

(a) (b)

Figure 20.
Figure 20. Local
Local analysis
analysis
(a) of ofthe
theabsolute
absolutetangential
tangential velocity
velocity distribution: (a)(b)
distribution: (a) representation of
representation of an
an arc
arc
and three
and three different
different distances
distances from
fromthe
thetrailing
trailingedge
edgeofofthe
theblade
blade(0,
(0,+a,
+a,+2a);
+2a); (b)
(b) absolute
absolute tangential
tangential
Figure 20. Local analysis of the absolute tangential velocity distribution: (a) representation of an arc
velocity
velocityalong
alongan anarc
arcat
atthree
threedifferent
differentdistances
distancesatatthe
themean
meanradius
radiusatatthe
thetrailing
trailingedge
edgeofofthe
the blade
blade at
at
and three different distances from the trailing edge of the blade (0, +a, +2a); (b) absolute tangential
the
the BEP.
BEP.
velocity along an arc at three different distances at the mean radius at the trailing edge of the blade at
the BEP.
The distances “a” = 0.5 cm and “2a” = 1 cm were chosen to be short enough to stay far from the
The distances “a“= 0.5 cm and “2a” = 1 cm were chosen to be short enough to stay far from the
influence
influenceof
of the
the curved
curved suction
suction elbow,
elbow,but
butsufficiently
sufficientlyclose
closeto
tothe
theleading
leadingedge
edgeto
to investigate
investigate the
the area
area
The distances “a“= 0.5 cm and “2a” = 1 cm were chosen to be short enough to stay far from the
affected by the flow deflection.
affected by the flow deflection.
influence of the curved suction elbow, but sufficiently close to the leading edge to investigate the area
As previously described, the slip factor gathers different sources of deflection, e.g., counter
affected by the flow deflection.
rotating vortex and pressure gradient at the trailing edge of the blade. Applying Euler’s equation
in order to compute the work both via 1-D theoretical model, Y1D,th , and CFD calculi, YCFD , at the
turbine BEP, the σturbine is equal to 0.967.
Water 2019, 11, x FOR PEER REVIEW 16 of 24

As previously described, the slip factor gathers different sources of deflection, e.g., counter
Water 2019,vortex
rotating 11, 565 and pressure gradient at the trailing edge of the blade. Applying Euler’s equation 16 of in
23

order to compute the work both via 1-D theoretical model, 𝑌 , , and CFD calculi, 𝑌 , at the
turbine BEP, the
Similarly 𝜎 theory
to the is equal to 0.967.
of pumps where the slip factor corrects the ideal work transferred by the
machine to the flow, the slip factor herewhere
Similarly to the theory of pumps the slip
introduced factor
might becorrects
appliedthe ideal
under work operating
turbine transferred by the
mode in
machine to the flow, the slip factor here introduced
order to correct the ideal work extracted by the PaT. might be applied under turbine operating mode
in order
Theto correct
slip factorthe
was ideal work extracted
computed by operating
at different the PaT. points, i.e., different mass flow rates, and the
The slip factor was
results are shown in Figure 21. computed at different operating points, i.e., different mass flow rates, and
the results
Then, are shown
a model of in
theFigure 21. was obtained by least-squares fitting, stated as follows
slip factor
Then, a model of the slip factor was obtained by least-squares fitting, stated as follows
 2  
Q Q
σturbine = 0.2365 𝑄 − 0.5537 𝑄 + 1.2846 (8)
σ = 0.2365Q BEP,P − 0.5537 Q BEP,P 1.2846 (8)
𝑄 , 𝑄 ,
Numerical simulations
Numerical simulations at
at different
different flow
flow rates
rates are
are suggested
suggested for
for the
the purpose
purpose of
of further
further exploring
exploring
the slip factor behavior.
the slip factor behavior.

Figure
Figure 21.
21. Slip
Slip factor
factor in
in turbine mode, σ
turbine mode, σturbine ,, vs
vsflow
flowrate
rateratio, 𝑄 ⁄𝑄BEP,T
ratio,Q/Q , ..

5. Theoretical
5. Theoretical 1-D
1-D Prediction
Prediction Model
Model
The results
The results of of numerical
numerical investigations
investigations were were applied
applied to to aa new
new 1-D
1-D model,
model, which
which was
was previously
previously
developed by
developed by Stefanizzi
Stefanizzi et et al.
al. [22].
[22]. TheThe work
work aims
aims at at adding
adding thethe slip
slip phenomenon
phenomenon effect
effect in
in order
order to
to
improve the model in terms of PaT performance prediction. The model is
improve the model in terms of PaT performance prediction. The model is accurately described in [22]accurately described in [22]
in order
in order toto predict
predict thethe entire
entire characteristic
characteristic of of aa PaT,
PaT, starting
starting from
from the
the knowledge
knowledge of of its
its geometrical
geometrical
parameters and
parameters and available
available related
related tests
tests inin turbine
turbine mode
mode operation.
operation.
Figure 22
Figure 22 schematically illustrates how
schematically illustrates how the
the model
model works: thanks to
works: thanks to the
the knowledge
knowledge of of detailed
detailed
geometrical data,
geometrical data,flow rate,𝑄,Q,
flowrate, and and rotational
rotational speed, 𝑛, itn,
speed, it is possible
is possible to accurately
to accurately calculate
calculate the
the correct
correct velocity
velocity triangles triangles
and the andtheoretical
the theoretical 𝐻 ,Hin
head,head, th , reverse
in reverse modeoperation.
mode operation.Afterwards,
Afterwards, volute
volute
(Zvolute,tot
(𝑍 ,
)) and
and runner
runner (Z (𝑍 runner,tot
,
)) losses
losses are
are modelled,
modelled, whereas
whereas the the theoretical
theoretical head
head isis reduced
reduced by by
means of
means of the
the slip
slip factor,
factor, σ𝜎turbine , ,proposed
proposedin inthis
this paper.
paper. All All these
these contributions
contributions are
are summed
summed up up to
to
finally predict the real PaT head, 𝐻turbine .
finally predict the real PaT head, H .
Some researchers,
Some researchers, suchsuch as as Barbarelli
Barbarelli et et al.
al. [17]
[17] and
and Gülich
Gülich [18],
[18], have
have developed theoretical
developed theoretical
approaches in order to predict PaT performance by taking into account velocity triangles, hydraulic
losses on simplified geometries rather than than using
using statistical
statistical and
and experimental
experimental correlations.
correlations. However,
However,
simplification of the geometry of impeller channel or volute could determine significant errors
the simplification
(±20%)in
(±20%) interms
termsof ofperformance
performanceprediction.
prediction.
WaterWater
2019,2019, 11, x FOR PEER REVIEW
11, 565 17 of1723of 24
Water 2019, 11, x FOR PEER REVIEW 17 of 24

Figure 22. Flow chart of the 1-D proposed model.


Figure
Figure 22. Flow chart
22. Flow chart of
of the
the 1-D
1-D proposed
proposed model.
model.
5.1. 5.1.
Volute Volute Losses
Losses
5.1. Volute Losses
For the For evaluation
the evaluation of volute
of volute losses,
losses, a generic
a generic doubledouble volute
volute constituted
constituted by inner
by the the inner
and and the outer
the outer
For
volute the evaluation
has been of volute
considered losses, a generic double volute constituted by the inner and the outer
volute has been considered and and divided
divided into into
threethree
parts, parts, as depicted
as depicted in Figure
in Figure 23: inlet
23: the the inlet convergent
convergent
volute has been
channel (part considered
1), half and divided
volute into(part
collector three2)parts,
and as depicted
the vane-less in Figure
space 23: the inlet
between theconvergent
throat volute
channel (part 1), half volute collector (part 2) and the vane-less space between the throat volute section
channel (part 1), half volute collector (part 2) and the vane-less space between the throat volute
andsection
the runner and inlet
the runner
sectioninlet
(partsection (part 3).
3). As stated in As stated (9),
Equation in Equation (9), thehead
the total volute total loss,
volute head loss,
Zvolute,tot ,
section
𝑍 and ,
the
, is runner inlet
obtained by section (part
summing up 3). As stated
different in Equation
contributes: (1) (9),
the the loss
head totalinvolute
the head
inlet loss,
convergent
is obtained by summing up different contributes: (1) the head loss in the inlet convergent nozzle
𝑍 nozzle , is(𝑍obtained by summing
, ); (2) theinvolute uploss
different
in(Zthecontributes:
collector (𝑍(1) the, );the head loss
(3)vanelessin the
the loss ininlet
the convergent
vaneless
,
(Zvolute,1 ); (2) the volute loss the collector volute,2 ); (3) the loss in space (Zvolute,3 ). space
nozzle
(𝑍 (𝑍 , ). ,
); (2) the volute loss in the collector (𝑍 , ); (3) the loss in the vaneless space
(𝑍 , ). = Zvolute,a
Zvolute,tot
𝑍 = 𝑍 = Zvolute,b
, = 𝑍 = Zvolute,1
, = 𝑍 + Zvolute,2𝑍+ Z volute,3𝑍.
, , , . (9) (9)
,
𝑍 , =𝑍 , =𝑍 , =𝑍 , 𝑍 , 𝑍 , . (9)

(b)
(a)
(b)
(a)
Figure 23. Division of the double volute: inner volute (a) and outer volute (b).
Figure 23. Division of the double volute: inner volute (a) and outer volute (b).
5.2. Runner Figure
Losses 23. Division of the double volute: inner volute (a) and outer volute (b).
5.2.
As Runner
stated in Losses
Equation (10), runner losses have different contributions. The first, Zrunner,1 , is
5.2. Runner Losses
connected to the friction
As stated in Equation loss inside the impeller
(10), runner losseschannel.
have different Furthermore,
contributions.with respect to Gülich’s
The first, 𝑍 , , is
approach,
As stated
connected the current
in the
to model(10),
Equation
friction introducesinsidenew
lossrunner the hydraulic
losses parameters,
have different
impeller channel. as depicted
contributions.
Furthermore, Thein Figure
with first, 𝑍24. to
respect A Gülich’s
,flow
, is
connected
incidence
approach, to the
loss is
the friction
contemplated loss at
current model inside
the the section
inlet
introduces impeller
new of channel.
the
hydraulic runner, Zrunner,2 , as
Furthermore,
parameters, with to
in depicted
order respect
take to Gülich’s
into
in Figure account
24. A flow
approach,
the difference
incidence theloss
current
between model introduces
the inlet
is contemplated bladeat the new section
congruent
inlet hydraulic
angle, parameters,
ofβthe2,B ,runner,
and the 𝑍flow as depicted
angle,
, , in β in.
order
2 Figure
Moreover,
to take24. A
a
into flow
blade
account
incidence
blockage loss
factor is contemplated
is used in order at the
to inlet
consider section
the of the
contractionrunner, of 𝑍
the
the difference between the inlet blade congruent angle, 𝛽 , , and the flow angle, 𝛽 . Moreover, inlet, , in order
cross-area. to take
This into account
modifies the a
the blade
inletdifference
velocity between
triangle, the
causing inlet
a blade
loss due congruent
to the angle,
variation
blockage factor is used in order to consider the contraction of𝛽 the
, , and
radial the flow
velocity angle,
component,𝛽
of the inlet cross-area. This . Z
Moreover, .a
modifies
runner,3
bladetheblockage factor istriangle,
inlet velocity used in order
causing to consider
a loss due thetocontraction
the variation of theofinlet
the cross-area.
radial velocityThis modifies
component,
the 𝑍inlet velocity triangle, Z
causing a
runner,tot =
loss Z due to
runner,1 + Z
the variation
runner,2 + Z of the
runner,3 radial velocity component, (10)
, .
𝑍 , .
𝑍 , =𝑍 , 𝑍 , 𝑍 , (10)
𝑍 , = 𝑍 , 𝑍 , 𝑍 , (10)
Water 2019,
Water11, 565
2019, 11, x FOR PEER REVIEW 18 of 24 18 of 23
Water 2019, 11, x FOR PEER REVIEW 18 of 24

(a)
(b)
(a)
(b)
24. Runner
FigureFigure losses:
24. Runner losses:(a)(a)the
theflow
flow incidence loss;
incidence loss; (b)(b) variation
variation ofradial
of the the radial component
component of the inlet
of the inlet
relative velocity.
relative
Figure velocity.
24. Runner losses: (a) the flow incidence loss; (b) variation of the radial component of the inlet
relative velocity.
5.3. Application of the
5.3. Application σturbineLeast-Square
σturbine
of the Fitting
Least-Square Fitting Curve
Curve to ato a 1-D
1-D ModelModel
5.3. Application of the σturbine Least-Square Fitting Curve to a 1-D Model
The 1-D
The model was
1-D model applied
was appliedto tothe
the machine numerically
machine numerically investigated
investigated in thein the previous
previous section 4.Section 4.
FigureFigure The
25 shows 1-Dthe
25 showsmodel
the was appliedbetween
comparison
comparison to the machine
between the numerically
theproposed
proposed investigated
model,
model, Gülich’s in themodel
model
Gülich’s previous section
and theand the4.current
current
Figure 25 shows
manufacturer’s the used
model comparison between the
in the prediction proposed
of the model,curve
characteristic Gülich’s model
during and mode.
turbine the current
All
manufacturer’s model used in the prediction of the characteristic curve during turbine mode.
manufacturer’s
predicted model
curves were used in the
compared withprediction
respect toofthat
theobtained
characteristic curveofduring
by means turbine
a numerical mode. All
simulation
All predicted curves were compared with respect to that obtained by means of a numerical simulation
predictedby
performed curves wereetcompared
Capurso al. [21]. with respect to that obtained by means of a numerical simulation
performed by Capurso
performed et al.et[21].
by Capurso al. [21].

Figure 25. Comparison of different prediction models for the PaT characteristic curve.
FigureFigure 25. Comparison
25. Comparison of differentprediction
of different prediction models
modelsfor thethe
for PaT characteristic
PaT curve.curve.
characteristic
As a result, the use of detailed geometrical information, the introduction of the slip factor and a
new
As As a result,
modelling
a result, the
theofuse use
theof of detailed
hydraulic
detailed geometrical
losses information,
have determined
geometrical thethe
a more
information, introduction
accurate of theof
slip
prediction
introduction offactor
the the and
slipPaT a and
factor
new modelling
performance under of the
design hydraulic
and losses
off-design have determined
conditions. Indeed, a
the more accurate
proposed
a new modelling of the hydraulic losses have determined a more accurate prediction of the PaT1-D prediction
model shows of the
the PaT
best
performance
prediction under
at the lowestdesign andhighest
and the off-design
flowconditions. Indeed, equal
rates (respectively the proposed
to −2.4%1-D
andmodel
0.1%)shows the best
with respect
performance under design and off-design conditions. Indeed, the proposed 1-D model shows the best
toprediction at the(respectively
Gülich’s model lowest and the highest
equal to 5%flow
andrates (respectively
18.1%) and to the equal to −2.4% and
manufacturer’s 0.1%)
model with respect
(respectively
prediction at
to the
2.4%lowest
to Gülich’s
equal model andasthe highest
(respectively
and 11%), equal
indicated inflow
to 5% rates
Table and (respectively
4. 18.1%) and to theequal to −2.4%model
manufacturer’s and 0.1%) with respect
(respectively
to Gülich’s
equalmodel
to 2.4%(respectively equal toin5%
and 11%), as indicated and4.18.1%) and to the manufacturer’s model (respectively
Table
equal to 2.4% and 11%), as indicated in Table 4.

Table 4. Comparison of different prediction models for PaT characteristic curve in terms of head
prediction errors.

Q/QBEP,T 0.67 0.82 1 1.17


Gülich 5.0% 17.5% 19.7% 18.1%
Proposed Model −2.4% 2.6% −1.5% 0.1%
Manufacturer 2.4% 4.8% 7.0% 11.0%

Figure 26 shows the contribution of each hydraulic loss considered in the proposed model in order
to predict the characteristic curve obtained by CFD. The introduction of volute and runner losses causes
an increasing of the turbine head, whereas the slip factor decreases the theoretical work calculated
Table 4. Comparison of different prediction models for PaT characteristic curve in terms of head
prediction errors. Q/QBEP,T 0.67 0.82 1 1.17
Gülich 5.0% 17.5% 19.7% 18.1%
Q/QBEP,T 0.67 0.82 1 1.17
Proposed Model −2.4% 2.6% −1.5% 0.1%
Gülich
Manufacturer 5.0%
2.4% 17.5%
4.8% 19.7%
7.0% 18.1%
11.0%
Proposed Model −2.4% 2.6% −1.5% 0.1%
Water 2019, 11, 565 Manufacturer 2.4% 4.8% 7.0% 11.0% 19 of 23
Figure 26 shows the contribution of each hydraulic loss considered in the proposed model in
order to predict the characteristic curve obtained by CFD. The introduction of volute and runner
Figure 26 shows the contribution of each hydraulic loss considered in the proposed model in
by meanslosses causes anequation,
of to
Euler’s increasing and
of thethus
turbine head, whereas
involves a betterthe slip factor for
prediction decreases the theoretical
the different work
operating
order predict the characteristic curve obtained by CFD. The introduction of volute and runner points.
calculated by means of Euler’s equation, and thus involves a better prediction for the different
Furthermore, the flow
losses causes incidence
an increasing loss
of the and the
turbine variation
head, of slip
whereas the the factor
radialdecreases
component of the relative
the theoretical work inlet
operating points. Furthermore, the flow incidence loss and the variation of the radial component of
calculated
velocitythe
show by
their means of
remarkable Euler’s equation,
contribution and thus
at lower involves
flow rate. a better prediction for the different
relative inlet velocity show their remarkable contribution at lower flow rate.
operating points. Furthermore, the flow incidence loss and the variation of the radial component of
the relative inlet velocity show their remarkable contribution at lower flow rate.

FigureFigure 26. Comparison


26. Comparison of differentsources
of different sources of
of losses
lossesfor
forthethe
performance prediction.
performance prediction.
Figure 26. Comparison of different sources of losses for the performance prediction.
Moreover,
Moreover, threethree different
different slipslip factordefinitions
factor definitions were
wereapplied
applied to tothethe
1-D1-D
model, as depicted
model, in
as depicted in
Figure
Figure 27. The 27. The results
resultsthree
show show that
that the the model with the assumption that the slip factor term is negligible
Moreover, different slipmodel
factor with the assumption
definitions were appliedthat
to the the1-D
slip factor
model, asterm is negligible
depicted in
(𝜎 = 1) overestimate the turbine head at the BEP with an error of 6.4%. On the other hand, both
= 1) overestimate
(σturbineFigure the turbine
27. The results show that the head
modelat thethe
with BEP with an that
assumption errortheofslip
6.4%. Onterm
factor theisother hand, both
negligible
the Busemann [23] and Stodola [24] formulations applied to the 1-D model provide good results at
(𝜎
the Busemann = 1) overestimate the turbine head at the BEP with an error of 6.4%. On the other hand, both
the BEP, [23] and Stodola the
but underestimate [24]turbine
formulations applied to
head at part-loads, thean1-D
with model
error provide
of −8%. Instead,good results at the
by applying
BEP, butthe Busemann [23] and Stodola [24] formulations applied to the 1-D model provide good results at
underestimate the turbine head at part-loads, with an error of − 8%. Instead,
the least-squares fitting of the CFD results, as described in Figures 21 and 26, errors are minimized by applying the
the BEP, but underestimate the turbine head at part-loads, with an error of −8%. Instead, by applying
least-squares fitting of the CFD results, as described in Figures 21 and 26, errors are minimized over a
over a wide range of flow rates.
the least-squares fitting of the CFD results, as described in Figures 21 and 26, errors are minimized
wide range
over aofwide
flowrange
rates.
of flow rates.

Figure 27. Comparison of different slip factor definitions (i.e., σturbine , σ proposed by Busemann [23], σ
proposed by Stodola [24] and σ = 1) with the CFD results.

6. Experimental Case Study


To further assess the importance of including the slip phenomenon at the outlet of the runner,
the model was applied to another double suction centrifugal pump, tested in both modes (nq = 26.46,
D2 = 432 mm). In this case, the machine is characterized by a similar specific speed number, nq ,
by the same number of blades and by similar outlet blade angle to those of the machine considered
in Section 4. This means that the least-square fitting curve of (Equation (8)) can be also used for
this machine. The proposed 1-D model has been applied to the machine by either considering or
neglecting the slip factor. For the former case, two slip formulations have been considered: the slip
factor evaluated by means of Equation (8) and the slip factor evaluated by Busemann. This comparison
6. Experimental Case Study
To further assess the importance of including the slip phenomenon at the outlet of the runner,
the model was applied to another double suction centrifugal pump, tested in both modes (𝑛 = 26.46,
𝐷 = 432 mm). In this case, the machine is characterized by a similar specific speed number, 𝑛 , by
the same number of blades and by similar outlet blade angle to those of the machine considered in
Water 2019, 11, 565 20 of 23
section 4. This means that the least-square fitting curve of (Equation (8)) can be also used for this
machine. The proposed 1-D model has been applied to the machine by either considering or
neglecting the slip factor. For the former case, two slip formulations have been considered: the slip
is illustrated in Figure 28, where the predicted curves are compared with respect to the experimental
factor evaluated by means of Equation (8) and the slip factor evaluated by Busemann. This
one. The case
comparison without the slip
is illustrated phenomenon
in Figure overestimates
28, where the all the
predicted curves experimental
are compared operating
with respect to the points,
whereasexperimental
the case with one.theThe
slipcase
factor by Busemann
without underestimates
the slip phenomenon them. Applying
overestimates the proposed slip
all the experimental
operating points,
factor formulation, thewhereas
model the case with
shows the slip factor behavior,
an intermediate by Busemann underestimates
resulting them.error
in a mean Applying
lower than
the proposed
other cases, slip factor formulation,
as summarized in Table 5. the model shows
Obviously, the an intermediatefitting
least-square behavior, resulting
curve in a mean(8) was
of Equation
error lower than other cases, as summarized in Table 5. Obviously, the least-square fitting curve of
obtained thanks to numerical investigations on a single machine. Hence, further investigations need
Equation (8) was obtained thanks to numerical investigations on a single machine. Hence, further
to be carried out in order to have a general applicability for machines characterized by different nq ,
investigations need to be carried out in order to have a general applicability for machines
different numbers of
characterized byblades,
differentand
𝑛 , outlet blade
different angles.
numbers of blades, and outlet blade angles.

Figure
Figure 28. 28. Comparison
Comparison of different
of different slipslip factordefinitions
factor definitions in
in the
the1-D
1-Dcurve
curveprediction model
prediction (i.e., σ(i.e.,
model turbine,σ
turbine ,
σ proposed
σ proposed by Busemann
by Busemann andand
σ =σ1)= 1) withthe
with theexperimental
experimental curve.
curve.

Table 5.Table 5. Comparison


Comparison of different
of different slipslip factordefinition
factor definition applied
appliedtoto
thethe
1-D1-D
proposed prediction
proposed model model
prediction
for PaT characteristic curve in terms of head prediction
for PaT characteristic curve in terms of head prediction errors. errors.

Q/QBEP,T 0.79 0.91 1 1.0


Q/QBEP,T 0.79 0.91 1 1.0
Proposed Model (No sigma) 4.6% 2.8% 0.8% 0.1%
Proposed Model (No sigma)
Proposed Model (𝜎 ) 4.6%2.9% 2.8% −2.9% 0.8%
−0.4% −3.7% 0.1%
Proposed Model (σ
Proposed Model ) Busemann)2.9%
turbine(σ −2.6% −0.4%−4.6%−2.9%
−3.2% −5.0% −3.7%
Proposed Model (σ Busemann) −2.6% −3.2% −4.6% −5.0%
7. Conclusions
7. Conclusions
This work shows numerical investigations of a centrifugal pump operating in reverse mode—
namely, working as a turbine—with emphasis on the flow field at the runner outlet. Based on the
This work shows numerical investigations of a centrifugal pump operating in reverse
results of these simulations, a correlation between slip factor and flow rate was introduced,
mode—namely,
improving working as of
the accuracy a aturbine—with emphasis
1-D model for the predictionon
of the
PaT flow field atcurves.
characteristic the runner outlet. Based
on the results of these simulations, a correlation between slip factor and flow rate was introduced,
improving the accuracy of a 1-D model for the prediction of PaT characteristic curves.
Firstly, possible sources of deflection of the flow at the outlet of PaT runners were discussed,
according to the mechanism evidenced by Ventrone [19] and Shi et al. [20]. Given that the flow
deflection at the runner outlet was associated with the finite number of blades, numerical investigations
were carried out in order to quantify this effect. Thus, numerical analyses of simplified single vane
geometries (no blade thickness and no rounded trailing edge influence) with different numbers of
blades, designed with a potential flow tool, were performed. 3D steady-state simulations confirmed
that lower the number of blades, higher the deflection.
After that, the flow field at the outlet of a commercial double suction centrifugal pump has
been investigated by solving the 3D U-RANS equations through the entire machine. The presence of
the slip phenomenon was pointed out by the evaluation of the tangential and axial absolute vectors
components obtained via the numerical simulations. These results have been compared with the
data resulting from the blade congruent flow theory. This analysis allowed the introduction of the
slip factor, σturbine , which has been helpful in order to improve a 1-D model developed to predict the
characteristic curve of a PaT starting from the pump geometry and its characteristic curves. Moreover,
Water 2019, 11, 565 21 of 23

a least-squares fitting curve, which describes the turbine slip factor, σturbine , at various flow rates, has
been proposed by the authors. Indeed, the proposed 1-D model, including the slip factor correlation,
shows good results over a wide range of mass flow rates, i.e., an improvement of the head prediction
at the BEP (−5% error compared to no-slip, σ = 1 = const, 1-D models) and at part-load (−8% error)
compared to the σ correlations proposed by Stodola [23] and Busemann [24].
In the future work, further investigations of a simplified geometry by varying its number of blade
and its flow rate will be carried out in order to draw a map of the PaT slip factor, which could be useful
for manufacturers.

Author Contributions: T.C. and M.S. wrote the paper; T.C. performed CFD analysis and M.S. developed the 1-D
prediction model; G.P., M.T. and S.R. co-operated with T.C. and M.S. on the theoretical aspects concerning this
topic. M.T., S.M.C. and B.F. supervised the work.
Funding: This research received no external funding.
Acknowledgments: The authors gratefully acknowledge Nuovo Pignone for their experimental support. During
the review process of this article, we were saddened by the death of our colleague and friend, Bernardo Fortunato,
who was one of the co-authors of this paper. Bernardo passed away on 31 January 2019, at the age of 68. Full
professor at the age of 36, he was a preeminent researcher at the Polytechnic of Bari, serving it in several academic
roles and as a responsible of research projects, always promoting studies on renewable energies. He supported
and encouraged us in carrying out activities in the field of hydraulic pumps and PaTs, which underlie this paper.
We had the pleasure and the honor to work with him and we hope to continue his work despite of his absence.
Conflicts of Interest: The authors declare no conflict of interest.

Abbreviations
BEP Best Efficiency Point
CFD Computational Fluid Dynamic
MHP Micro Hydropower Plant
PaT Pump as Turbine
WDN Water Distribution Network

Nomenclature
ac [m/s2 ] Coriolis acceleration
aw [m/s2 ] Relative acceleration
C [Nm] Torque
c [m/s] Absolute velocity
D [m] Diameter
fq [-] Number of impeller entries
G [kg/s] Mass flow rate
g [m/s2 ] Gravitational acceleration
H [m] Head
l [m] Blade height
n [rpm] Rotational speed
Specific speed number (n in rpm,
nq [-]
Q in m3 /s and H in m)
p [Pa] Pressure
Q [m3 /s] Flow rate
r [m] Radius
u [m/s] Tangential velocity
w [m/s] Relative velocity
Y [J/kg] Specific work
Zblades [-] Number of blades
Zrunner [m] Volute loss
Zvolute [m] Runner loss
Water 2019, 11, 565 22 of 23

Greek Symbols
α [rad] Absolute flow angle
β [rad] Relative flow angle
η [-] Efficiency
θ [rad] Angular coordinate
ρ [kg/m3 ] Density
σ [-] Slip factor
ω [rad/s] Angular velocity

Subscripts and Superscripts


1 Outlet section of the runner
2 Inlet section of the runner
B Blade
m Velocity meridional component
P Pump
T Turbine
th Theoretical
u Velocity tangential component
tot Total

References
1. Laux, C.H. Reverse-running pumps as energy recovery turbines. Sulzer Tech. Rep. 1982, 2, 23–27.
2. Apfelbacher, R.; Etzold, F. Energy-Saving, Shock-Free Throttling with the Aid of a Reverse Running
Centrifugal Pump. KSB Tech. Rep. 1988, 24e, 33–41.
3. Binama, M.; Su, W.T.; Li, X.B.; Li, F.C.; Wei, X.Z.; An, Z. Investigation on pump as turbine (PAT) technical
aspects for micro hydropower schemes: A state-of-the-art review. Renew. Sustain. Energy Rev. 2017, 79,
148–179. [CrossRef]
4. Barbarelli, S.; Amelio, M.; Florio, G. Experimental activity at test rig validating correlations to select pumps
running as turbines in microhydro plants. Energy Conversion and Management 2017, 107, 781–797. [CrossRef]
5. Derakhshan, S.; Nourbakhsh, A. Experimental study of characteristic curves of centrifugal pumps working
as turbines in different specific speeds. Exp. Therm. Fluid Sci. 2008, 32, 800–807. [CrossRef]
6. Nautiyal, H.; Kumar, A.; Yadav, S. Experimental investigation of centrifugal pump working as turbine for
small hydropower systems. Energy Sci. Technol. 2011, 1, 79–86.
7. Pugliese, F.; De Paola, F.; Fontana, N.; Giugni, M.; Marini, G. Experimental characterization of two Pumps as
Turbines for hydropower generation. Renew. Energy 2016, 99, 180–187. [CrossRef]
8. Rossi, M.; Righetti, M.; Renzi, M. Pump-as-Turbine for energy recovery applications: The case study of an
aqueduct. Energy Procedia 2016, 101, 1207–1214. [CrossRef]
9. Singh, P.; Nestmann, F. An optimization routine on prediction and selection model for the turbine operation
of centrifugal pumps. Exp. Therm. Fluid Sci. 2010, 34, 152–164. [CrossRef]
10. Shi, G.; Liu, X.; Wang, Z.; Liu, Y. Conversion relation of centrifugal pumps as hydraulic turbines based on
the amplification coefficient. Adv. Mech. Eng. 2017, 9, 1–8. [CrossRef]
11. Yang, S.S.; Derakhshan, S.; Kong, F.Y. Theoretical, numerical and experimental prediction of pump as turbine
performance. Renew. Energy 2012, 48, 507–513. [CrossRef]
12. Tan, X.; Engeda, A. Performance of centrifugal pumps running in reverse as turbine: Part II—Systematic
specific speed and specific diameter based performance prediction. Renew. Energy 2016, 99, 188–197.
[CrossRef]
13. Carravetta, A.; Fecarotta, O.; Ramos, H.M. A new low-cost installation scheme of PATs for pico-hydropower
to recover energy in residential areas. Renew. Energy 2018, 125, 1003–1014. [CrossRef]
14. Balacco, G. Performance Prediction of a Pump as Turbine Sensitivity Analysis Based on Artificial Neural
Networks and Evolutionary Polynomial Regression. Energies 2018, 11, 3497. [CrossRef]
Water 2019, 11, 565 23 of 23

15. Rossi, M.; Renzi, M. A general methodology for performance prediction of Pumps-as-Turbines using Artificial
Neural Networks. Renew. Energy 2018, 128, 265–274. [CrossRef]
16. Stefanizzi, M.; Torresi, M.; Fortunato, B.; Camporeale, S.M. Experimental investigation and performance
prediction modeling of a single stage centrifugal pump operating as turbine. Energy Procedia 2017, 126,
589–596. [CrossRef]
17. Barbarelli, S.; Amelio, M.; Florio, G. Predictive model estimating the performance of centrifugal pumps used
as turbines. Energy 2016, 107, 103–121. [CrossRef]
18. Gülich, J.F. Centrifugal Pumps, 2nd ed.; Springer: Berlin, Germany, 2008.
19. Ventrone, G. Deviazione della corrente relative nelle giranti delle turbine Francis. L’Energia Elettr. 1972, 9,
569–574.
20. Shi, G.; Liu, X.; Yang, J.; Miao, S.; Li, J. Theoretical research of hydraulic turbine performance based on slip
factor within centripetal impeller. Adv. Mech. Eng. 2015, 7. [CrossRef]
21. Capurso, T.; Stefanizzi, M.; Torresi, M.; Pascazio, G.; Caramia, G.; Camporeale, S.M.; Fortunato, B.;
Bergamini, L. How to improve the Performance Prediction of a Pump as Turbine by Considering the
Slip Phenomenon. Proceedings 2018, 2, 683. [CrossRef]
22. Stefanizzi, M.; Capurso, T.; Torresi, M.; Pascazio, G.; Ranaldo, S.; Camporeale, S.M.; Fortunato, B.;
Monteriso, R. Development of a 1-D prediction model of a PaT performance under design and off-design
conditions. Proceedings 2018, 2, 682. [CrossRef]
23. Busemann, A. Das Forderhohenverhaltnis Radialer Kreiselpumpen mit Logarithmisch-Spiraligen Schaufeln.
ZAMM J. Appl. Math. Mech. 1928, 8, 372–384. [CrossRef]
24. Stodola, A. Dampf und Gasturbinen; Springer: Berlin, Germany, 1924.
25. D’Agostino, L.; Pasini, A.; Valentini, D. A reduced order model for preliminary design and performance
prediction of radial turbopumps. In Proceedings of the 47th AIAA/ASME/SAE/ASEE Joint Propulsion
Conference & Exhibit, San Diego, CA, USA, 31 July–3 August 2011. [CrossRef]
26. Menter, F.R.; Kuntz, M.; Langtry, R. Ten Years of Industrial Experience with the SST Turbulence Model.
Turbul. Heat Mass Transf. 2003, 4, 625–632.
27. Huang, P.; Bardina, J.; Coakley, T. Turbulence Modeling Validation, Testing, and Development; NASA Technical
Memorandum; NASA Ames Research Center: Moffett Field, CA, USA, 1997; Volume 110446.
28. Capurso, T.; Torresi, M.; Bergamini, L. Design and CFD performance analysis of a novel impeller for double
suction centrifugal pumps. Nuclear Eng. Design 2018, 341, 155–166. [CrossRef]
29. Adams, T.; Grant, C.; Watson, H. A Simple Algorithm to Relate Measured Surface Roughness to Equivalent
Sand-grain Roughness. Int. J. Mech. Eng. Mechatron. 2012, 1, 66–71. [CrossRef]
30. Nilsson, H. Evaluation of OpenFOAM for CFD of turbulent flow in water turbines. In Proceedings of the
23rd IAHR Symposium, Yokohama, Japan, 17–21 October 2006.

© 2019 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/).

You might also like