Energy production assessment in a complex hydropower development

Due to massive and/or unexpected energy production from renewable sources, especially the fluctuant wind and solar power generation, conventional hydropower developments and pumped storage power developments became crucial in the stabilization of electric grids (Anagnostopoulos & Papantonis 2012; Boamba et al. 2014; Constantin et al. 2015; Pérez-Díaz et al. 2015). To fulfil this special purpose, hydropower units are often running at partial loads, instead of running at their best efficiency point, to cover a wide range of operating parameters variation. So, decision makers are interested to determine the appropriate variation range of those parameters to ensure safe operation of the system. Also, hydropower producers are interested to establish a strategy for bidding on the day-ahead electricity market, based on a decision support system that relies on the water availability (Bozorg Haddad et al. 2008; Jamshid Mousavi & Shourian 2010; Tica et al. 2017). Their goal is to settle and/or to optimize the amount of power that can be sold/bought in the spot market, as well as the amount of power that can be kept for ancillary services (Chazarra et al. 2016).

The present study focuses on a complex hydropower development (HPD), namely a pumped storage hydropower open system. In the attempt to assess the operation mode of that system, for different scenarios linked to the variation of water level in the reservoirs, the hydraulic analysis was performed based on an equivalent numerical model. The results, analysed in terms of the balance between the energy production and the energy consumption, can be implemented within a decision support system (Popa et al. 2008) to enhance the overall efficient operation of the hydropower system.

The numerical model attached to the studied HPD is described in the next section, which includes the HPD real scheme and the simplified model, the HPD operation scenarios, and the modelling approach. The results are presented and discussed in the third section, while the concluding remarks are drawn at the end of the paper.

Gâlceag HPD consists of three reservoirs, one high head hydropower plant (HPP), and one pumping station (PS). The resulting multi-reservoir system represents the upstream part of a hydropower cascade built on the Sebeş River in Romania (Popa et al. 2017). The above hydropower cascade contains three successive high head HPPs and ends by a low head (toe dam type) HPP (Hidroconstructia 2005). Gâlceag HPP is a peak load HPP, operating 2–3 h a day. The present study focuses exclusively on Gâlceag HPD, for which the hydraulic analysis is conducted during 1 h, assuming that the water exiting the HPP is retained in the nearby downstream reservoir. This assumption is equivalent with keeping shutdown, over 1 h, the other HPPs located downstream on the hydropower cascade.

The scheme of Gâlceag HPD is displayed in Figure 1 (left frame).

The main components of Gâlceag HPD are briefly described below, based on available technical data (Hidroconstructia 2005; Dunca et al. 2010, 2018; Bucur et al. 2014):

• Gâlceag HPP (year of commissioning: 1980) − a high head underground HPP of 150 MW installed capacity and 260 GWh average annual energy production; it is powered by two Francis turbines, of 75 MW each, at 465 m head and 22.8 m³/s flow rate.

• Gâlceag PS (year of commissioning: 2003) − a PS with two double-entry two-stage centrifugal pumps, each of 10 MW at 260 m pumping head and 3 m³/s discharge; PS elevation is 980 m above sea level (a.s.l.).

• Oaşa Reservoir − the biggest upstream reservoir of the entire hydropower cascade on the Sebeş River, with an active storage capacity of 136.2 MCM at 1,255 m a.s.l. normal retention level; water is retained by Oaşa Dam (a rockfill dam 91 m height, with a reinforced concrete mask; year of commissioning: 1972).

• Cugir Reservoir − a storage reservoir, used as PS suction reservoir (0.85 MCM active storage capacity, at the normal retention level of 1,007 m a.s.l.); water is retained by Cugir Dam (an arch dam 48 m height).

• Tău Reservoir − the downstream reservoir within Gâlceag HPD scheme, with an active storage capacity of 21.3 MCM, at 790 m a.s.l. normal retention level; water is retained by Tău Dam (an arch dam 78 m height).

• Gâlceag HPP hydraulic circuit − a headrace tunnel (3.7 m diameter and a total length of 8,456 m) connects Oaşa Reservoir to the penstock (2.8 m diameter and a length of 742 m, positioned at 32° from the vertical); at the downstream part of the headrace, there is a double-chamber surge tank; the valve house at the penstock inlet contains a butterfly valve with a diameter of 2.8 m; at the penstock outlet, a distributor conveys the water inside the powerhouse towards the turbines (before each turbine, a heavy-duty spherical valve is used for shutoff purposes); a tailrace tunnel (4 m diameter and a length of 610 m) conveys the water to Tău Reservoir.

