28th IAHR Symposium on Hydraulic Machinery and Systems

Simulation of energy recovery on water utility

networks by a micro-turbine with counter-rotating
runners.
L Andolfatto1 , E Vagnoni1 , V Hasmatuchi2 , C Münch-Alligné2 and F
Avellan1
1
École polytechnique fédérale de Lausanne, Laboratory for Hydraulic Machines, Avenue de
Cour 33 bis, 1007 Lausanne, Switzerland
2
University of Applied Sciences Valais, Route du Rawyl 47, 1950 Sion, Switzerland
E-mail: loic.andolfatto@epfl.ch

Abstract. Wherever relief valves and other energy dissipation devices are installed to limit
the pressure, water utility networks provide unexploited hydropower potentials. This is mainly
due to a lack of economically viable technologies for energy recovery in the pico and micro
hydropower range below 100 kW. Micro-turbine with counter-rotating runners proved suitable
to harvest these potentials with limited investments and almost no environmental impact. An
appropriate command strategy must therefore be applied to maximize the recovered energy.
This paper deals with the construction of a Virtual Energy Recovery Station (VERS) model
to simulate the energy recovery on a given installation site. It includes models of the turbine,
of the water consumption and it allows to implement various command strategies. The VERS
can serve various purposes. The fine tuning of the command algorithm for a specific installation
site is demonstrated in the paper.

1. Introduction
In hilly or mountainous area, the difference of altitude between spring catching areas or fresh
water reservoirs and water consumption areas sometimes impose to regulate the pressure in the
consumption area. A relief valve is installed to keep the pressure at preg in the consumption
area, as pictured in Figure 1.
Then, with a discharge Q in the pipe, the hydraulic power Ph dissipated by the valve is given
by:
Q2 preg
 
Ph = ρ Q gH − keq − (1)
2 ρ
with H the gross head, keq the equivalent loss coefficient of the head water infrastructure, ρ the
density of the fluid and g the gravity. Instead of dissipating this power, it can be recovered to
produce electricity with an appropriate power plant, as pictured in Figure 1. The environmental
impact related to the energy produced by this mean can reach negligible levels as it uses already
existing infrastructures.
Many successful application of energy recovery on water utility networks have been
reported [1]. But the limited revenues generated by the power station allow to reach economical

1701
28th IAHR Symposium on Hydraulic Machinery and Systems

Figure 1. Water utility network with hydro power potential.

viability only when the available power is higher than a given threshold related to the capital
expenditure for the installation and its operating cost.
The concept of micro-turbine with counter-rotating runners pictured in Figure 2 and
presented in depth in [2, 3] is a candidate for energy recovery on drinking water network, even
for an available power in the range of 5 kW to 25 kW. It features a compact axial architecture
ensuring a lean in-line installation on existing facilities, therefore limiting the capital expenditure
and the environmental impact of the infrastructure. Such a machine operated at variable speed
can cover a wide operating domain to fit with the need of installations with uncontrolled
consumer-driven flow discharge. Using a variable speed ratio between the two runners also
enhances the overall efficiency of the machine [4], increasing the expected revenues.

Figure 2. Axial micro-turbine

with counter-rotating runners.

In order to maintain the machine cost and the maintenance cost as low as possible, the micro-
turbine is commanded with a Maximum Power Point Tracking (MPPT) algorithm. It allows to
avoid monitoring the operating conditions of discharge and available head and thus saving on
the associated sensors. The principle of the MPPT is described in section 2.
To find the optimum values of the MPPT parameters, the entire Energy Recovery Station
system is modeled and simulated. A stochastic process modeling the consumer-driven discharge
of an instrumented pilot site is built. Discharge trajectories can be generated to estimate

1702
 28th IAHR Symposium on Hydraulic Machinery and Systems

NA

DN t P
B+
P
NB
NB
0
t A- P
P

A+
P P NA
Dt0 Dt DN
B+ B-
P t P
0 B-
P A+ P
P A-
P

Figure 3. Chronogram of a perturb and Figure 4. Estimation of the steepest direction

observe cycle of the MPPT. according to the observed output power.

the energy recovered with the micro-turbine with various MPPT settings. The simulation is
presented in section 3 and the results obtained are detailed in section 4.

2. Maximum Power Point Tracking command

