## Abstract

Previous studies suggested variable speed operation (VSO) of Francis turbines as a measure to improve the efficiency at off-design operating conditions. This is, however, strongly dependent on the hydraulic design and, for an existing turbine, improvements can be expected only with a proper redesign of the hydraulic surfaces. Therefore, an optimization algorithm is proposed and applied to the runner of a low specific speed Francis turbine, with an optimization strategy specifically constructed to improve the variable speed performance. In the constrained design space of the reference turbine, the geometry of the replacement runner is parametrically defined using 15 parameters. Box–Behnken method was used to populate the design space with 421 unique samples, needed to train fully quadratic response surface models of three characteristic efficiencies defined by the proposed objective function. Computational fluid dynamics (CFD) was used to calculate the responses for each sample. The parametric study showed that the anticipated variation of the shape of the hill chart, needed to improve the variable speed performance of the turbine, is limited within a narrow range. The presented method is general and can be applied to any specific speed in the Francis turbine range, for both synchronous speed and variable speed optimization tasks.

## 1 Introduction

In Norway, as well as throughout most of Europe, the hydropower sector is facing an increasing demand for additional operational flexibility. Supporting ambitious directives and targets set for sustained deployment of nondispatchable renewable energy and reduction of carbon emissions [1–3], hydraulic turbines are now operated far more roughly than before and with extended periods of off-design operation [4–6]. Most of the hydraulic turbines are Francis turbines, and are operated at constant rotational speeds. In such case, the operating conditions are dictated by the available head and the discharge through the turbine, eventually leading to a mismatch of the rotational speed and more losses [7,8].

Owning to the pioneering work done in Japan on variable speed operation (VSO) of reversible pump-turbines [9–12], several researchers have explored the idea to apply this technology to improve off-design efficiency of conventional Francis turbines [13–15]. However, both old and recent studies have shown that for low specific speed Francis turbines operated close to their rated head, such as most of the installed turbines in Norway are, enabling VSO might result in minor efficiency improvements [13,16,17]. This outcome is found to be closely related to the hydraulic design of the runner [16], suggesting that improvements might be expected with a proper redesign/optimization.

With the increased performance of today's computers, as well as the advances in modeling and numerical algorithms done over the past decade, the flow inside the complex geometry of Francis turbines can be numerically analyzed and accordingly optimized. Risberg et al. [18] used computational fluid dynamics (CFD) to train approximation model of an optimization objective that was used for an inexpensive search after the optimal configuration in the design space. Lyutov et al. [19] and Nakamura and Kurosawa [20] combined CFD with population-based genetic algorithm to evolve initial design population toward the optimal design in several generations/offspring. Thum and Schilling [21] used methods to estimate the gradients for the objective functions that are then used for ascend/descent “walk” to approach the optimal solution. So far, these methods were successfully used to improve the performance of synchronous speed Francis turbines, however, without any evidence for VSO applications yet.

To investigate the possibility for improvement of the VSO performance of Francis turbines, an optimization strategy based on numerically trained metamodels is proposed and studied. The VSO performance is regarded from the efficiency point of view only, where the speed factor $nED$ is adjusted to ensure operation along a specific curve in the $nED\u2212QED$ plane, which corresponds to the hill profile as projected on the $\eta h\u2212QED$ plane. Other operational aspects, such as cavitation and/or pressure pulsations, are not considered in the optimization procedure; however pressure distributions were manually checked once an optimal design is found. The presented case study is a model of a low specific speed turbine where a replacement runner was optimized with no changes on the remaining parts of the turbine. The geometry of the replacement runner is described parametrically using 15 free parameters, resulting in total of 421 design samples needed to estimate the coefficients of fully quadratic response surfaces. For each design sample, CFD was used to calculate the hill chart in the zone where the efficiency curve for VSO can be estimated.

The remaining of the paper is structured as follows: Sec. 2 describes the theoretical background regarding the optimization strategy for improvement of the variable speed performance. A case study is presented in Sec. 3, with details on: the VSO performance of the reference turbine, the parameterization method used, the CFD details for the numerical calculation of the objectives and the analysis using the response surface models. A discussion regarding the possibilities for VSO improvement is presented in Sec. 4 and the final conclusions are disclosed in Sec. 5.

## 2 Theoretical Background of the Optimization Procedure

Regarding the hydraulic efficiency of the turbine, two operating scenarios can be defined where speed adjustments might be necessary: (1) operation under large head variation; and (2) off-design operation close to the design head of the turbine. In the former scenario, there is no significant challenge to improve the efficiency of an existing turbine with VSO, i.e., the optimal value of the speed factor $nED=n\xb7D2/g\xb7Hn$ can be achieved by proportionally matching of the speed $n$ to the net head $Hn$. However, due to equipment and reservoir restrictions, not many sites exist with such operating conditions. For the latter, and perhaps more common scenario, efficiency improvement with VSO appears to be limited by the hydraulic design of the turbine, or at least that is the case for low specific speed machines. Therefore, the optimization task hereafter is focused on improving the off-design efficiency benefits from VSO close to the design head of the Francis turbine.

### 2.1 Optimization Strategy for Variable Speed Operation.

*objective function*. The object/process that is subject to the optimization is normally described with a set of inputs, or

*free parameters*, where the goal of optimization is to minimize/maximize an objective function. In most cases, the

*design space*is constrained, meaning that each input parameter has a restricted range of their values. Typically, an optimization will be performed as a tradeoff between two or more objective functions, also called

*multi-objective optimization*. For optimization of Francis turbines, the hydraulic efficiency $\eta h$ at one or more operating conditions (

*Pareto optimal*) is usually set to be maximized, given by the general mathematical definition