• Gâlceag PS hydraulic circuit − a headrace (2.8 m diameter and a length of 6,500 m) conveys water from the Cugir Reservoir to Gâlceag PS; at the downstream part of the headrace, before the PS powerhouse, there is a surge tank; the PS discharge pipe (1.3 m diameter and a length of 320 m) is connected to the upper part of the HPP's penstock at an elevation of 980 m a.s.l. (this unusual junction allows injection of the pumped water directly into the penstock); throttling valves (inside the PS) control the discharge.

Gâlceag HPD model is presented in Figure 1 (right frame). The hydraulic analysis conducted in this paper does not involve transients, thus no surge tank is included in the equivalent scheme. The main components of Gâlceag HPD equivalent scheme built-in EPANET are listed in the following (Dunca et al. 2018):

• Gâlceag HPP with Francis turbines T1 and T2 − the turbines are modelled in EPANET by general purpose valves (GPVs) that introduce a given head loss-flow rate curve (Rossman 2000); to mimic a turbine, each GPV operates based on the turbine head-flow rate curve; the characteristic curves of the Francis turbines from Gâlceag HPP, namely the head-flow rate curve H_t = H_t(Q_t) and the efficiency-flow rate curve η_t = η_t(Q_t), are displayed in Figure 2 (left frame); the head values H_t (in m), flow rate values Q_t (in m³/s), and efficiency values η_t (in %) were measured in situ (Bucur et al. 2012).

• Gâlceag PS with pumps P1 and P2, operating upon the pumping head-flow rate curve H_p = H_p(Q_p) and efficiency-flow rate curve η_p = η_p(Q_p), is displayed in Figure 2 (right frame); the pumping head values H_p (in m), discharge values Q_p (in m³/s), and efficiency values η_p (in %) were measured in situ (Dunca et al. 2007).

• Three reservoirs (variable head water storages) denoted from R1 to R3, modelled in EPANET by storage tanks − an upper tank R1 (Oaşa Reservoir), a storage tank R2 (Cugir Reservoir), and a lower tank R3 (Tău Reservoir); for each tank, the water level z (in m a.s.l.) varies upon the storage capacity V (in MCM), according to the capacity curve of each reservoir, as in Figure 3.

• Hydraulic circuit − pressurized pipes labelled from j = 1 to j = 5, namely HPP headrace tunnel (j = 1), PS headrace (j = 2), PS discharge pipe (j = 3), HPP penstock (j = 4), and HPP tailrace tunnel (j = 5); main junction N; a butterfly valve inside the valve house, modelled in EPANET by a throttling control valve (TCV), kept fully open.

Gâlceag HPD can operate according to the following scenarios, denoted from ① to ③, as shown in Figure 4.

① Normal operation scenario, where the HPP is fully operational, while the PS is shutdown − in this case, the HPP is working as a conventional HPP, thus the water flows only from Oaşa Reservoir (tank R1) to Tău Reservoir (tank R3), passing through the turbines T1 and T2.

② Simultaneous operation scenario, where both the HPP and the PS are fully operational − in this case, the water pumped from Cugir Reservoir (tank R2) is injected into the penstock, being added to the water flowing from Oaşa Reservoir (R1); by increasing the flow rate through the turbines, the energy output of the HPP can be increased.

③ Pumped storage scenario, where the HPP is shutdown, while the PS is fully operational − in this case, the water is pumped from Cugir Reservoir (R2) to Oaşa Reservoir (R1) to increase the storage of the upstream reservoir.

Both normal operation and simultaneous operation scenarios correspond to peak load. The pumped storage scenario corresponds to off-peak when the electrical energy production exceeds the consumption within the National Power Grid.

