Optimization of Francis Turbines for Variable Speed Operation Using Surrogate Modeling Approach

(1)

Igor Iliev

¹

Waterpower Laboratory, Norwegian University of Science and Technology (NTNU), Alfred Getz’ Vei 4, Trondheim 7034, Norway e-mail: [email protected]

Erik Os Tengs

EDR & Medeso AS, Lysaker Torg 45, Lysaker 1366, Norway e-mail: [email protected]

Chirag Trivedi

Mem. ASME Rainpower Norge AS, Instituttveien 8, Kjeller 2027, Norway e-mail: [email protected]

Ole Gunnar Dahlhaug

Waterpower Laboratory, Norwegian University of Science and Technology (NTNU), Alfred Getz’ Vei 4, Trondheim 7034, Norway e-mail: [email protected]

Optimization of Francis

Turbines for Variable Speed Operation Using Surrogate Modeling Approach

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.[DOI: 10.1115/1.4047675]

1 Introduction

In Norway, as well as throughout most of Europe, the hydropower sector is facing an increasing demand for additional opera- tional 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, even- tually 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, with- out 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 factorn_EDis adjusted to ensure operation along a specific curve in then_EDQ_EDplane, which corresponds to the hill profile as projected on theg_hQ_EDplane. Other opera- tional 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.

1Corresponding author.

Contributed by the Fluids Engineering Division of ASME for publication in the JOURNAL OF FLUIDS ENGINEERING. Manuscript received January 6, 2020; final manuscript received June 2, 2020; published online August 5, 2020. Assoc. Editor:

Philipp Epple.

Downloaded from http://asmedigitalcollection.asme.org/fluidsengineering/article-pdf/142/10/101214/6554501/fe_142_10_101214.pdf by NTNU Universitets Biblioteket user on 06 August 2020

(2)

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 factornED¼nD2=pffiffiffiffiffiffiffiffiffiffiffiffigH_n

can be achieved by propor- tionally matching of the speednto the net headH_n. However, due to equipment and reservoir restrictions, not many sites exist with such operating conditions. For the latter, and perhaps more com- mon 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.

Any optimization problem requires a performance measure based on output parameters, also termed as anobjective function. The object/process that is subject to the optimization is normally described with a set of inputs, orfree 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 calledmulti-objective optimization.

For optimization of Francis turbines, the hydraulic efficiencyg_hat one or more operating conditions (Pareto optimal) is usually set to be maximized, given by the general mathematical definition

maximizeyð Þ;x xRⁿ;x¼ðx1;…;x_nÞsubject tox^l_ix_ix^u_i; i¼1;2;3;…;n

(1) 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 andu denote the lower and upper values, respec- tively, in the constrained range of theith 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 genera- tion 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 asg_BEP.

Turbine operation at off-design conditions can be classified into two categories, one being part load (PL) when the operating value ofQ_EDis lower than the optimal, the other being high load (HL) when the operating value ofQEDis 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

efficiency at HL and PL will be reduced by theDglosses with the subscriptsI;Sreferring to the incidence and swirl losses, respec- tively. At this point, it becomes obvious that the combination of Dg_IþDg_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 differentnEDvalue. In fact, for any Francis turbine, there exists an operating curve in thenEDQED

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 theperfect balancebetween the partial efficiency losses for a constantQED 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 thegQEDplane.

To achieve that, two more objectives are defined as: maximize the efficiency for two discharge factor levels Q^ðPLÞ_ED andQ^ðHLÞ_ED , irre- spective of the speed factorn_ED. For that purpose, the minimum and maximum values in the operating range ofQ_EDcan be used.

In the following sections, these two efficiencies are noted asg_PL andg_HL(see also Fig.1(b)), and for the case study presented in Sec.3, the two discharge factor levelsQ^ðPLÞ_ED andQ^ðHLÞ_ED are set to be at 50% lower and 20% higher than the optimalQ^ðBEPÞ_ED , respec- tively. 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 stretchingin 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.Perhaps one of the simplest ways to predict and model the behavior of any output parameter defined byn-input parameters is to use polynomial approximation models, a technique also known as response surface modeling (RSM). This is particularly suitable for modeling of the behavior of the VSO objectives defined in Sec.2.1, which are based on the free parameters used to describe the geometry of the turbine (see Sec.3.2). The idea is to use statistical tools and fit a multidimen- sional hypersurface of some degree to a data collected from a series of CFD calculations (samples) performed on the object/