where $y$ is the optimization objective, $n$ is the number of free parameters, $x$ is the design vector of the free parameters. The superscripts $l\u2009\u2009\u2009and\u2009\u2009u$ denote the lower and upper values, respectively, in the constrained range of the *i*th design parameter.

In contrast to the synchronous speed operation (SSO), VSO provides continuous variation of the rotational speed during operation of the turbine. This is typically achieved by either using full size frequency converters (FSFC) or doubly fed induction machines to enable decoupling between the power grid and the production unit. Both technologies have certain limitations and introduce additional losses and complexity in the power generation process [17]. In some cases, especially for FSFC operating close to the peak efficiency of the turbine, it is more efficient to switch back to SSO and avoid the additional losses. Due to this, one assumption is that a runner optimized for VSO must have a synchronous speed as a design value, i.e., a predefined location of the best efficiency point (BEP) in the performance hill chart. Hence, the first objective (i.e., goal) of the VSO optimization is to—maintain the peak efficiency of the turbine and the BEP position in the hill chart calculated with CFD. In this following sections, this efficiency is noted as $\eta BEP$.

Turbine operation at off-design conditions can be classified into two categories, one being part load (PL) when the operating value of $QED$ is lower than the optimal, the other being high load (HL) when the operating value of $QED$ is higher than the optimal. On Fig. 1(a) shown are the positions of HL and PL operating points for SSO in the hill chart area of a typical Francis turbine. The *zero-swirl* curve represents the location with near optimal flow conditions at the outlet of the runner, while the *zero-incidence* curve represent the location with near optimal flow at the inlet of the runner. Moving away from those two curves will obviously worsen the flow conditions and reduce the efficiency of the turbine. The peak efficiency, i.e., the BEP, is typically found close to the intersection point between these two curves, where the flow conditions at both the inlet and outlet of the runner is optimal. The efficiency at HL and PL will be reduced by the $\Delta \eta $ losses with the subscripts $I,S$ referring to the incidence and swirl losses, respectively. At this point, it becomes obvious that the combination of $\Delta \eta I+\Delta \eta S$ at any operating point (i.e., power output) along the SSO line may not necessarily be the most optimal one and that a better tradeoff might exist at a different $nED$ value. In fact, for any Francis turbine, there exists an operating curve in the $nED\u2212QED$ plane that corresponds to the hill profile shown on Fig. 1(b) and this is the highest possible efficiency that can be achieved only with adjustable speed configuration. In other words, VSO seeks for the *perfect balance* between the partial efficiency losses for a constant $QED$ and this is dependent on the dominant losses in a turbine, i.e., the specific speed and the hydraulic design.

At this point, it is reasonable to conclude that an improved VSO performance means a flatter hill profile in the $\eta \u2212QED$ plane. To achieve that, two more objectives are defined as: maximize the efficiency for two discharge factor levels $QED(PL)$ and $QED(HL)$, irrespective of the speed factor $nED$. For that purpose, the minimum and maximum values in the operating range of $QED$ can be used. In the following sections, these two efficiencies are noted as $\eta PL$ and $\eta HL$ (see also Fig. 1(b)), and for the case study presented in Sec. 3, the two discharge factor levels $QED(PL)$ and $QED(HL)$ are set to be at 50% lower and 20% higher than the optimal $QED(BEP)$, respectively. The VSO optimization then seeks for designs that have at least the same BEP efficiency as the SSO representative, while simultaneously increasing the HL and PL efficiencies defined above. This effect has an explicit correspondence to a hill chart *stretching* in a suitable direction which is expected to improve the VSO performance, as it is suggested in Ref. [13] as well.

### 2.2 Response Surface Modeling.

*goodness-of-fit*must be checked and controlled to confirm that the accuracy of the approximated model is reliable for the optimization purpose. Two different criteria are typically used, the

*root-mean-squared error*$\sigma e$ defined as

where $y\xaf$ is the mean of the simulated output parameters.

### 2.3 Population of the Design Space Using Design of Experiments.

Appropriate estimation of the regression coefficients is crucial for successful application of the method described in Sec. 2.2. Therefore, design of experiments methods are typically used to populate the design space and reduce the amount of computational resources needed [22,23]. The efficiency of the RSM fitting relies on the task of getting the maximum amount of information with minimum number of design samples. The idea here is to choose $m\u2265p$ samples that are well distributed within the entire design space and can provide enough information on the expected nonlinearities in the response. Normally, the range of each parameter will be discretized in three levels and this is advantageous over the cheaper two levels since it permits estimation of the quadratic and interaction terms of the RSM. While one might be tempted to use either *three-level full factorial* or *random* sampling distribution, it is shown in Ref. [22] that the efficiency will be quite poor when the number of free parameters is $n\u22655$. Additionally, there is a lack of *rotatability* in the three-level full factorial designs, where the variance of the predicted values $y$ is also a function of the direction the point lies from the center point.

Since the number of simulations for the three-level full factorial designs grows rapidly with the number of free parameters, i.e., $m=3n$, a three-level incomplete factorial design is proposed in Ref. [24]. Commonly known as the *Box–Behnken designs (BBD)*, this sampling method has been specifically constructed for the use together with second-order RSMs and is able to produce *nearly rotatable* samples. A comparison between BBD and full factorial designs is shown on Fig. 2, where the number of free parameters is $n=3$. The construction of the design BBD matrix requires $m=2n(n\u22121)+CP$ number of samples, where $CP$ is the number of center points (sometimes more than one is needed for uniform estimate of the prediction variance). Most importantly, the design avoids computations in the extreme conditions, i.e., the vertices of the cube on Fig. 2(a), where the numerical results might be erroneous or difficult to converge.