Primarily, for the proposed HPD, the numerical model was built in EPANET and run as an extended period simulation (EPS) (Van Zyl et al. 2006; Giustolisi et al. 2008; Todini 2011), starting from different initial levels in the upstream reservoir. The extended period simulations were performed using a hydraulic time step δt = 1 min over a total time duration t_max = 60 min. The 1-min time step was considered small enough with respect to the free surface area of the reservoirs (Giustolisi et al. 2012), so that the steady-state approach used by EPANET in EPS gives accurate results. Due to the fact that turbines were in fact modelled as valves, no efficiency-flow rate curve could be introduced in EPANET, so the calculations of power and energy production/consumption were performed in MS Excel, based on the flow rate and head values obtained from EPANET. This procedure proved useful when investigating only the influence of the initial level in the upstream reservoir (tank R1) on the HPD operation (Dunca et al. 2018). To provide a more complete analysis of the influence of the initial levels in the reservoirs on the production/consumption of energy in the HPD, the above-mentioned procedure was inadequate. In a first attempt to solve this problem, the EPANET model was imported in Bentley WaterGEMS, a commercial computer software that has a turbine object already defined. Unfortunately, this did not solve the problem, as for the turbine a head-flow rate curve can be inserted, but the efficiency of the turbine is considered constant, which would lead to the same procedure as the one used in EPANET for the computation of power and energy production. Finally, a GNU Octave script file was developed for this specific application, as presented in the sequel.

Tanks R1 ÷ R3 are set to be filled and/or emptied with respect to their capacity curves (volume curves), denoted z_R1 = f(V_R1), z_R2 = f(V_R2), and z_R3 = f(V_R3), plotted in Figure 3. The minimum level and maximum level setting values for each tank equal the minimum and normal retention levels of the selected reservoir. The initial level of each tank is set for each run, within the following ranges: for the upper tank R1, the initial level decreases, starting from the normal retention level of 1,255 m a.s.l., with a step of 1 m, down to 1,235 m a.s.l. (resulting in 21 values); for the storage tank R2, the initial level is kept constant and equal to the normal retention level of 1,007 m a.s.l. (being a small storage tank); for the lower tank R3, the initial level decreases with a step of 1 m, from the intermediate level of 780 m a.s.l., to the level of 760 m a.s.l. (resulting also in 21 values). Thus, for each of the first two operation scenarios, 441 runs can be performed, while the third scenario can be analysed upon 21 runs, yielding a total of 903 runs.

Selecting as initial conditions, at time t⁽⁰⁾ = 0, different initial water-level values z_R1(t⁽⁰⁾), z_R2(t⁽⁰⁾), and z_R3(t⁽⁰⁾) in tanks R1 ÷ R3, the updated water levels z_R1(t^(k)), z_R2(t^(k)), and z_R3(t^(k)) are computed at the end of each time step, i.e. at time t^(k) = (t^(k−1) + δt), for k = 1 ÷ 60, where t⁽⁶⁰⁾ = t_max.

The code allows running a steady-state hydraulic analysis on each time step δt, yielding a snapshot solution at the beginning of the time step, at time t^(k−1) (here k = 1 ÷ 60). So, the EPS involves computing successive snapshot solutions, connected from time t^(k−1) to t^(k) through mass balance computations. This procedure allows to decrease or increase the water volume V^(k−1) in any tank, and the water level z^(k−1) = f(V^(k−1)), upon the direction of the flow rate Q^(k−1). The volume variation from time t_k−1 to t_k can be approximated as δV^(k) = (V^(k) − V^(k−1)) = ±Q^(k−1)·δt, where the negative variation corresponds to the outflow (tank emptying), while the positive variation corresponds to inflow (tank filling).

A system of nonlinear equations can be derived for Gâlceag HPD, at time t^(k−1) (for k = 1 ÷ 60), according to each operation scenario (Figure 4). The system consists of energy balances between two successive nodes, and continuity equations at nodes, e.g. as described within the system (1) of Giustolisi et al. (2012), where one must consider that pumps run at a constant speed, and turbines are modelled by minor losses. Efficiency curves can be attached to the pumps to complete the cited equations system. Such a system of nonlinear equations can be solved in EPANET using the Global Gradient Algorithm (Todini & Pilati 1988).

Potential instabilities may appear in EPANET computations, when dealing with EPS and multiple tanks that are not nearly far enough (Todini 2011). In our model, tanks R1 ÷ R3 are quite far apart from one another, and their storage capacity is huge with respect to the water volume variation due to the inflow/outflow, so EPANET can provide trustworthy solutions (Dunca et al. 2018), but as already mentioned, each solution requests post-processing outside EPANET. The solutions presented in this paper are computed using the function fsolve of GNU Octave (Eaton et al. 2019), also available in MATLAB^® (MathWorks 2019).