2.1. Principle
The principle of Maximum Power Point Tracking has been initially proposed by Cherdak and
Douglas [5] to control the output voltage of power sources in order to maximize the output
power. It consists in applying a perturbation to the controllable parameters of the system while
observing the associated alteration of the output power, what is sometimes refereed to as perturb
and observe. The controller moves in the steepest gradient direction to improve the output
power. It has been widely implemented for the control of photovoltaic system [6]. This solution
has been derived for small hydro applications, such as the one proposed by Borghetti et al. [7]
for instance.
As far as the micro-turbine is concerned, the controllable parameters are the two runners
rotational speed, denoted NA and NB respectively. The output power P is the sum of the power
PA and PB delivered by the generators coupled to the two runners. The rotational speeds are
perturbed by ∆N in the positive and negative directions and the associated output power are
recorded, according to the chronogram of Figure 3. The steepest direction ∇P is determined
using eq. (2).
∂P
 

∇P =  ∂N A (2)

∂P 
∂NB
Then the new set point for the time step i is computed using the direction di as defined in
eq. (3):
 
NA
Ni = = N i−1 + α · di = N i−1 + α · (∇P i−1 + β · di−1 ) (3)
NB i

1703
28th IAHR Symposium on Hydraulic Machinery and Systems

Table 1. Parameters of the MPPT algorithm.

Parameter Unit Description
∆t [s] Duration of the period at current operating set point N i
∆t0 [s] Duration of the period at each perturbed operating set point
∆N [min−1 ] Magnitude of the rotational speed perturbation
α [W−1 min−2 ] Magnitude of the set point step in the direction di
β [−] Conjugate gradient parameter for the direction di computation

To keep the implementation simple and limit the resources required for an embedded
controller, all the parameters summarized in Table 1 are considered to be fixed once for all. The
effect of each parameter is investigated in the following sections. Apart from these parameters,
the other degree of freedom in the MPPT design investigated in this paper is the method used
to compute the steepest direction ∇P . First and second orders approximations are explained in
subsections 2.2 and 2.3 respectively.

2.2. First order MPPT

Considering the output power as locally linear is a simple and straightforward approach to
estimate its gradient. The central difference scheme of eq. (4), with the notations from Figure 3,
directly provides the steepest direction.
 A+
P − A−P

1
∇P = (4)
2∆N B+P − B−P

2.3. Second order MPPT

As the gradient is likely to be estimated in neighborhood of local optima, as in Figure 4 for
the runner A, the first order approach may lead to incorrect estimations. A second order
approximation as defined in eq. (5) is therefore also investigated.
   
T a T c e/2
P (N + ∆N ) = P0 + ∆N · + ∆N · ·∆N (5)
b e/2 d

The coefficient e cannot be identified as the perturbation ∆N only involves one runner at a
time. The others coefficients are identified by solving the least square problem of eq. (6).
X 2
k
ε= P − P (N +k∆N ) , k ∈ {A+, A−, B+, B−} (6)

The gradient is finally obtained according to the expression of eq. (7).

 
a
∇P = (7)
b

This generic formulation requires more computing efforts than the first order but it is also more
robust in the neighborhood of the operating domain boundaries where all the k∆N are not
accessible and some alternative perturbations must be used.

1704
 28th IAHR Symposium on Hydraulic Machinery and Systems

Command

 Q2   N A , N B   PA , PB   preg 
 Q, gH  keq 
(Q, gH )  2   Q, 
  
Head water Head water
Micro-turbine Relief valve Consumers
reservoir duct pipe
Energy Recovery Station

Figure 5. Energy Recovery Station within its environment.

3. Simulation of the Energy Recovery Station

3.1. Description of the system
The primary function of the water utility network consists in delivering water to consumers.
This function is characterized by a target pressure preg at the tail water side of the relief valve
and a discharge Q that is consumer-driven. The water is conveyed from a head water reservoir to
the relief valve chamber through a duct pipe characterized by its equivalent losses coefficient keq .
Then, the command for micro-turbine runners rotational speeds NA and NB must be computed
in order to maximize the sum of the output power PA and PB from the two generators. To keep
the pressure higher or equal to preg at the inlet of the relief valve, the maximum specific energy
that can be extracted by the micro-turbine Eav is defined by:

Q2 preg
−
Eav = gH − keq (8)
2 ρ
This is summarized in Figure 5. The net head H and the equivalent losses coefficient keq
can be analytically computed according to the characteristics of the installation. Alternatively,
they can be identified from on site measurements as in [4]. In order to simulate the system,
the proposed analytic model of the micro-turbine and stochastic model of the consumer-driven
discharge are presented in the following subsections.

3.2. Analytic model of the micro-turbine

