Numerical analysis of an initial design of a counter-rotating pump-turbine

Renewable sources of energy are on the rise and will continue to increase the coming decades [1]. A common problem with the renewable energy sources is that they rely on effects which cannot be controlled, for instance the strength of the wind or the intensity of the sunlight. The ALPHEUS Horizon 2020 EU project has the aim to develop a low-head hydraulic pump-turbine which can work as a grid stabilising unit. This work presents numerical results of an initial hub-driven counter-rotating pump-turbine design within ALPHEUS. Computational fluid dynamics simulations are carried out in both prototype and model scale, for pump and turbine modes, and under steady-state and unsteady conditions. The results indicate that the initial design have a hydraulic efficiency of roughly 90 % in both modes and for a wide range of operating conditions. The unsteady simulations reveal a complex flow pattern downstream the two runners and frequency analysis show that the dominating pressure pulsations originates from the rotor dynamics. Given the promising high efficiency, this initial design makes an ideal platform to continue the work to optimise efficiency and transient operations further.


Introduction
Hydropower plays a key role to provide a stable and flexible electrical grid. With an increasing amount of energy from renewable sources, such as wind, solar, etc., the need of complementary controllable energy sources increases. By storing large amount of water when there is excess power in the grid, and later utilise the stored water when there is a lack of power, hydropower is a stabilising unit for the electrical grid [2,3]. The ALPHEUS (Augmenting Grid Stability Through Low Head Pumped Hydro Energy Utilization and Storage) Horizon 2020 EU project has the aim to develop a low-head to ultra low-head sea-water based hydraulic machine that can work both as an energy storage by pumping water to a reservoir, and to extract the energy by reversing the pump to a turbine. The goal is that the pump-turbine should have a round-trip efficiency of 70 -80 % and a switching time of about 120 s [4]. The final machine should have a scalable design with capacity ranges from 300 kW up to 10 MW depending on location. The aim is to have the peak efficiency at a head of 10 m, but heads down to 2 m are considered. Three pump-turbine concepts are to be investigated, a counter-rotating hub-driven, a counter-rotating rim-driven, and a positive-displacement alternative. The latter as a fish friendly option. A 30th IAHR Symposium on Hydraulic Machinery and Systems IOP Conf. Series: Earth and Environmental Science 774 (2021) 012066 IOP Publishing doi: 10.1088/1755-1315/774/1/012066 2 rigorous optimisation process will be applied to maximise the round-trip efficiency for a wide range of operating conditions.
The counter-rotating propeller concept dates back to the early 19th century and is traditionally associated with aeronautical and maritime propulsion systems [5]. Furukawa et al. [6] claims that an axial counter-rotating pump can have the advantages of being more compact, having a lower rotational speed, preferable head characteristics, and a wider range of operations at a high efficiency than a conventional single rotor pump configuration. The counter-rotating pump-turbine in a hydropower context has recently been investigated both numerically and experimentally in a series of studies [7,8]. The work has focused on increasing the performance of a model counter-rotating pump-turbine.
This work presents a numerical analysis, of an initial design of a low-head counter-rotating hub-driven pump-turbine, in prototype and model scale. The blade geometries, shown in Figure  1a, were designed by the Advanced Design Technology Ltd (ADT) company. The rotational speed and direction of each runner is controlled individually, maintaining a close to axial flow both before and after the machine, in both pump and turbine modes. Computational Fluid Dynamics (CFD) simulations are made in steady-state on a stage model (blue and red parts) seen in Figure 1b. The steady-state simulations are done in both prototype and model scale. Unsteady simulations are carried out in model scale on the full domain, gray area (with extensions) in Figure 1b. The simulations made in model scale will be validated with experimental data that will be generated as part of ALPHEUS.
(b) Single blade-passage (stage) computational domains, colored blue and red (for runners 1 and 2), and center part of the full computational domain, colored gray. The full domain is extended with straight pipes at both ends (not shown).