process that needs to be optimized. Since building of the RSM is not an optimization method on its own, the search for the opti- mum configuration is done with a suitable optimizer [22,23], or by selection of a tradeoff candidate from a set of Pareto optimal solutions in a multi-objective case. For accurate prediction of the nonlinearities, the unknown function yðxÞ will be approximated with a second-order polynomial that gives the values of the pre- dictoryðxÞ^

^yð Þ ¼x b₀þXⁿ

i¼1

b_ixiþXⁿ

i¼1

Xⁿ

j1

b_ijxixjþXⁿ

i¼1

b_iix²_i (2)

wherebare the regression coefficients that need to be determined, x_i—linear terms, x²_i—pure quadratic terms, and x_ix_j—mixed/

interaction terms of the free parameters. Depending on the number of free parameters, i.e., degrees-of-freedom of the design, the minimum number of samples needed to estimate the unknown coefficients isp¼ ðnþ1Þðnþ2Þ=2. The approximation will produce a normally distributed random errorethat represents the difference between the CFD calculated response and the prediction by the fit- ted model

yð Þ ¼x ^yð Þ þx e;xRⁿ (3)

(3)

Because the system is overdetermined, for the given number of samplesmp, the regression coefficients are estimated by a least squares minimization. This, in matrix notation is written as

y¼Xbþe; b¼ðX^TXÞ¹X^Ty (4)

wherebis a column vector of sizep1 containing the coefficients,Xis ampdata matrix of themsamples on thepterms and is given by

X¼

1 x^{ð Þ}₁¹ x^{ð Þ}_n¹ x^{ð Þ}₁¹x^{ð Þ}₂¹ x^{ð Þ}_n1¹ x^{ð Þ}_n¹ x^{ð Þ}₁¹

2

x^{ð Þ}_n¹

2

⯗ ⯗ .. .

⯗ 1 x^{ð Þ}₁^m x^{ð Þ}n^m x^{ð Þ}₁^mx^{ð Þ}₂^m x^{ð Þ}_n1^m x^{ð Þ}n^m x^{ð Þ}₁^m

2

x^{ð Þ}n^m

2

2 66 66 4

3 77 77

5 (5)

The 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, theroot-mean-squared errorr_edefined as

r_e¼

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X^m

i¼1

e^{ð Þ}ⁱ ð Þ²=m s

; e^{ð Þ}ⁱ ¼

y^^{ð Þ}ⁱ y^{ð Þ}ⁱ y^{ð Þ}ⁱ

(6)

and theadjusted R²coefficient of determination, which is adjusted for the sample size and the number of regression coefficients in the model and defined with the expression

R²_adj¼1 P^m

i¼1

y^{ð Þ}ⁱ y^^{ð Þ}ⁱ

2

P^m

i¼1

y^{ð Þ}ⁱ y^{ð Þ}ⁱ

2m1

mp (7)

whereyis 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 choosempsamples 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 eitherthree-level full factorialorrandom sampling distribution, it is shown in Ref. [22] that the efficiency will be quite poor when the number of free parameters isn5.

Additionally, there is a lack ofrotatabilityin the three-level full factorial designs, where the variance of the predicted valuesyis 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¼3ⁿ, a three-level incomplete factorial design is proposed in Ref. [24]. Commonly known as theBox–Behnken designs (BBD), this sampling method has been specifically constructed for the use together with second-order RSMs and is able to producenearly rotatablesamples. A comparison between BBD and full factorial designs is shown on Fig. 2, where the number of free parameters isn¼3. The construction of the design BBD matrix requiresm¼2n nð 1Þ þ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

Fig. 1 Variable speed operation for maximized efficiency. (a) Predominant losses in the turbine resulting from a mismatch in the rotational speed and (b) the VSO curve identifies the most balanced state between the partial losses in the turbine.