To perform computations using the built-in function fsolve, the nonlinear system of equations attached to Gâlceag HPD will be simplified in the sequel, to reduce the number of equations. The system can be simplified, assuming that both turbines run with identical loads when the HPP is fully operational, and both pumps run with identical duty points when the PS is fully operational. To obtain the solution with fsolve, the nonlinear system will be inserted in a user-defined function, with two arguments (Georgescu & Georgescu 2014), namely a vector containing the unknowns of the problem, and a vector containing three parameters, namely the initial water levels in tanks R1 ÷ R3.

The system will be derived for the second scenario, related to simultaneous operation of HPP and PS, being the most complex one. Simplifications will be mentioned further for the first and third scenarios.

All four regression curves defined by (2) and (3) are displayed in Figure 2, together with the corresponding R² values, and the limits of the total uncertainty bandwidth.

The dependencies H_t = H_t(Q_t) from (2) and H_p = H_p(Q_p) from (3), considered at time t^(k−1), are inserted into the system (1). Due to the friction factor formula, the resulting system is highly nonlinear. The solution of (1), consisting of the flow rates through one turbine, ⁠, and one pump, ⁠, yields the flow rates on all five pipes at time t^(k−1).

Iterations continue up to k = 60.

As expressed in (15), within the pumped storage operation, the upper tank R1 is filled.

Due to the lack of experimental data, the uncertainties induced by the proposed mathematical model over the energy assessment cannot be estimated.

When the Gâlceag HPD model was firstly created and run in EPANET (Dunca et al. 2018), it became obvious that the pumped storage operation (scenario ③) is possible only over a limited time period since the capacity of the PS suction reservoir (storage tank R2) is far too small with respect to the capacity of the HPP upstream reservoir (upper tank R1) − the ratio between the active storage capacities of those reservoirs is of 1/160. As highlighted by Dunca et al. (2018), emptying the tank R2 by pumping, starting from the normal retention level (where the reservoir is full), will account only for about 1.3% of the volume of the upper tank R1 corresponding to the average annual minimum water level.

This is the reason for which all computations related to scenarios ② and ③ (where the PS is fully operational) were performed for a unique value of the initial water-level z_R2(t⁽⁰⁾) in the storage tank R2, namely the maximal (normal retention) level.

Results for computations started with the initial water level equal to the normal retention level in both HPP upstream and PS suction reservoirs, and equal to 770 m a.s.l. in HPP's downstream reservoir, show that (Dunca et al. 2018):

• the power delivered only by the turbines (no pumps working) is about 115 MW;

• the overall power delivered to the electric grid when all pumps and turbines are working (obtained as the difference between the power provided by the turbines and the power consumed by the pumps) is of about 113 MW;

• the power required to pump water from the PS suction reservoir to the HPP upstream reservoir is about 15.5 MW.

The relatively small difference between the first scenario (normal operation) and the second scenario (simultaneous operation) provided a clue as to what could be the real purpose of such a combined facility. A comparison among the values obtained for the overall energy production in the scenarios ① and ②, when the initial level z_R1 in the HPP upstream reservoir was varied between the normal retention level and the average annual minimum level, while the other initial levels were set to z_R2 = 1,007 m a.s.l. and z_R3 = 770 m a.s.l., is shown in Figure 5. It is clear from Figure 5, that for an initial water level in the downstream reservoir (770 m a.s.l.), for some values of the water level in the HPP upstream reservoir, the overall quantity of energy provided by the system when all pumps and turbines are working, exceeds the quantity that would be provided if only the turbines would work. This occurs for water levels in the HPP upstream reservoir between 1,252 and 1,238 m a.s.l. Differences are relatively small, so in order to increase the accuracy, we present in Figure 6, the gain in energy, obtained when pumps and turbines work together − i.e. the differences between the overall energy produced in scenario ② and the energy produced in scenario ① (Dunca et al. 2018). Obviously, positive values represent a gain in the overall energy delivered by the system. The maximum gain in energy is of about 2,650 kWh, obtained when the initial water level in the HPP upstream reservoir is of 1,242 m a.s.l.