Design of the pump-turbine
The basis of a well-designed hydraulic counter-rotating pump-turbine is that the machine has axial flow before and after the two runners. The principle is that the first runner swirls the flow and the second runner de-swirls the flow, making the flow axial after the second rotor. The initial design is derived in prototype scale with ADT's TURBOdesign suite. It is a software with an inverse design method that uses inviscid flow to quickly generate and evaluate blade geometries based on a number of input parameters. The main constraints to determine the initial blade design are power, head, flow rate, work coefficient at leading and trailing edge, and the blade loading coefficient at four streamwise positions [9]. The initial blade design shown in Figure 1a has a diameter of roughly 6 m in prototype scale. We here consider cases where the first (blue) runner has a rotational speed of 50 RPM and the second (red) runner rotates in the opposite direction at 90 % of the speed of the first runner. The initial prototype design is 3 simulated with CFD as two stages with one blade per runner in each stage, to verify that the performance at prototype scale is in accordance with the requirements set by ALPHEUS.
The numerical results are later to be validated against experimental data, generated in the hydraulic laboratory at TU Braunschweig, on an optimised counter-rotating pump-turbine design. To accommodate with requirements in the laboratory an appropriate model scale is a runner with a diameter of 27 cm. By applying scaling laws for hydraulic machines on the prototype scale, suitable operating conditions are found in model scale [10]. In model scale, the first runner has a rotational speed of 1453 RPM in pump mode and 832 RPM in turbine mode, the second runner rotates at 90 % of the speed of the first runner in each mode.

Computational Domain
The steady-state simulations of the prototype and model scales are done over the two runner stages, shown as blue and red in Figure 1b

Simulation set-up and mesh
Numerical CFD simulations are made using the computational domains shown in Figures 1b and 2. The flow is from left to right in pump mode and from right to left in turbine mode. Steady-state simulations are computed over the two stages with one blade per runner and assuming an infinite hub. MRF zones, a mixing plane interface, and periodic GGI interfaces are used to mimic the complete counter-rotating runners. The unsteady simulations include a full 360°domain, two rotating runners, support struts, a finite hub, and contraction/expansions that are extended by straight pipes. The total length of the full domain is 12.8 times the diameter of the runners. Three AMI interfaces are adopted, before, between, and after the runners. FOAMextend 4.1 and OpenFOAM-v1912 are used for the stage and full simulations, respectively. The mixing plane is only available in the former version. The steady-state simulations made in prototype scale are compared with results produced by the proprietary CFX CFD code. In all the simulations, a total and static pressure is applied for inlet and outlet individually, the flow rate is thus calculated by the solver. This type of boundary conditions have successfully been verified against experimental test data for a model scale Francis hydro turbine [11].
The computational mesh is block-structured for the stage simulations, and a combination of block-structured and unstructured for the full simulations. The mesh for the full simulations has a number of refinement zones to capture the flow field in detail in areas of interest. The mesh is finest at the runners and is gradually coarser further away from the runners. The blockstructured mesh of the runners is illustrated in Figure 3a

Numerics in OpenFOAM
The numerical simulations utilise the finite volume method to discretise and solve the incompressible Navier-Stokes equations on a computational mesh [12]. To account for turbulence and closure of the Navier-Stokes equations the two-equation eddy-viscosity k-ω SST RANS model is applied for the steady-state, stage, simulations [13]. The model has been proved to capture the flow field to a great extent for hydro turbines [7,14]. For the full unsteady simulations the additional Scale Adaptive Simulation (SAS) modifications are applied to the k-ω SST model. The supplementary SAS modifications is an improved version of the RANS model for unsteady flows. It resolves the larges scales by allowing the flow to turn unsteady by decreasing the turbulent viscosity where appropriate [13]. The convection terms of the momentum equations are discretised using the second-order linear upwind scheme and the Linear-Upwind Stabilised Transport (LUST) scheme, for the stage and full simulations respectively. The convection terms of the k and ω equations are discretised with the first-order upwind scheme. Temporal discretisation is made with the second-order backward scheme [15]. The time step is chosen to 5 × 10 −5 s, as a trade-of between accuracy and the total computational time. The Courant number is below 2.2 in 90 % of the domain, and the maximum Courant number is less than 35.
For the stage simulations the BiCGStap solver is applied with the DILU preconditioner for all variables. This is done since a solver for asymmetric matrices is required when employing a mixing plane interface. The full simulations use the smoothSolver along with Gauss seidel for velocity and transported variables. For pressure the GAMG solver is utilised together with the DICGaussSeidel smoother [15,16]. To couple pressure and velocity the SIMPLE algorithm is used for the stage simulations and the PIMPLE algorithm for the full simulations. The later is a combination of the PISO and SIMPLE algorithms allowing larger time steps [15,17]. The PIMPLE loop is set to operate with two inner loops, one non orthogonal corrector, and a maximum of ten outer corrector loops or until convergence is reached within the time step.