(4)

hypercube, andquasi-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.2is 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.The model of the reference turbine is installed in the Waterpower laboratory at the Norwegian University of Science and Technology (NTNU) in Trondheim, Norway. The geometry is open, and it has been a part of extensive research through a series of projects and workshops regarding high-head Francis turbines [25,26]. The test rig can provide operating conditions and experimental precision in compliance with the international IEC standard for model acceptance tests [27]. The distributor of the Francis turbine has 14 stay vanes and 28 guide vanes, while the splitter-bladed runner is equipped with 15 full-length and 15 splitter blades. On Fig. 3 shown is a cross section through the entire passage of the turbine, as well as an isometric view of the hydraulic surfaces of the iso- lated turbine runner. The experimentally measured hill chart is shown on Fig.4, with a comparison on the respective efficiencies

between the VSO and SSO at optimal head. More details on the efficiency and pressure pulsations during VSO at different operating conditions can be found in Ref. [16]. The data were collected in an open-loop configuration, with more than 200 operating points measured in a wide operating range of the turbine model, and are presented using dimensionless factors, defined as [27]

n_ED¼ nD2

ffiffiffiffiffiffiffiffi gH_n

p ; Q_ED¼ Q D²₂ ffiffiffiffiffiffiffiffi

gH_n

p ; g¼ Tx

qgQH_n (8)

wherenðs¹Þis the rotational speed,qkg m³

is the density of the water at the operating temperature,gðm s²Þis the local gravitational acceleration, H_nð Þm is the net head of the turbine, Qðm³s¹Þ is the discharge through the turbine,TðNmÞis the Fig. 2 A sampling example of a three-level design on n53. (a) The black dots are the

Box–Behnken designs and are obtained by excluding the (b) gray dots from a full factorial design. The hollow marker is the center point (0;0;0) of the design.

Fig. 3 Cross section showing the geometry of the reference turbine and a detail of the splitter-bladed turbine runner

Fig. 4 Experimental results showing a comparison of the efficiencies between VSO and SSO at the optimal head for the reference turbine

(5)

torque of the runner,xðrad s¹Þis the angular velocity of the runner, andD2ð Þm is the outlet diameter of the runner. The characteristic dimensions and performance values are given in Table1. 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 610% 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].

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.To perform a VSO optimization with reasonable computational resources, the runner must be defined parametrically with minimum number of free parameters. In contrast to this, the resulting geometry must also be defined in a nonrestrictive fashion to ensure that an optimal design can be reached. Several studies suggest different approaches using parametric curves and surfaces, with degrees-of- freedom corresponding to the specifications of the optimization tasks [29–31]. Following the classical definition, the three- dimensional (3D) chamber surface of the runner blade can be described by combining onemeridional projectionof the runner with several conformal maps of thestream surfaces. To obtain the pressure and suction sides of the blade, a varyingblade thickness distribution is applied in both directions of the surface normal vectors. In this paper, all these steps are parametrized using quadratic/cubicBezier curves, generally defined as [32]

pð Þ ¼t Xⁿ

i¼0

Bⁿ_ið ÞPt _i; Bⁿ_ið Þ ¼t n i

tⁱð1tÞⁿⁱ (9)

wherepð Þt is the position vector of a point lying on the curve,Pis the vector comprised of the control points,Bⁿ_ið Þt are theBernstein polynomials,

n i

are thebinomial coefficients,t½0;1is a parameter of the curve,i½0;nis an integer number of the sum,nis 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.).

In the meridional view, the hub profile is defined with a cubic Bezier curve, while the profile of the trailing edge (TE) is defined with a quadratic Bezier curve. The allowed directions of move- ment for each of the control points (hollow circles of the poly- gons) are shown on Fig. 5. Fixed and reused from the existing runner are the profile of the shroud, the leading edge, and the starting point of the hub curve that was necessary to preserve the inlet widthb1. The position of the midpoint in the TE polygon, marked with a hollow square, is defined to move with respect to the position of the two endpoints of the TE curve and it is not a free parameter of the design. The blade is represented by four spanwise sections—streamlines, where for each streamline calculated are initial blade angles b^ðjÞ_1;2;j¼1–4 at the inlet—1 and the outlet—2 using design values forH_n,Q, andnfrom the existing turbine. For simplification, approximated meridional velocities jcmjare calculated for axisymmetric potential flow with neglected influence from the runner blades, given by the mathematical definition