A wide palette of design of experiments methods exists, such as *central composite design*, *fractional factorial design*, *Latin hypercube*, and *quasi-random designs* [23]; however, for the purpose in this paper, BBD seems to be the most appropriate tradeoff for balanced accuracy versus CPU time.

## 3 A Case Study

To investigate the possibilities for hydraulic optimization of variable speed turbines, the presented idea from Sec. 2 is used to design a replacement runner of a relatively low specific speed Francis turbine ($nq=22.35$). As it is shown in Sec. 3.1, the observed efficiency difference between the SSO and VSO for the reference turbine is relatively small when operated at the design head, which makes a perfect ground for application and testing of the optimization procedure. The strategy is to explore the constrained design space of the reference runner in order to find a hydraulic design that has “stretched” hill chart which will improve the VSO performance.

### 3.1 A Brief Description of the Reference Turbine.

where $n\u2009(s\u22121)$ is the rotational speed, $\rho \u2009(kg\u2009m\u22123)$ is the density of the water at the operating temperature, $g\u2009(m\u2009s\u22122)$ is the local gravitational acceleration, $Hn\u2009(m)$ is the net head of the turbine, $Q\u2009(m3\u2009s\u22121)$ is the discharge through the turbine, $T\u2009(N\xb7m)$ is the torque of the runner, $\omega \u2009(rad\u2009s\u22121)$ is the angular velocity of the runner, and $D2\u2009(m)$ is the outlet diameter of the runner. The characteristic dimensions and performance values are given in Table 1. As can be seen from Fig. 4, the VSO performance of the turbine is very similar to that of the SSO at the optimal head, improving the hydraulic efficiency by no more than $0.5%$ in the entire operating range. The necessary speed variation for the VSO is within the range of $\xb110%$ of the optimal rotational speed and can be achieved with both FSFC and doubly fed induction machines methods. However, the VSO improvements of the hydraulic efficiency are not enough to compensate for the additional losses resulting from the variable speed devices and will obviously result in an overall reduction of the turbine efficiency. The situation is slightly better if the turbine is operated under larger variations of the operating head, where a reduction of the pressure pulsation amplitudes from the rotor–stator interaction (RSI) can be also expected [16].

Parameter | Value |
---|---|

Net head, $Hn$ | $12\u2009(m)$ |

Discharge, $Q*$ | $199\u2009(l/s)$ |

Rotational speed, $n*$ | $323\u2009(rpm)$ |

Speed factor, $nED*$ | $0.173\u2009(\u2013)$ |

Discharge factor, $QED*$ | $0.151\u2009(\u2013)$ |

Guide vane, $\alpha GV*$ | $9.8\u2009(deg)$ |

Peak efficiency, $\eta *$ | $92.6\u2009(%)$ |

Inlet diameter, $D1$ | $0.620\u2009(m)$ |

Outlet diameter, $D2$ | $0.349\u2009(m)$ |

Blade angle, $\beta 1(shroud)$ | $82.4\u2009(deg)$ |

Blade angle, $\beta 2(shroud)$ | $18.2\u2009(deg)$ |

Inlet width, $b1$ | $59.6\u2009(mm)$ |

Max blade thickness, $h1$ | $\u223c6.9\u2009(mm)$ |

Min blade thickness, $h2$ | $\u223c1.8\u2009(mm)$ |

Parameter | Value |
---|---|

Net head, $Hn$ | $12\u2009(m)$ |

Discharge, $Q*$ | $199\u2009(l/s)$ |

Rotational speed, $n*$ | $323\u2009(rpm)$ |

Speed factor, $nED*$ | $0.173\u2009(\u2013)$ |

Discharge factor, $QED*$ | $0.151\u2009(\u2013)$ |

Guide vane, $\alpha GV*$ | $9.8\u2009(deg)$ |

Peak efficiency, $\eta *$ | $92.6\u2009(%)$ |

Inlet diameter, $D1$ | $0.620\u2009(m)$ |

Outlet diameter, $D2$ | $0.349\u2009(m)$ |

Blade angle, $\beta 1(shroud)$ | $82.4\u2009(deg)$ |

Blade angle, $\beta 2(shroud)$ | $18.2\u2009(deg)$ |

Inlet width, $b1$ | $59.6\u2009(mm)$ |

Max blade thickness, $h1$ | $\u223c6.9\u2009(mm)$ |

Min blade thickness, $h2$ | $\u223c1.8\u2009(mm)$ |

The hydraulic design of the runner has a significant impact on the overall performance of the turbine [28], hence designing a replacement runner only, while keeping the remaining turbine parts fixed, is perhaps the simplest optimization task for VSO. The following sections are then focusing mainly on the search for an optimal runner design only.

### 3.2 Parametric Definition of the Runner.

*meridional projection*of the runner with several conformal maps of the

*stream surfaces*. To obtain the pressure and suction sides of the blade, a varying

*blade thickness distribution*is applied in both directions of the surface

*normal vectors*. In this paper, all these steps are parametrized using quadratic/cubic

*Bezier curves*, generally defined as [32]

where $p(t)$ is the position vector of a point lying on the curve, $P$ is the vector comprised of the control points, $Bin(t)$ are the *Bernstein polynomials*, $(ni)$ are the *binomial coefficients*, $t[0,1]$ is a parameter of the curve, $i[0,n]$ is an integer number of the sum, $n$ is the degree of the Bernstein polynomials, representing also the degree of the Bezier curve (i.e., $n=2$ for quadratic, $n=3$ for cubic, etc.).

*stream function*in polar coordinates, defined as a

*scalar field solution*to the partial differential equation in Eq. (10) [33]. The neglected influence from the blades will influence the resulting design, meaning that additional adjustments will have to be done on the design values of $Hn$ and $Q$ to correct the blade angles globally. Hence, for fixed rotational speed $n$, the head $Hn$ and discharge $Q$ become parameters of the design that are free to vary in a predefined range. Individual corrections of the blade angles at the outlet are also enabled for three streamlines