The behaviour of the system can be explained as follows. On one hand, the flow rate through the turbines is increased when pumps and turbines work together (scenario ②), with respect to the case when only turbines are used (scenario ①). On the other hand, the flow through the pumps is also increased as the level in the HPP upstream reservoir decreases.

For some intermediate values of the initial level in the upper tank R1, the flow increase shifts the values of the efficiency of both turbines and pumps towards bigger values, leading to the more efficient production of energy in the turbines and more energy introduced by the pumps in the flow.

According to Figure 6, for the initial water level of 1,242 m a.s.l. in the tank R1, the gain in energy reaches a maximum. This energy gain does not necessarily mean that the system is more efficient.

In this respect, it may seem like a good idea to use the pumps when the initial water level in the HPP upstream reservoir is between 1,252 and 1,238 m a.s.l., as more energy is provided to the electric grid.

As already stated, 21 equally spaced values of the initial water level were set (successively and decreasingly) for each of the upper tank R1 and lower tank R3, starting from the normal retention level of R1, and from an intermediate level of R3.

The results obtained for the energy production E_HPP|₍₁₎ in Gâlceag HPP for scenario ① (only HPP in operation), and the global energy production E_HPD|₍₂₎ in Gâlceag HPD for scenario ② (both HPP and PS in operation), are presented in Figures 7 and 8, as 3D surfaces and corresponding 2D contour plots, with respect to the variation of the initial water levels z_R1 in the tank R1 and z_R3 in the tank R3. The resulting energy gain ΔE_HPD, defined by (13) for the entire HPD, to compare the simultaneous operation (scenario ②) and the normal operation (scenario ①), is presented in Figure 9, as 3D surface and attached 2D contour plot, as function of the initial water levels z_R1 and z_R3.

The variation of the energy gain ΔE_HPD with respect to the initial water level z_R1 in tank R1, is plotted in Figure 10, for three values of the initial level in tank R3, namely z_R3(t⁽⁰⁾) = {770; 765; 760} m a.s.l.

The difference between the water levels in the upper tank R1 and the lower tank R3 is the gross head H_g of the HPP. The computed data shows that the energy gain (Figure 9) is always positive for some specific values of the initial gross head, e.g. for H_g∈[469; 477] m, attached to the initial water levels z_R1(t⁽⁰⁾) ≤ 1,237 m a.s.l. and z_R3(t⁽⁰⁾) ≤ 766 m a.s.l. Other specific values of H_g attached to positive ΔE_HPD can be retrieved from the narrow white band on the right frame of Figure 9.

This paper presents the numerical modelling of a HPP coupled to a PS, in a complex HPD existing in Romania. Measured characteristic curves of the pumps and the hydraulic turbines were used in the model.

Three main scenarios were tested: ① only the HPP in operation, as a conventional HPP; ② both HPP and PS in operation (where the pumped water is discharged directly into the penstock, increasing the flow rate towards the HPP); and ③ only PS in operation (where water is pumped from the PS suction reservoir to supply the HPP upstream reservoir, through the HPP penstock and headrace). The initial water levels in the HPP upstream and downstream reservoirs were varied (over 20 m, by 20 equally spaced steps), starting from the normal retention level, and an intermediate level, respectively.

The simulations give the energy production and/or consumption in the studied cases. The results show that for some predefined ranges of water levels in the HPP upstream and downstream reservoirs, more energy is generated if turbines and pumps work simultaneously (scenario ②).

For a more complete analysis, further work is necessary to assess the partial operation of the HPD (e.g. only one turbine in operation and/or only one pump, or one turbine and both pumps).

Although the accuracy of the model is difficult to ascertain without measured values of energy production and consumption at the HPD, in close correlation to water levels at the reservoirs, the presented methodology, i.e. using on-site measured curves for pumps and turbines in a numerical model of a pipe network, yields usually better results than using catalogue curves (performances of hydraulic machineries decrease in time, while catalogue curves are determined for new items). Moreover, with the advent of renewable energies that have a more unsteady production in the energy mix, the role of hydropower as a regulator of the power system has increased. For such purposes, the analysis could prove useful.

The results can be used within a decision support system, to assess the overall operation of the HPD upon the water availability (water levels in reservoirs), in order to evaluate the amount of power that can be used for power generation.