jcmj ¼1 r

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

@w

@z

2

þ @w

@r

2

s

@²w

@r²1 r@w

@r þ@²w

@z² ¼0 8>

>>

><

>>

:

(10)

wherewðr;zÞis thestream functionin 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 ofH_nandQto correct the blade angles globally. Hence, for fixed rotational speedn, the headH_nand dischargeQbecome 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

b^ðjÞ_2corr¼b^ðjÞ_2iniþDb^ðjÞ; j¼1–3 (11) whereDb^ðjÞare the free parameters constrained in a specific range for each streamline.

The two-dimensional planar streamlines are expanded to the 3D space of the runner by applying an inverse design rule

du¼ 1

rtanbdLs (12)

whereuis thewrapping angleof the blade at any location and dL_s 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 u_SLðjÞ¼Ð2

1 1

rtanbdL_s. The overall length of the blade is controlled by allowing u_SL1 to be a free parameter within a constrained range, while the remainingu_SLðjÞforj¼2–4 are controlled by a leaning angleu_TEwith a predefined TE shape as seen in thetop viewof the runner on Fig.5. Additionally, a leaning of the LE is enabled through the free parameterUthat rotates the streamlines in the circumferential direction, with its effect best seen on the side view of the runner on Fig.5.

A symmetric thickness distribution is defined using two con- catenated cubic Bezier curves that can be controlled by seven parameters [34]. The number of parameters is further reduced to only one free parameterh1by setting a fixed value for the minimum thickness h2 and using modified Gottingen-410 airfoil.

Based on the required maximum thickness of the blade h1, the Table 1 Dimensions and specifications of the model at BEP

Parameter Value

Net head,Hn 12 mð Þ

Discharge,Q 199 l=sð Þ

Rotational speed,n 323 rpmð Þ

Speed factor,n_ED 0:173ð Þ–

Discharge factor,Q_ED 0:151ð Þ–

Guide vane,a_GV 9:8 degð Þ

Peak efficiency,g 92:6 %ð Þ

Inlet diameter,D1 0:620 mð Þ

Outlet diameter,D2 0:349 mð Þ

Blade angle,b_1ðshroudÞ 82:4 degð Þ

Blade angle,b_2ðshroudÞ 18:2 degð Þ

Inlet width,b1 59:6 mmð Þ

Max blade thickness,h1 6:9 mmð Þ

Min blade thickness,h2 1:8 mmð Þ

(6)

airfoil is scaled and elongated to fit the length of the streamlines together with a TE modification for reduction of the vortex shed- ding amplitudes. On Fig. 6 shown is a random runner design,

defined by 15 parameters in total. To keep the number of parameters as low as possible, the replacement runner is defined to have only full-length blades. Even though the same parametrization can be applied to include splitter blades as well, however, this will double the number of free parameters and increase the cost of optimization to an unacceptable level. Also, the number of runner blades is not included in the list of free parameters. Instead, the number of blades is selected in the typical range ofZ_r¼16–19 and to avoid excitation ofmodal shapeswith less than 5nodal diameters by the fundamentaland second harmonic of theRSI according to Ref. [35]

mZ_g6k¼nZ_r (13)

Fig. 5 Display of the parametric definition of the runner. The meridional view is controlled with six parameters marked with the arrows; the global length of the blade is controlled with the wrapping angle of the hub streamlineu_SL1; the leaning of the TE is controlled withu_TE; the leaning of the leading edge is controlled with the angleU.

Fig. 6 A random runner design after a blade thickness distribution is applied. The blade thickness is controlled with h1

only, which scales a predefined Gottingen-410 to fit each section of the chamber surface accordingly.

Table 2 Significant RSI excitation patterns

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

Table 3 Constraints of the design