where $\Delta \beta (j)$ are the free parameters constrained in a specific range for each streamline.

where $\phi $ is the *wrapping angle* of the blade at any location and $dLs$ is a length along the streamline in the meridional view. To apply the equation, intermediate blade angle distributions are defined for the four streamlines using quadratic Bezier curves. Integrating from inlet to outlet of the blade, the total wrapping angles for each streamline can be defined parametrically $\phi SL(j)=\u222b121rtan\beta dLs$. The overall length of the blade is controlled by allowing $\phi SL1$ to be a free parameter within a constrained range, while the remaining $\phi SL(j)$ for $j=2\u20134$ are controlled by a leaning angle $\phi TE$ with a predefined TE shape as seen in the *top view* of the runner on Fig. 5. Additionally, a leaning of the LE is enabled through the free parameter $\Phi $ that rotates the streamlines in the circumferential direction, with its effect best seen on the side view of the runner on Fig. 5.

*modal shapes*with less than 5

*nodal diameters*by the

*fundamental*and

*second harmonic*of the

*RSI*according to Ref. [35]

where $m=1,2,3\u2026$ is the harmonic number, $Zg=28$ is the number of guide vanes, $n=1,2,3\u2026$ is an arbitrary integer, and $k$ is the number of nodal diameters with negative values referring to a rotation of the modal shape in the opposite direction compared to the runner. In Table 2, given are the lowest numbers of nodal diameters that can be excited by each of the first five harmonics of the RSI, for the typical range of number of blades. The most appropriate number of blades is selected to be $Zr=17$. The free parameters used to define the geometry of the runner are summarized in Table 3, with the corresponding range of variation for each parameter used for the optimization procedure. The coordinates for the meridional view are given in a global coordinate system that is positioned to have the $z=0$ plane passing through the middle of the channel with width $b1$.

$abs(kmin)$ | $Zr=16$ | $Zr=17$ | $Zr=18$ | $Zr=19$ |
---|---|---|---|---|

$m=1$ | 4 | 6 | 8 | 9 |

$m=2$ | 8 | 5 | 2 | 1 |

$m=3$ | 4 | 1 | 6 | 8 |

$m=4$ | 0 | 7 | 4 | 2 |

$m=5$ | 4 | 4 | 4 | 7 |

$abs(kmin)$ | $Zr=16$ | $Zr=17$ | $Zr=18$ | $Zr=19$ |
---|---|---|---|---|

$m=1$ | 4 | 6 | 8 | 9 |

$m=2$ | 8 | 5 | 2 | 1 |

$m=3$ | 4 | 1 | 6 | 8 |

$m=4$ | 0 | 7 | 4 | 2 |

$m=5$ | 4 | 4 | 4 | 7 |

Parameter | Range |
---|---|

$r$ of point $P1\u2009(mm)$ | $[170,\u2009200]$ |

$r$ of point $P2\u2009(mm)$ | $[50,\u200980]$ |

$z$ of point $P2\u2009(mm)$ | $[\u221280,\u2009\u221240]$ |

$z$ of point $P3\u2009(mm)$ | $[\u2212250,\u2009\u2212180]$ |

$r$ of point $T1\u2009(mm)$ | $[50,\u200980]$ |

$z$ of point $T2\u2009(mm)$ | $[\u2212171,\u2009\u2212151]$ |

$\phi SL1\u2009(rad)$ | $[1.1,\u20091.6]$ |

$\Phi \u2009(deg)$ | $[\u22121,\u20093]$ |

$\phi TE\u2009(deg)$ | $[0,\u200930]$ |

$\Delta \beta 2(1)\u2009(deg)$ | $[0,\u200912]$ |

$\Delta \beta 2(2)\u2009(deg)$ | $[0,\u20096]$ |

$\Delta \beta 2(3)\u2009(deg)$ | $[0,\u20093]$ |

$Hn/Hn*\u2009(%)$ | $[0.95,\u20091.05]$ |

$Q/Q*\u2009(%)$ | $[0.75,\u20091]$ |

$h1\u2009(mm)$ | $[15,\u200925]$ |

Parameter | Range |
---|---|

$r$ of point $P1\u2009(mm)$ | $[170,\u2009200]$ |

$r$ of point $P2\u2009(mm)$ | $[50,\u200980]$ |

$z$ of point $P2\u2009(mm)$ | $[\u221280,\u2009\u221240]$ |

$z$ of point $P3\u2009(mm)$ | $[\u2212250,\u2009\u2212180]$ |

$r$ of point $T1\u2009(mm)$ | $[50,\u200980]$ |

$z$ of point $T2\u2009(mm)$ | $[\u2212171,\u2009\u2212151]$ |

$\phi SL1\u2009(rad)$ | $[1.1,\u20091.6]$ |

$\Phi \u2009(deg)$ | $[\u22121,\u20093]$ |

$\phi TE\u2009(deg)$ | $[0,\u200930]$ |

$\Delta \beta 2(1)\u2009(deg)$ | $[0,\u200912]$ |

$\Delta \beta 2(2)\u2009(deg)$ | $[0,\u20096]$ |

$\Delta \beta 2(3)\u2009(deg)$ | $[0,\u20093]$ |

$Hn/Hn*\u2009(%)$ | $[0.95,\u20091.05]$ |

$Q/Q*\u2009(%)$ | $[0.75,\u20091]$ |

$h1\u2009(mm)$ | $[15,\u200925]$ |

### 3.3 Description of the Numerical Model and Estimation of the Numerical Error.