Post-processing
The engineering quantities that are compared between the stage and full simulations are given by [10]: Where P is power, T is torque, ω is rotational speed, H is head, ∆p is total pressure drop over the runners, ρ is density, g is gravity, η is efficiency, and Q is the flow rate. Subscripts R, P and T stand for runner, pump mode, and turbine mode, respectively. The density of water is ρ = 998 kg/m 3 and the gravity is g = 9.81 m/s 2 .

Prototype scale
To ensure that the required performance is according to ALPHEUS, and to show that the results with CFX and FOAM-extend-4.1 are consistent, the results are compared in prototype scale as shown in Figure 4. The solver settings and boundary conditions for CFX are as far as possible similar to those in FOAM-extend. By comparing the efficiency in pump mode as a function of flow rate in Figure 4a it is found that the two softwares produce very similar results. Small differences are noted when comparing the power and head in Figures 4b and 4c. The reason for this may be due to differences in solver design, numerical schemes, mesh, wall functions, etc. The differences are, for the time being, considered small, and the results from the two solvers are sufficiently similar. Lacking experimental validation data, the similar results from two different CFD solvers makes the results trustworthy.
The ALPHEUS project has the goal to design a low-head pump-turbine with a round-trip efficiency of 70 − 80 % and a power of roughly 10 MW. The results for the initial design in prototype scale surpass the requirements according to the performance presented in Figure 4. Both the best efficiency, and the wide range of high efficiency are promising given that this is the initial design before any optimisation.

Model Scale
The model scale is derived based on requirements of the hydraulic laboratory at TU Braunschweig, where validation experiments will be done at a later stage. Figure 5 shows the computed efficiency, power, and flow rate, as functions of head for both pump and turbine modes and for both the stage (steady) and full (unsteady) simulations. It is shown that the results of the stage and full simulations coincide to a satisfying extent. Some of the differences originate in the difference in the computational domains. The full simulations use OpenFOAM-v1912, and the stage simulations utilise FOAM-extend-4.1. This is due to the fact that the mixing plane interface, required for the stage simulations, is only available in FOAM-extend. Note that the full simulations are performed for one operating condition for each mode, while the stage simulations are done for a range of operating conditions. By comparing the results in prototype scale, Figure 4, and model scale, Figure 5, one also has a third verification that the results are consistent between solvers and mesh resolution. This is because in prototype scale FOAM-extend-4.1 is verified against CFX on a stage model, and in model scale the full simulations done with OpenFOAM-v1912 is compared to stage simulations made with FOAM-extend-4.1. The mesh used for the full simulations is slightly coarser than for the stage simulations, but the results are almost identical. Hence, one can conclude that the results are mesh and solver consistent. In total three different meshes have been evaluated, one for each solver and simulation type.