The simulation of the system is based on a model of the micro-turbine that predicts the output
power P – sum of PA and PB – as a function of the discharge Q and of the velocity set point N .
Additionally,the extracted specific energy is also estimated in order to check that the constraint
on the tail water pressure is satisfied.
An experimental campaign – explained in depth in [4] – is conducted on a laboratory test rig in
order to build such a model of a micro-turbine prototype. The entire operating range from Nmin
to Nmax on both runners is explored at n+1 levels of extracted specific energy E = {E0 , . . . , En }.
Analytic models based on Hermite polynomials interpolation are build accordingly:
• for the discharge passing through the micro-turbine : {QEi (N )}Ei ∈E .
• for the output power of the first runner : {PA,Ei (N )}Ei ∈E ;
• for the output power of the second runner : {PB,Ei (N )}Ei ∈E ;
The discharge models have their values sorted according to the specific energy level, as defined
in eq. (9).
∀N ∈ [Nmin , Nmax ]2 , Ei > Ej ⇐⇒ QEi (N ) > QEj (N ) (9)
This property is exploited in the step (a) of the algorithm detailed in Figure 6. The estimated
level of extracted specific energy Eest is estimated thanks to a linear interpolation in step (b)
and (c). The output power PA and PB are then predicted accordingly in step (d).

1705
28th IAHR Symposium on Hydraulic Machinery and Systems

Read Q, N
if Q  QE1 ( N ), i  1
Find i (a) 
else i QEi1 ( N )  Q  QEi ( N )
Q  QEi ( N )
Compute a (b) a 
QEi1 ( N )  QEi ( N )

Compute Eest (c) Eest  a  Ei 1  (1  a )  Ei

 PA  a  PA, Ei1 ( N )  (1  a )  PA, Ei ( N )

Compute PA, PB (d)  Figure 6. Prediction of the output power
 PB  a  PB , Ei1 ( N )  (1  a )  PB , Ei ( N )
PA and PB according to the discharge Q
Return PA, PB, Eest and of the velocity set point N .

A similar approach is also used to correct the ascent direction di of the MPPT in order to
satisfy the constraint Eest < Eav . Nonetheless, this is not detailed in this paper.

3.3. Stochastic model of the discharge

To complete the proposed framework for the simulation of the energy recovery, the consumer-
driven discharge is also modeled so that realistic time series trajectories {Q(k·δt)}k∈{1,...,T } can
be generated.
The model of time varying discharge must cope with fluctuations in the consumption at
various time scale: seasonal, intra-week, intra-day and instantaneous. As the simulations aims
at estimating integral quantities such as the energy production over a given period for instance,
the average of the simulated discharge time series must also be unbiased.
Developing a multi-agent model considering each individual consumer or batches of consumers
as proposed by Davis [8] is theoretically possible but requires a huge amount of information about
the network. The global approach considering the discharge as a stochastic process is preferred,
as it only requires some measured time series to be tuned. Even if auto-regressive moving-
average models (ARMA) are not reported as the most efficient for monthly prediction of the
discharge [9], they are selected for their versatility when it comes to short δt period of around
one second as in this work.
ARMA is a popular model for stationary hydrologic time series [9]. The eq. (10) provides
the formulation defined by Whittle [10].
p
X q
X
Q(tj ) = εj + ϕi · Q(tj−i ) + θi · εj−i (10)
i=1 i=1
The model is characterized by its auto-regressive order p, its moving-average order q, its
coefficients ϕi and θi . The random variable εk is a white noise. The identification of the model
is made by maximizing the corrected Akaike information criterion AICc [11]:

2k(k + 1
AICc = 2k − 2 ln L + (11)
n−k−1
with k = p + q + 1, n is the number of points in the measured time series and L is the likelihood
of the modeled time series with respect to the measured one.
To capture the discharge fluctuation both at the daily scale and at the second scale, two
sets of measurements have been performed on an instrumented pilot site. The first data set
Q15 = {Q15,k }k∈{1,...,350 040} covers one year from June 2014 to June 2015 and provides the
average discharge on the 35’040 periods of 15 minutes. A sample along 14 days is pictured

1706
 28th IAHR Symposium on Hydraulic Machinery and Systems

103 × Q [m3 ·s−1 ]

14

12

10

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
Days

Figure 7. Sample of the 15 minutes averaged measured discharge during 14 days.

δ Q2
δ Q1

δ t1 δ t2

Figure 8. Sample of the instantaneous discharge measured over 15 minutes.

in Figure 7. The second data set Q1 = {Q1,m }m∈{1,...,860 400} covers one day in June 2015 and
provides the instantaneous discharge with a 1 Hz acquisition frequency. A sample along 15
minutes is pictured in Figure 8. Two separate models have been build and merged to yield the
stochastic model of the consumer-driven discharge.
The average value of any contiguous subset of the measured average discharge Q15 is not
constant as some periodical fluctuations are experienced. It is therefore not stationary and can
not be modeled directly with an ARMA model. In order to overcome this hurdle, the eq. (12)
provides a multiscale model accounting for the periodical fluctuations.

Q15 (t) = Q + Qs (t) + Qd (t) + Qh (t) + ws (t) · wd (t) · wh (t) · q15 (t) (12)
The discharge trend functions Qs , Qd , and Qh are periodic with a mean value equal to zero over
their period and they are constant over each of their sub-periods. The weight functions ws , wd ,
and wh are periodic with a mean value equal to one over their period and they are constant
over each of their sub-periods. Characteristics of these functions are detailed in the Table 2.
They are iteratively identified – from the seasonal trend to the hourly trend – to match mean
values and standard deviations of the measured discharge Q15 . The residual q15 finally follows
a stationary time series. It is modeled with an ARMA process.

1707
28th IAHR Symposium on Hydraulic Machinery and Systems

Table 2. Characteristics of the periodical discharge trend and weight functions.

Trend Functions Period Sub-period
Seasonal Qs , ws 1 year 1 week
Daily Qd , wd 1 week 1 day
Hourly Qh , wh 1 day 1 hour

The second data set Q1 is used to identify the values of the time series within each period of
15 minutes previously defined. A piecewise linear trend Qpl,k (t) is identified on each subset of
Q1 covering 15 minutes. The trend and its four parameters are pictured in Figure 8. The mean
value of Qpl,k equals the mean discharge Q15,k obtained thanks to eq. (12) and the continuity
with the previous period is ensured by imposing Qpl,k (0). A composed distribution C for the
parameters δt1,k , δt2,k , δQ1,k and δQ2,k is identified according to Q1 . Then, the stochastic model
of the discharge is finally given by eq. (13).

Q(t) = Qpl,k (t) + ·Qf (t) (13)

The residual instantaneous fluctuations Qf can be considered as a stationary time series. It
is therefore possible to model it with an ARMA process.
Finally, the generation of a time series for simulation purpose is sequential. First, a future
trajectory generated from the ARMA model of q15 is transformed thanks to eq. (12) providing
mean values over 15 minutes for generated time series Qg . Then, for each time step, a set of
trend parameters δt1,k , δt2,k , δQ1,k and δQ2,k is drawn from the composed distribution C to form
the trend and the final time series is obtained thanks to the generation of a future trajectory
from the ARMA model of Qf , according to eq. (13).

3.4. Implementation
All the elements described in the previous subsections have been implemented in Python [12].
The statistical analysis and the identification of the ARMA models relies on the OpenTURNS
library [13]. The time step for the simulation has been set to one second. Simulation results are
shown in Figure 9, where the rotational speeds are following the discharge fluctuation to track
the best output power setting.

4. Results
For both first and second order gradient method, the optimum MPPT parameters have been
searched considering.
• β = 0 to evaluate the relevance of a direct use of the steepest direction;
• β = 0.15 to estimate the impact of introducing the conjugate gradient.
For the four cases addressed, ∆t0 , ∆t, ∆N and α have been optimized in order to maximize the
mean daily energy production W . It has been estimated according to simulation running on one
year, with the same trajectory of consumer-driven discharge.
The results are presented in Table 3.
Even though the conjugate gradient and the second order provide the best performance, the
differences are almost negligible. In terms of MPPT parameters selection, the four cases also
led to similar results.

1708
 28th IAHR Symposium on Hydraulic Machinery and Systems

11
103 × Q [m3 ·s−1 ]

Q
10
9
8
7
6
0 500 1000 1500 2000 2500 3000 3500
t [s]
3500
3000 NA NB
N [min−1 ]

2500
2000
1500
1000
500
0 500 1000 1500 2000 2500 3000 3500
t [s]
2000
PA PB P
1500
P [W]

1000

500

0
0 500 1000 1500 2000 2500 3000 3500
t [s]

Figure 9. Simulation of the first 3600s after a start at N = (2000, 2000).

Table 3. Optimal parameters of the MPPT and associated mean daily energy production W .
Order β ∆t0 ∆t ∆N α W
[−] [−] [s] [s] [min−1 ] [W−1 min−2 ] [kWh]
1 0 8.935 7.833 104.690 259.631 19.610
2 0 8.934 7.833 104.690 259.631 19.617
1 0.15 8.937 7.832 104.690 259.631 19.611
2 0.15 8.936 7.833 104.690 259.631 19.618

5. Conclusion
The design of a system dedicated to energy recovery on existing infrastructure requires to
consider typical constraints that are not experienced in the case of conventional hydropower
infrastructures. This paper addresses the energy recovery in drinking water networks, where the
primary function of the infrastructure is to deliver water to consumers. The primary function
must be preserved at all time and the Energy Recovery Station must be designed consequently.
A complete simulation framework is proposed in order to validate and to tune the command

1709
28th IAHR Symposium on Hydraulic Machinery and Systems

strategy of the Energy Recovery Station. It is fueled by a discharge time series generator
to reproduce the behavior of consumers with specific statistical characteristics. A MPPT
command algorithm to control a micro-turbine with counter-rotating runners is presented. The
performance of the system is estimated and tuned according to the characteristics of an identified
pilot site.
The successful use of the proposed approach for the dedicated tuning of a station for an
identified site should not hide the actual potential of the proposed approach. A similar framework
can also be applied to simulate the operation of several station on a water utility network.
Other components, such as intermediate reservoirs or pumps, can also be modeled on the same
principal, allowing to simulate a cluster of small hydropower devices connected to a smart energy
network.

Acknowledgment
The research leading to the results published in this paper has received funding from the Swiss
Competence Center for Energy Research Supply of Electricity (SCCER SoE) granted by the
Swiss Commission for Technology and Innovation (CTI); the Ark, the foundation for innovation
of Valais Canton and from the Swiss Commission for Technology and Innovation as part of
the DuoTurbo project number 17197.1 PFEN IW. The authors would also like to thank the
implementation partners Jacquier-Luisier SA, Valelectric SA and Telsa SA for their commitment
and their support in the DuoTurbo project.

References
[1] Hintermann M 1994 L’eau potable génératrice d’électricité: Inventaire et étude du potentiel des
usines électriques sur l’alimentation en eau potable en suisse SFOE DIANE 10 publication URL
http://www.bfe.admin.ch/kleinwasserkraft/03834/04171/index.html?lang=en
[2] Münch-Alligné C, Richard S, Meier B, Hasmatuchi V and Avellan F 2014 Numerical simulations of a counter-
rotating micro-turbine Advances in Hydroinformatics Springer Hydrogeology ed Gourbesville P, Cunge J
and Caignaert G (Springer Singapore) pp 363–373 ISBN 978-981-4451-41-3
[3] Andolfatto L, Euzenat C, Vagnoni E, Münch-Aligné C and Avellan F 2015 A mixed standard/custom
design strategy to minimize cost and maximize efficiency for picohydro power potential harvesting 5th
International Youth Conference on Energy (IYCE), 2015 pp 1–8
[4] Andolfatto L, Delgado J, Vagnoni E, Münch-Aligné C and Avellan F 2015 Analytical hill chart towards the
maximisation of energy recovery on water utility networks with counter rotating micro-turbine 36th IAHR
World Congress, The Hague, The Netherlands iSBN: 978-90-824846-0-1
[5] Cherdak A S and Douglas J L 1971 Maximum power point tracker uS Patent 3,566,143 URL
https://www.google.com/patents/US3566143
[6] Salas V, Olas E, Barrado A and Lzaro A 2006 Solar Energy Materials and Solar Cells 90 1555 – 1578 ISSN
0927-0248
[7] Borghetti A, Di Silvestro M, Naldi G, Paolone M and Alberti M 2008 Maximum Efficiency Point Tracking
for Adjustable-Speed Small Hydro Power Plant 16th Power Systems Computation Conference (PSCC’08)
[8] Davis D 2000 Computers, Environment and Urban Systems 24 173 – 190 ISSN 0198-9715
[9] Wang W C, Chau K W, Cheng C T and Qiu L 2009 Journal of Hydrology 374 294 – 306 ISSN 0022-1694
[10] Whittle P 1951 Hypothesis testing in time series analysis vol 4 (Almqvist & Wiksells)
[11] Cavanaugh J E 1997 Statistics & Probability Letters 33 201 – 208 ISSN 0167-7152
[12] Python Software Foundation Python language reference, version 2.7 Available at http://www.python.org
[13] OpenTURNS: an Open source initiative for the Treatment of Uncertainty, Risks’N Statistics
http://doc.openturns.org/openturns-latest/sphinx/contents.html, Consulted in January 2016.

1710