Numerical prediction of the performance is done using ansyscfx by solving the steady-state Reynolds averaged Navier–Stokes equations coupled with the standard $k\u2212\epsilon $ model for the turbulence quantities. The calculation is done for a broad range of operating points using reduced/simplified geometry of the turbine, which is comprised of: (1) the full guide vanes cascade, (2) the entire runner with excluded friction losses and leakage losses through the labyrinth seals, and (3) the entire draft tube. The spiral casing and the stay vanes are excluded from the simulated domain and are assumed to have an insignificant impact on the off-design performance trends of the turbine. However, the head losses in the excluded passages are accounted for in the calculation of the net head of the turbine and are modeled as a function of the discharge over the entire simulated range, i.e., $\Delta pSC+SV=k\xb7Q2$. Therefore, any secondary flows and possible nonuniformity of the flow field originating from the excluded passages are not modeled, meaning that uniform angle of the velocity vector is imposed on the guide vanes inlet surface in both circumferential and spanwise directions. The pipe constant $k$ was determined by simulating the BEP operating point for the entire geometry, including the spiral casing and stay vanes, and this value is later used as a constant to calculate over the entire hill chart. For all operating points in the hill chart, total-pressure-inlet boundary condition with prescribed flow direction in cylindrical coordinates is used at the inlet of the domain, while static-pressure-outlet condition is used for the outlet of the draft tube. No-slip wall condition is used for all wetted surfaces, while “frozen rotor” interface is used to connect the nonconformal meshes from both stationary and rotating domains [36]. A high-resolution (second-order) advection scheme is used for all equations being solved. The length scale is evaluated using the cube root of the domain volume and then multiplied with an iteration-dependent timescale factor (linearly increasing). The number of iterations per operating point was set to 600, achieving at least three (and in most cases up to five) orders of magnitude decrease in the RMS residuals for all equations being solved.

Setting for the discharge $Q$ to be an output from the simulation, the described numerical setup is well suited for simulation purposes when the performance of the turbine is not known. In contrast to the numerical model where the discharge is prescribed from prior experimental data (e.g., mass-flow-inlet), the presented model gives somewhat slower and poorer convergence rate, especially for calculations of off-design operating conditions of the turbine. In order to estimate the spatial discretization error, the *Richardson extrapolation-based grid convergence index* was used [37,38]. Since the turbine efficiency is the focus, and due to the lack of conservation property in the definition of the grid convergence index, the discretization error is calculated for the efficiency only and not for the torque and discharge separately. Three structured hexahedral meshes with different resolution were generated (coarse, medium, and fine), initially starting from the coarsest mesh that can achieve acceptable mesh quality with $20\u2264y+\u2264200$ and then increasing the element count throughout a structured refinement in all directions (excluding the height of the cells attached to the walls). Out of the total element count for each mesh, approximately $\u224850%$ of the elements are in the runner, $\u224830%$ are in the draft tube, and the remaining $\u224820%$ are in the guide vanes ring. This ratio was found to be suitable in order to achieve: (i) as conformal as possible transition of the meshes on the connecting interfaces, (ii) balanced spatial resolution according to the contribution/influence of the domains on the overall performance of the turbine, and (iii) comparable convergence trends of all three domains.

The coarse mesh used is shown on Fig. 7. The systematic refinement was done with a mesh refinement factor $r\u22651.35$ in all three directions. The meshes were simulated for three operating points, namely: (1) $\alpha GV=6\u2009deg$, (2) $\alpha GV=10\u2009deg$, and (3) $\alpha GV=13\u2009deg$, and until asymptotic convergence of the RMS residuals was achieved or the change in the resulting parameters are less than $0.05%$ in every 200 new iterations. This requirement was achieved with less than 3000 iterations per calculation. With such setup, the iterative error is assumed very small and negligible, providing a good estimate of the apparent order of the method $p$. The results are summarized in Table 4. Considering the low discretization error that the fine mesh has, it can be assumed that a fully mesh independent solution over the entire hill chart has been achieved. The medium mesh gives balanced results in terms of accuracy versus computational effort. However, for the optimization purposes in this paper, using the fine and the medium mesh turns out to be very expensive and lengthy, especially when the number of parameters being optimized is more than ten. Therefore, in order to limit the computational effort to a reasonable level, as well as to achieve an efficient parallelization by fitting a single calculation on one cluster node only, the coarse mesh was selected for further usage. In that case, the estimated computational workload for the entire optimization process is approximately $7\xd7\u2009105$ CPU hours.

$\alpha GV=6\u2009deg$ | $\alpha GV=10\u2009deg$ | $\alpha GV=13\u2009deg$ | |
---|---|---|---|

$N1,\u2009N2,\u2009N3\u2009[106]$ | 23.19, 9.09, 3.72 | 23.04, 9.42, 3.77 | 23.05, 9.42, 3.71 |

$r21,\u2009r32$ | 1.3556, 1.3569 | 1.3473, 1.3568 | 1.3474, 1.3641 |

$\eta 1,\u2009\eta 2,\u2009\eta 3$ | 92.478, 92.345, 91.891 | 95.133, 95.097, 94.784 | 94.070, 94.081, 93.838 |

$p$ | 4.025 | 7.134 | 10.215 |

$\eta ext21$ | 92.534 | 95.137 | 94.069 |

$GCIfine\eta $ | 0.07% | 0.006% | 0.0007% |

$GCImedium\eta $ | 0.25% | 0.05% | 0.015% |

$GCIcoarse\eta $ | 0.87% | 0.46% | 0.31% |

$\alpha GV=6\u2009deg$ | $\alpha GV=10\u2009deg$ | $\alpha GV=13\u2009deg$ | |
---|---|---|---|