Flow investigations of the full simulations in model scale
The results from the unsteady simulations, shown in Figure 6, resolve the unsteady wakes of the struts and blades, and the interaction between the different parts of the machine. The downstream runner cuts through the wakes of the upstream runner, causing a complex flow pattern after the pump-turbine. The wakes and vortex shedding from the struts also contribute to the rotor-stator-interaction, and reveals that the outgoing flow is rather axial in pump mode, while there is a slight remaining swirl in turbine mode. The support strut upstream the runners cause a wake region in front of the runners. This has a rather small impact on the performance of the machine since the engineering quantities, shown in Figure 5, are consistent between the stage and full simulations. The stage simulations do not include any support struts.  As can be seen in Figure 7, the time averaged velocity profiles on the lines defined in Figure  2, downstream the runners for each mode, reveal interesting characteristics for the two modes. The axial velocity component, u z , is quite asymmetrical in the pump mode, Figure 7a, and the radial and tangential components, u r and u θ , fluctuates around zero. The shape of the axial component is caused by the position of the support struts in relation to the location of the velocity line. Because there is practically no swirl, there is noting that alters the axial flow profile. The small swirl component is a proof that the machine is operating at a high efficiency as the second runner has managed to effectively de-swirl the flow caused by the first runner.
The velocity profiles downstream the runners in turbine mode has a different behavior than in pump mode. The axial velocity in Figure 7b is more symmetrical over the center line despite that the struts are positioned in the same manner as for the pump mode. This is explained by the fact that the tangential component, u θ , is everywhere negative, which causes the axial velocity profile to rotate tangentially in a counterclockwise direction and the wake of the strut cannot be seen on the probe line. The remaining swirl after the runners yields a less efficient machine since the upstream runner does not produce enough angular momentum for the downstream runner to effectively de-swirl.   The axial, u z , radial, u r , and tangential, u θ , velocity components are defined by the coordinate system shown in Figure 2. The tangential velocity is positive according to the right hand rule around the z-axis, and the radial velocity is positive outwards.
A frequency analysis of the signals from the pressure probes before, between, and after the runners, for both pump and turbine modes, is shown in Figure 8. Pressure probe P3 is located between the runners, while probes P1 and P5 are located before or after the runners depending on mode, as shown in Figure 2. All the dominating frequencies correspond to various combinations of the runners rotational speed and the number of runner blades. Recall that in pump mode the first runner has a rotational speed of 152 rad/s or f R1 = 24.2 rev/s (Hz), and the second runner rotates at 90 % the speed of the first runner, f R2 = 21.8 Hz. The first runner has eight blades and the second seven blades. The blade passing frequency of the first runner is determined by 8 × f R1 = 193 Hz, and the corresponding value for the second runner is 7 × f R2 = 152 Hz. According to Figure 8a, the frequency caused by runner two (downstream) has a larger impact on the pressure pulsations between the runners than the upstream runner. This is indicated by the higher peak at 152 Hz than at 193 Hz. The higher peak is caused by the stagnation pressure from the leading edge of the second runner that has a larger impact when it cuts the wakes from the first runner, than the wakes from the first runner itself. There is also a distinct peak at the blade passing frequency correlated to the first runner, 193 Hz, that is caused by the wakes of the first runner.
In turbine mode the rotational speed is lower. Runner one (downstream) rotates at 87.1 rad/s, with a frequency of f R1 = 13.9 Hz, and the upstream runner two has a rotational frequency of f R2 = 12.5 Hz. The dominant frequencies between the runners are in turbine mode 8 × f R1 = 110 Hz, and 7 × f R2 = 87 Hz. The two peaks, shown in Figure 8b, of 110 and 87 Hz have almost identical magnitude. This means that the wakes from the upstream runner and stagnation pressure from the downstream runner are equally responsible for the pressure pulsations between the runners in turbine mode and under these operating conditions.
(a) Pump mode, P1 is upstream, and P5 downstream the runners.
(b) Turbine mode, P5 is upstream and P1 downstream the runners. Figure 8: Frequency analysis of pressure probes before, between, and after the runners. P3 is located between the runners.
Most of the other peaks present in Figure 8 originate from harmonics of the blade passing frequencies from the runners. It is noted that at P1 and P5 the frequency of 346 Hz and 198 Hz for pump and turbine modes respectively have a small peak. Those frequencies originate from the sum of the first and the second rotor blade passing frequencies, 8 × f R1 + 7 × f R2 . Hence, combinations of the main frequencies also plays a role for the pressure pulsations in the system.

Conclusions
The present study has demonstrated a CFD analysis of an initial design of a counter-rotating hub-driven pump-turbine, as part of the ALPHEUS Horizon 2020 EU project. Given that this is the initial design and that the hydraulic efficiency is about 90 % for a wide range of operating conditions in both pump and turbine modes, the results are promising for the ALPHEUS project. The favorable performance provides an ideal platform to optimise the design and operation of the machine further. Similarities in result between the different solvers (CFX, Foam-extend 4.1, and OpenFOAM-v1912), mesh sizes, and simulation type (stage or full) confirms that the results are to a great extent numerically consistent and trustworthy. Later in the ALPHEUS project, experimental investigations of an optimised model scale pump-turbine will be performed and used to validate the simulation strategy.
The full unsteady simulations revealed that there is a complex flow pattern produced by the runners, as the downstream runner cuts the wakes of the upstream runner. The placement of support struts have a clear impact on how the flow field develops after the runners, as the support struts cause large unsteady wakes and vortex shedding. The pump-turbine is less dependent on the support struts before the machine, judging by the similarity of results between the stage and full simulations. Frequency analyses of the pressure fluctuations at the runners revealed that the dominating pressure pulses originate from the rotor dynamics caused by the blade passing frequency. The largest pressure pulsations are obtained between the runners, but some peaks are recognised before and after the runners as well. This initial work provides a solid platform to continue investigating unsteady and transient operations of the counter-rotating pump-turbine.