Parameter Range

rof pointP1ðmmÞ ½170;200

rof pointP2ðmmÞ ½50;80

zof pointP2ðmmÞ ½80;40

zof pointP3ðmmÞ ½250;180

rof pointT1ðmmÞ ½50;80

zof pointT2ðmmÞ ½171;151

u_SL1ðradÞ ½1:1;1:6

UðdegÞ ½1;3

u_TEðdegÞ ½0;30

Db^ð1Þ₂ ðdegÞ ½0;12

Hn=H_nð Þ% ½0:95;1:05

Q=Q ð Þ% ½0:75;1

h₁ðmmÞ ½15;25

(7)

wherem¼1;2;3… is the harmonic number,Zg¼28 is the number of guide vanes,n¼1;2;3… is an arbitrary integer, andkis 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 beZr¼17. The free parameters used to define the geometry of the runner are summarized in Table3, 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 thez¼0 plane passing through the mid- dle of the channel with widthb1.

3.3 Description of the Numerical Model and Estimation of the Numerical Error.Numerical prediction of the performance is done usingANSYS CFXby solving the steady-state Reynolds aver- aged Navier–Stokes equations coupled with the standard ke 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., Dp_SCþSV¼kQ². 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 constantkwas 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 dischargeQto 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 indexwas used [37,38]. Since the turbine efficiency is the focus, and due to the lack of conservation property in the definition of the grid Fig. 7 A 3D view of the coarse mesh of the entire guide vanes ring, runner and draft tube, con-

sisting of approximately 3.7 M structured hexahedral elements in total. On the right, a detail is given of the mesh near the guide vane—runner interface for the three guide vane openings.

Table 4 Discretization error for the reference case

a_GV¼6 deg a_GV¼10 deg a_GV¼13 deg

N1;N2;N3½10⁶ 23.19, 9.09, 3.72 23.04, 9.42, 3.77 23.05, 9.42, 3.71

r21;r32 1.3556, 1.3569 1.3473, 1.3568 1.3474, 1.3641

g₁;g₂;g₃ 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

g²¹_ext 92.534 95.137 94.069

GCI^g_fine 0.07% 0.006% 0.0007%

GCI^g_medium 0.25% 0.05% 0.015%

GCI^g_coarse 0.87% 0.46% 0.31%

(8)

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 y^þ200 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 50% of the elements are in the runner, 30% are in the draft tube, and the remaining20% 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 factorr1:35 in all three directions. The meshes were simulated for three operating points, namely: (1)a_GV¼6 deg, (2)a_GV¼10 deg, and (3)a_GV¼13 deg, 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 methodp. 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 10⁵CPU hours.

3.4 Validation of the Numerically Obtained Hill Chart.

The numerical method described in Sec.3.3is 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 ofnED¼0:15–0:20 anda_GV¼4 deg–14 deg. 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 theg_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 Fig. 8 Comparison between experimental and numerical results for the reference turbine. A triangulated grid of 33 operating points, plus the simulated BEP, were used.

Table 5 Indicators for the RSM goodness-of-fit of the VSO objectives

R²_adjPL R²_adjBEP R²_adjHL rePL reBEP reHL

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

(9)

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 ofnED¼0:17–0: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 byDg¼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 ofe_n_ED¼ⁿ^ED;exp_nⁿ^ED;sym

ED;exp 100 8% for the speed factor ande_Q_ED ¼^Q^ED;exp_Q^Q^ED;sym

ED;exp 100

4: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ð1;0;1Þ for the normalized ranges of the 15 free parameters between1 and 1, resulting in total ofm¼421 designs. To estimate the VSO objectives, the hill charts of each design sample was simulated, which required 421 34¼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 inter- polated efficiency surface is smooth, hasC²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 b 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 Table5.

Ideally, a perfect fit will result inR²_adj¼1 andr_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. Compar- ing 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.

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 Fig. 9 Magnitude and distribution of the residuals for the three RSM models after the first refinement of the regression coefficients. 20 bins are used for the histogram representation

Fig. 10 Two orthographic projections of the VSO objectives. With 200,000 random samples, the model already identifies a saturation pattern that forms a front of potentially optimal solutions.