$N1,\u2009N2,\u2009N3\u2009[106]$ | 23.19, 9.09, 3.72 | 23.04, 9.42, 3.77 | 23.05, 9.42, 3.71 |

$r21,\u2009r32$ | 1.3556, 1.3569 | 1.3473, 1.3568 | 1.3474, 1.3641 |

$\eta 1,\u2009\eta 2,\u2009\eta 3$ | 92.478, 92.345, 91.891 | 95.133, 95.097, 94.784 | 94.070, 94.081, 93.838 |

$p$ | 4.025 | 7.134 | 10.215 |

$\eta ext21$ | 92.534 | 95.137 | 94.069 |

$GCIfine\eta $ | 0.07% | 0.006% | 0.0007% |

$GCImedium\eta $ | 0.25% | 0.05% | 0.015% |

$GCIcoarse\eta $ | 0.87% | 0.46% | 0.31% |

### 3.4 Validation of the Numerically Obtained Hill Chart.

The numerical method described in Sec. 3.3 is used to calculate the zone of the hill chart where the VSO curve is expected to be found. For that purpose, a nearly uniform triangulation grid of 33 operating points is selected from the experimental data in the range of $nED=0.15\u20130.20$ and $\alpha GV=4\u2009deg\u201314\u2009deg$. As it is shown in Refs. [39] and [40], this sampling offers a good balance between accuracy and number of operating points considered. Additionally, the point where the simulated BEP of the reference turbine is located was added in the computations to improve the accuracy of the $\eta BEP$ objective. The boundary conditions for the CFD model were iteratively matched to produce the exact head, or up to a centimeter of accuracy, to what was measured for each operating point. On Fig. 8 shown is a comparison between the experimental and simulated performance curves for the reference turbine, together with the 34 operating points of the hill chart. The efficiency curves are normalized with the respective BEP efficiencies found in the measurements and in the simulations. It can be concluded that in the zone of the VSO curve, i.e., for the range of $nED=0.17\u20130.19$, the numerical method is correctly capturing the efficiency trends of the reference turbine, with an error of less than $1%$ between the normalized efficiencies for the entire range of GV openings. In absolute terms, the simulations have overestimated the peak efficiency by $\Delta \eta =2.184%$, together with a nearly systematic underestimation of the discharge $Q$ for all operating points simulated (see vertical downshift of the simulated points in the hill chart on Fig. 8). Also, the simulation predicts the location of the BEP with a relative error of $\epsilon nED=nED,exp*\u2212nED,sym*nED,exp*\xb7100\u2248\u22128%$ for the speed factor and $\epsilon QED=QED,exp*\u2212QED,sym*QED,exp*\xb7100\u2009\u2248\u20094.6%$ for the discharge factor. The same numerical method, as well as the same operating points, is used to calculate hill charts of the new designs in the optimization.

### 3.5 Accuracy and Refinement of the Response Surface Models.

The BBD sampling was done in three levels $(\u22121,0,1)$ for the normalized ranges of the 15 free parameters between $\u22121$ and $1$, resulting in total of $m=421$ designs. To estimate the VSO objectives, the hill charts of each design sample was simulated, which required $421\u2009*\u200934=14,314$ calculations in total. Since the efficiency is only known in discrete locations, the VSO objectives are evaluated by a bi-harmonic spline interpolation between the 34 operating points. A spline was selected to ensure that the interpolated efficiency surface is smooth, has $C2$ continuity, and passes through the unstructured data points.

With 15 degrees-of-freedom for each RSM of the three characteristic efficiencies, the number of unknown regression coefficients $\beta $ in Eq. (2) is $p=136$ and were estimated using least squares minimization in matlab. To check the accuracy of the fitting, 30 new design samples (i.e., 1020 more calculations) were simulated and compared against predicted values from the models, which indicated an excellent accuracy and reliability of the RSMs. The additional 30 samples were than used for an additional refinement of the regression coefficients and the results from the goodness-of-fit tests for both runs RSM-1,2 are given in Table 5. Ideally, a perfect fit will result in $Radj2=1$ and $\sigma e=0$ and since both criteria for each objective are very close to the ideal values, the trend-capturing accuracy of the models is confirmed. Comparing the values for both runs, only minor improvement is achieved in the fitting accuracy so that further refinement of the models is assumed to be unnecessary.

$Radj2$− PL | $Radj2$ − BEP | $Radj2$ − HL | $\sigma e$ − PL | $\sigma e$ − BEP | $\sigma e$ − HL | |
---|---|---|---|---|---|---|

RSM-1 | 0.99822 | 0.98895 | 0.99518 | 0.00017 | 0.00029 | 0.00046 |

RSM-2 | 0.99801 | 0.98969 | 0.99541 | 0.00017 | 0.00028 | 0.00045 |

$Radj2$− PL | $Radj2$ − BEP | $Radj2$ − HL | $\sigma e$ − PL | $\sigma e$ − BEP | $\sigma e$ − HL | |
---|---|---|---|---|---|---|

RSM-1 | 0.99822 | 0.98895 | 0.99518 | 0.00017 | 0.00029 | 0.00046 |

RSM-2 | 0.99801 | 0.98969 | 0.99541 | 0.00017 | 0.00028 | 0.00045 |

The residuals between the calculated and predicted values are shown on Fig. 9. As can be seen, the prediction of the model throughout the entire sample range is frequently accurate to the first decimal of the turbine efficiency expressed in percentage. Few peaks can be spotted on the residual plots, especially visible in the modeling of the HL objective, which is found to be a result from several factors, such as numerical noise, occasional convergence problems, and hill chart interpolation error. The histogram representation with 20 bins is used to display the distribution of the residuals, which confirms that the error is clustered around the zero mean with a very small occurrence rate for samples with $|\epsilon |>1\xd7\u200910\u22121$. On Fig. 10, the tested samples are shown in two orthographic projections of the three-dimensional response space, together with 200,000 randomly distributed designs that are predicted with the trained RSMs. The characteristic efficiencies of the VSO objective are normalized with the respective efficiencies of the reference turbine that is found at $(1,1)$ in both charts, which splits the projections into four characteristic quadrants. Most of the tested samples are close to the reference design, meaning that the range of parameters is wide enough to include the reference and that the potential-flow-based design gives a good starting point for the optimization. The hollow curve represents the front of the performance improvement and is found by brute-force prediction with the models on more than 200 million random designs. The models predict that the peak efficiency of the reference design can be further improved and that it is not a Pareto optimal according to the defined VSO objectives.

## 4 Discussion on the Possibilities for an Optimal Design

From Fig. 10, if the design is found in the first quadrant on both charts at the same time, then it is superior to the reference turbine in all VSO aspects. Similarly, if the design is found in the third quadrant of both charts at the same time, then it is inferior to the reference turbine in all VSO aspects. Finally, if the design is found in at least one of the second and fourth quadrants of both graphs, it suggests that the design has been only partially improved. The best candidate for VSO will obviously be present in the first quadrant of both graphs and as far as possible from the reference point with coordinates $(1,1,1)$.

Due to this, the optimal design must be found as a tradeoff between the three VSO objectives, for which all design parameters will have different contribution. On Fig. 11, the models are used to assess the sensitivity around the center point for the three VSO objectives of the design space. Surprisingly, the hydraulic profile of the hub has very small effect on the improvement of the objectives and this holds true for any location on the response surfaces. As described in Sec. 3.2, one reason for this could be that initial blade angles are recalculated and corrected each time a change of the meridional view is done. Moderate sensitivity can be observed for the parameters controlling the TE shape in the meridional view and the leaning of the LE. This suggests that for low specific speed designs, the design parameters of low and moderate influence can be excluded from the hydraulic optimization of the runner and will have more influence on the overall weight and price, as well as the structural integrity and dynamic characteristics. The most influential parameters are found to be (1) the overall length of the blade, (2) the design values of the head and discharge, which actually repositions the location of the BEP of the design, (3) the leaning of the TE, (4) the correction factors of the outlet blade angles, and (5) the maximum thickness of the blade. However, nothing is excluded in the VSO optimization process presented in this paper and all 15 design parameters are optimized simultaneously to exploit their joint effects on the objectives. Despite the different sensitivity of the free parameters, some also have contradictory effects to the VSO objectives. For instance, reducing the overall length of the blade through the parameter $\phi SL1$ only will increase the efficiency for the $\eta BEP$ and $\eta HL$ objectives but, it will also reduce the efficiency at $\eta PL$, revealing the need for balancing between the design parameters.

To narrow down the analysis only on the designs that fulfill the BEP objective, on Fig. 12 plotted are the modeled responses of approximately 3 million random designs that have at least equal or higher efficiency than the reference turbine at the VSO BEP location. A grayscale colormap is used to color the designs according to the normalized BEP performance, where white color represents designs that have equal $\eta BEP$ objective with the reference design, while darker colors represent designs that are better. Due to the white background of the plot, some partially improved designs with $\eta BEP/\eta BEP,ref=1$ are not visible. Once again, Pareto optimal solutions can be found on the front seen in the first quadrant, suggesting that the BEP objective of the reference can be improved by 0.3% absolute efficiency at best. Two designs *D*1 and *D*2 that have equal efficiency with the reference at the BEP objective are selected on the opposite ends of the *VSO range* that represents the possible stretching of the hill chart. Compared to the performance of the reference at both PL and HL objectives, the D1 design underperforms by approximately 0.5%, while D2 outperforms by 0.2% in relative terms. The differences of their hill charts are shown on the same figure, visualizing the anticipated stretching/modification of the hill chart in the $QED$ direction, accompanied by a position shift of their respective BEPs that creates a difference between the speed numbers. For clarity, the corresponding efficiency levels are marked with $l1\u22128$, while the constant GV opening curves and the contour markers are not displayed. This proves that having only one objective is not enough because many different designs with improved off-design operation might coexist.

Designs which have $\eta BEP/\eta BEP,ref>1$ (darker colored on Fig. 12) are also found lying on the Pareto front next to D2, showing that better designs exists with an improvement in all three VSO objectives. A final selection of the optimal design according to all objectives was done based on linear search in a list of more than 200 million random designs that were evaluated with the trained RSMs. To confirm the model-predicted performance of the selected optimal design, a CFD calculation was done using the same method described in Sec. 3.4. The comparison of the CFD results between the VSO-optimal and the reference design is shown on Fig. 13, together with three characteristic views of both runners. As expected, the location of the peak efficiency for the optimized runner is found to be very close to the optimization requirements. The optimization leads to improvement of the efficiency in the entire hill chart area, suggesting that the optimal design also outperforms the reference at SSO as well. However, the results also show that a significant flattening of the VSO efficiency curve was not achieved, and this is mainly because the reference design was already a state of the art in its performance according to today's standards. This is also confirmed by the location of the reference design that is very close to the Pareto front predicted by the trained surrogate models, see Figs. 10 and 12. The narrow VSO range found between the D1 and D2 designs on Fig. 12 reveals an insensitive and somewhat predefined performance characteristic of the low specific speed turbine. Consistent and simultaneous improvement of both PL and HL objectives was achieved only with two parameters marked $T1r$ and $T2z$, and they have moderate-to-low contribution to the objectives, while the rest of parameters have contradictory effects (see Fig. 11). The optimized runner made a very small change on the guide vane opening characteristics and preserved the original discharge capacity of the turbine.