(10)

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 jej>110¹. 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

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

Fig. 12 Two characteristic designsD1 andD2 predicted to have the same efficiency at the BEP objective but different efficiencies at the PL and HL objectives. The values of the PL and HL objectives for both designs are on the opposite side of the tradeoff range. The comparison of the predicted hill charts of the corresponding designs reveal the anticipated “hill-chart stretching,” plus a position shift of the respective BEPs, i.e., a change of the specific speeds.

Fig. 11 One-dimensional contours of the 15-dimensional RSM showing the parameter sensitivity at the center point. Different parameters have different, and sometimes contradictory, effect for the three VSO objectives/responses.

(11)

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 parameteru_SL1only will increase the efficiency for theg_BEPand g_HLobjectives but, it will also reduce the efficiency atg_PL, reveal- ing the need for balancing between the design parameters.

To narrow down the analysis only on the designs that fulfill the BEP objective, on Fig.12plotted 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 equalg_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 withg_BEP=g_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 D1 andD2 that have equal efficiency with the reference at the BEP objective are selected on the opposite ends of theVSO rangethat 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 theQ_EDdirection, accompanied by a position shift of their respective BEPs that cre- ates a difference between the speed numbers. For clarity, the corresponding efficiency levels are marked with l₁₈, while the

constant GV opening curves and the contour markers are not dis- played. This proves that having only one objective is not enough because many different designs with improved off-design operation might coexist.

Designs which have g_BEP=g_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 and12.

The narrow VSO range found between the D1 and D2 designs on Fig.12reveals 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 T₁^r andT₂^z, 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.

Fig. 13 Comparison of the top, bottom, and side views between the reference and the VSO optimal design. Numerically calculated hill charts of both designs, as well as their respective efficiency curves at VSO, are plotted on top of each other.

(12)

ments of the rotational speed. Speed variation can be done to correct the speed factorn_EDand operate along the hill chart profile found in thegQEDplane. To verify that the full potential is uti- lized 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 rede- signed 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.710⁶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 speednand the ratio of the inlet and outlet diameters of the runner. This conclusion is drawn with- out 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,² 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)

Nomenclature b¼widthðmmÞ

B¼Bernstein polynomials c_m¼meridional velocityðm s¹Þ

D¼diameterðmÞ

g¼gravitational accelerationðm s²Þ

h¼minimum or maximum blade thicknessðmmÞ H_n¼net headðmÞ

i¼integer j¼selection order

k¼number of nodal diameters; kinetic energy of the turbu- lenceðm²s²Þ; pipe constant including all lengths, areas, area ratios and friction coefficients

l¼iso-contour level

L_s¼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 n_ED¼speed factor

n_q¼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ðm³s¹Þ Q_ED¼discharge factor

r¼cylindrical coordinate; Refinement factor R¼set of real numbers

R²_adj¼adjustedR²coefficient of determination t¼parameter of a curve

T¼torqueðNmÞ; control point of the trailing edge polygon x¼Cartesian coordinate

x¼vector of design parameters x_i¼design parameter

y¼Cartesian coordinate

y¼mean of simulated responses y^þ¼dimensionless wall distance yð Þ ¼x vector of simulated responses

^

yð Þ ¼x vector of modeled responses z¼cylindrical coordinate Z_g¼number of guide vanes Z_r¼number of runner blades Greek Symbols

aGV¼guide vane opening angleðdegÞ

b¼regression coefficients; blade angleðdegÞ e¼error; energy dissipation rateðm²s³Þ g_h¼hydraulic efficiency

q¼density of waterðkg m³Þ r_e¼root-mean-square error

U¼angle of the leading edge leanðdegÞ u_SL¼wrapping angle of a quasi-streamlineðdegÞ u_TE¼angle of the trailing edge leanðdegÞ

w¼stream function

x¼angular velocityðrad s¹Þ

Subscripts and Superscripts i¼counter

I¼incidence losses l¼lower value S¼swirl losses u¼Upper value

2https://www.hpc.ntnu.no/display/hpc/Vilje