## 5 Conclusions

At off-design operating conditions of Francis turbines, the predominant losses are originating from the inlet and outlet of the runner itself and they can be balanced with appropriate adjustments of the rotational speed. Speed variation can be done to correct the speed factor $nED$ and operate along the hill chart profile found in the $\eta \u2212QED$ plane. To verify that the full potential is utilized from application of variable speed technology to an existing turbine, a case study is presented with an attempt for further improvement of the variable speed performance by means of brute-force optimization. A low specific speed Francis turbine was used as a reference, where a replacement runner was redesigned based on numerical simulations and response surface modeling.

The variable speed objective can be defined to be a tradeoff between three characteristic efficiencies, namely, one for the desired position of the peak efficiency and two objectives corresponding to discharge factor that is 20% higher and 50% lower than the optimal. These objectives are dependent on the design parameters and are found to have a nonlinear character that can be accurately approximated with a fully quadratic function. For numerical optimization purposes, a nonrestrictive parametrization of the runner can be achieved with 15 parameters using Bezier curves. The performance trends in the zone where the variable speed operation curve is expected to be found can be accurately predicted with simple steady-state simulations on a reduced domain and approximately 3.7 × 10^{6} elements.

Based on detailed examinations of the surrogate models, most of the design parameters are found to have contradictory effects on the objectives, meaning that the variable speed performance can be altered in a narrow range. These parameters contribute mainly to the position of the best efficiency point in the hill chart and the level of the peak efficiency. Only the parameters that control the shape and location of the trailing edge in the meridional view contributed to a simultaneous improvement of the high load and part load objectives. From this perspective, it is reasonable to conclude that the off-design performance is mainly governed by the design values of the rotational speed $n$ and the ratio of the inlet and outlet diameters of the runner. This conclusion is drawn without hard evidence because those parameters were not included in the optimization procedure. According to the results from numerical simulations, the selected optimal design has improved the efficiency by approximately 0.2% in the entire operating range, however better outcome may be expected for turbines of medium-to-higher specific speeds.

The proposed optimization method does not include the cavitation performance in the loop, which can be defined with an additional objective function, and instead of that the optimal runner was manually checked for the pressure distributions on the blades. Additionally, simulating the entire hill chart is computationally expensive, meaning that including more parameters in the optimization procedure will be impractical. On the other hand, the proposed method is general and, with small adjustments and variations, can be applied to improve the efficiency of synchronous speed Francis turbines as well.

## Acknowledgment

This research was supported in part with computational resources at the Norwegian University of Science and Technology (NTNU) provided by the “Vilje” cluster,^{2} where most of the numerical simulations were carried out.

## Funding Data

The work was funded by the Norwegian Research Center for Hydropower Technology—HydroCen (Project No. 257588) that is hosted by the Norwegian University of Science and Technology (NTNU).

## Nomenclature

- $b$ =
width $(mm)$

- $B$ =
Bernstein polynomials

- $cm$ =
meridional velocity $(m\u2009s\u22121)$

- $D$ =
diameter $(m)$

- $g$ =
gravitational acceleration $(m\u2009s\u22122)$

- $h$ =
minimum or maximum blade thickness $(mm)$

- $Hn$ =
net head $(m)$

- $i$ =
integer

- $j$ =
selection order

- $k$ =
number of nodal diameters; kinetic energy of the turbulence $(m2\u2009s\u22122)$; pipe constant including all lengths, areas, area ratios and friction coefficients

- $l$ =
iso-contour level

- $Ls$ =
length along a streamline $(m)$

- $m$ =
number of samples; harmonic number

- $n$ =
rotational speed $(rpm)$; number of design parameters; degree of Bernstein polynomials; degree of a Bezier curve; integer

- $N$ =
element count of a computational mesh

- $nED$ =
speed factor

- $nq$ =
specific speed

- $p$ =
number of regression coefficients; vector of points along a curve; pressure $(Pa)$; apparent order

- $P$ =
control point of the Hub polygon

- $P$ =
vector of control points

- $Q$ =
discharge $(m3\u2009s\u22121)$

- $QED$ =
discharge factor

- $r$ =
cylindrical coordinate; Refinement factor

- $R$ =
set of real numbers

- $Radj2$ =
adjusted

*R*^{2}coefficient of determination - $t$ =
parameter of a curve

- $T$ =
torque $(N\xb7m)$; control point of the trailing edge polygon

- $x$ =
Cartesian coordinate

- $x$ =
vector of design parameters

- $xi$ =
design parameter

- $y$ =
Cartesian coordinate

- $y\xaf$ =
mean of simulated responses

- $y+$ =
dimensionless wall distance

- $y(x)$ =
vector of simulated responses

- $y\u0302(x)$ =
vector of modeled responses

- $z$ =
cylindrical coordinate

- $Zg$ =
number of guide vanes

- $Zr$ =
number of runner blades

### Greek Symbols

- $\alpha GV$ =
guide vane opening angle $(deg)$

- $\beta $ =
regression coefficients; blade angle $(deg)$

- $\epsilon $ =
error; energy dissipation rate $(m2\u2009s\u22123)$

- $\eta h$ =
hydraulic efficiency

- $\rho $ =
density of water $(kg\u2009m\u22123)$

- $\sigma e$ =
root-mean-square error

- $\Phi $ =
angle of the leading edge lean $(deg)$

- $\phi SL$ =
wrapping angle of a quasi-streamline $(deg)$

- $\phi TE$ =
angle of the trailing edge lean $(deg)$

- $\psi $ =
stream function

- $\omega $ =
angular velocity $(rad\u2009s\u22121)$

### Subscripts and Superscripts

## Footnotes

## References

*Official Journal of the European Union, 5.6.2009,*pp. 16–62.