1 Introduction

The draft tube is an important part of a hydroelectric power installation, as its function is to decelerate, in an orderly fashion, the flow leaving the runner in such a way as to convert as much as possible of the remaining kinetic energy back into static pressure, thus increasing the effective head the runner is subjected to. However, there are complications which affect its performance including viscous turbulent effects (such as wall separations) generated by the draft tube geometry. To complicate matters further, for most of its operation, the hydraulic turbine operates at off-design conditions, which the design of the draft tube must accommodate to ensure a consistent performance. The design of the draft tube therefore needs to be thought out carefully to achieve the best possible compromise of hydraulic efficiency and construction costs. Traditionally, the design of a draft tube has been tempered through experimental observations and semi-empirical formulae derived from established geometries (notably: Gubin 1973). To explore potential new designs, Computational Fluid Dynamics (CFD) has proved to be a powerful tool for the engineer, allowing for comprehensive analysis of complex flowfields where experimental work provides limited insight (Abbas and Kumar 2017). CFD becomes especially appealing when combined with an optimisation method which may significantly reduce the number of evaluations during the design cycle.

1.1 Hölleforsen–Kaplan draft tube geometry

The draft tube geometry considered in this work is a 1:11 scaled model of the Hölleforsen–Kaplan draft tube, constructed in 1949. This design has served extensively as a benchmark test case for both experimental and numerical studies in the literature—largely through the European Research Community On Flow, Turbulence And Combustion (ERCOFTAC) Turbine-99 Workshop series (Gebart et al. 2000; Engström et al. 2001; Cervantes et al. 2005). This design was selected because of the wealth of literature studying it, and particularly the experimental and CFD data available for validation and design comparison. Additionally, the underlying engineering design principles and fluid dynamics for this case are well understood, giving us a solid baseline for comparison and understanding any new designs. A schematic of the draft tube geometry is shown in Fig. 1.

Fig. 1
figure 1

Schematic of 1:11 scaled model of the Hölleforsen–Kaplan sharp-heeled draft tube Andersson (2008). All dimensions are in mm. In this paper, we focus on optimising the inflow cone (red), elbow (green) and diffuser (blue) sections of the geometry. (Color figure online)

1.2 Optimisation of the Hölleforsen–Kaplan draft tube—previous studies

General optimisation strategies involve modification of design variables with a view to improve appropriate performance goals (the objective functions). Often the objective functions are non-linear and multi-modal, and hence polynomial time algorithms for solving such problems may not be available. In some cases, Evolutionary Algorithms (EAs) have been used to explore the global optimum (Zingg et al. 2008). However, the high computational costs associated with evaluating a large number of objective functions prevent applications of EAs to many practical engineering design problems. In order to cut the prohibitive costs, a surrogate model may be used to reduce the number of required objective function evaluations (see Queipo et al. 2005). However the regions of design parameter space of greatest interest are unlikely to be known in advance and the a-priori surrogate model may not be particularly accurate in these regions. Bayesian optimisation is a process in which the computer iteratively ‘learns’ the best surrogate model for the design space as part of the strategy to find the optimum solution. The results of this article are one element of an overall research project applying Bayesian optimisation to engineering design problems—in this particular contribution, the performance of a hydraulic turbine draft tube.

Automated shape optimisation of the of the Hölleforsen–Kaplan draft tube has been reported by a few authors. Marjavaara and Lundström carried out investigations into the redesign the outer-heel using circular Marjavaara and Lundström (2006) and elliptic Marjavaara (2006) inserts. The Response Surface Method was the chosen surrogate modelling strategy, while the CFD was performed using the commercial code CFX4.4. Overall, the authors report that only a small improvement was made to the pressure recovery for the optimum heel design in their investigation. To the best of our knowledge, no study has been conducted for the shape optimisation for the other components of this draft tube design, which is one of the novelties of our work presented here. For alternate elbow-draft tube designs, parameter-based optimisation has largely been conducted through fine-tuning dimensions of fixed geometrical features (e.g. so-called “Swedish” parameterisation (McNabb et al. 2014)) or by the dimensions of cross-sections along the draft tube (e.g. Eisinger and Ruprecht 2002). At the same time, these studies have largely been carried out using a single-objective optimisation, with the few multi-objective studies showing little analysis of justification for the relations between the objectives. To the best of the authors’ knowledge, an appropriate multi-objective investigation of the draft tube design is yet to be presented in the literature, let alone a parameter-based optimisation without any conceptual barriers to explore potential new designs. Clearly there are benefits to being able to include human understanding into the design process, not least in being able to substantially reduce the design space to search—and in the past, when only a few CFD runs could be undertaken, this was of paramount importance—and in providing additional practical constraints such as might be related to retrofitting improvements to existing installations. However, doing this does also introduce human bias into the optimisation, and now with comparatively cheaper CFD there is also an argument in favour of allowing a more unrestricted search through design space, with the possibility of throwing up designs which have not previously been considered. Since the objective here is to demonstrate the Bayesian Optimisation process on a realistic engineering problem, rather than necessarily to identify new solutions for the hydroelectric industry, this will be the focus of the present work. A Python-based framework combining multi-objective Bayesian optimisation with CFD was developed for this purpose, and has already been demonstrated (see Daniels et al. 2018 or Daniels et al. 2019b); the current paper represents an attempt to apply this methodology and framework to a more complex, real-world case.

1.3 Paper overview

Following on from the above review we have identified the following topics to be addressed in this paper:

  1. 1.

    To propose a new methodology for parametrisation of the draft tube components for each evaluation;

  2. 2.

    To perform an efficient multi-objective optimisation using Bayesian optimisation methods for pressure recovery and wall-frictional losses in all three draft tube components simultaneously;

  3. 3.

    To evaluate the appropriateness of the draft tube performance quantities for multi-objective shape optimisation.

The structure of this paper reflects the stages of work undertaken towards achieving the above goals. López (2006) provides a 4-step methodology for the optimization process in turbomachinery, consisting of (1) Draft tube parameterization, (2) CFD methodology and validation, (3) Optimization algorithm and (4) Objective function. We have discussed i and ii in an earlier paper (Daniels et al. 2020), so Sect. 2 reviews the methodology for evaluating the draft tube performance which has been presented in more detail there, including a discussion of the objective functions (Sect. 2.1) and the CFD methodology (Sect. 2.2). Section 3 describes the methodology for multi-objective optimisation, as well as the parameterisation for deforming the geometry for each evaluation—thus addressing the 1st topic of this work. The details of possible objective functions and their impact on the flow were also discussed in the earlier paper. The results and discussion of the optimisation are discussed in Sect. 4, focusing on the 2nd and 3rd topic of this paper. Finally, in Sect. 5, the observations, and suggestions for future work are summarised.

2 Performance evaluation of the draft tube

In this section we set out the theoretical framework for evaluating the performance of an arbitrary draft tube in terms of pressure recovery and friction factors. Then in Sect. 2.2 we present an overview of the CFD methodology being used to determine the flow properties, and hence the performance of the system.

2.1 Draft tube performance measures

The efficiency of a draft tube is primarily measured by the pressure recovery factor:

$$\begin{aligned} C_{p} = \frac{1}{\frac{1}{2}\rho \left( \frac{Q}{A_{in}} \right) ^{2}}\left[ \frac{1}{A_{out}} \int _{A_{out}}{\overline{p}}_{out} dA_{out} - \frac{1}{A_{in}} \int _{A_{in}}{\overline{p}}_{in} dA_{in}\right] , \end{aligned}$$
(1)

where \(\rho\) is the density of the fluid, \({\overline{p}}\) is the static pressure, Q is the volumetric discharge, and A denotes the cross-sectional area for the inflow (in) and outflow (out) boundaries respectively. Additionally, the total energy-loss factor, \(\zeta\), is also employed to measure the draft tube performance:

$$\begin{aligned} \zeta _{1} = \frac{1}{\frac{1}{2}\rho \left( \frac{Q}{A_{in}} \right) ^{2}}\left[ \frac{1}{A_{in}}\int _{in}P_{t,in} dA_{in} - \frac{1}{A_{out}}\int _{out}P_{t,out} dA_{out}\right] , \end{aligned}$$
(2)

where \(P_{t}\) is the total pressure, i.e., \(P_{t} = {\overline{p}} + 0.5\rho (U_{i}^{2})\). Alternatively, an approximation for \(\zeta _{1}\) may also be derived by an algebraic manipulation of Eqs. 1, and 2:

$$\begin{aligned} \zeta _{1} \approx 1 - C_{p} - \frac{A_{in}^{2}}{A_{out}^{2}}, \end{aligned}$$
(3)

which, to a certain degree, implies \(C_{p}\) and \(\zeta _{1}\) have an anti-correlated relation under the following conditions:

  1. 1.

    The flowfield can be described by a fully potential flow, providing an adequate account of the flow pattern (streamlines) for the interior of the draft tube under consideration.

  2. 2.

    The velocity distribution across the draft tube, and deviations from the potential flow approximation, are reasonably captured by a single multiplicative constant to the kinetic energy term in Bernoulli’s equation.

With regard to the first premise, it should be noted that the potential flow approximation within the draft tube is adequate to represent the flow pattern as the inertial time scale is generally much smaller than the viscous time scale for transient (rapidly accelerating) flows. Regarding the second premise, near-uniform velocity distributions along a diffuser has previously been shown to improve the pressure recovery, with a high-correlated linear relationship between the two performance measures (Daniels et al. 2019a). However, with the presence of swirl flow, turbulence and the incurring unsteady flow phenomena, some deviation from the intrinsic relation between \(C_{p}\) and \(\zeta _{1}\) is expected. Furthermore, it should be noted that this relation varies considerably if the optimisation approach allows for cross-sectional areas of the inflow and outflow to vary, in which case the last term of Eq. 3, \(A_{in}^{2}/A_{out}^{2}\), becomes dominant. In the present work, the cross-sectional areas of the inflow and outflow remain constant throughout the optimisation procedure. Therefore, an anti-correlation between \(C_{p}\) and \(\zeta _{1}\) is expected in the results. Finally, in the literature, the energy loss of the draft tube has been expressed in an alternate form:

$$\begin{aligned} \zeta _{2} = \frac{1}{\frac{1}{2}\rho \left( \frac{Q}{A_{in}} \right) ^{2} U_{i}\cdot n}\left[ \frac{1}{A_{in}}\int _{in}P_{t,in} U_{i}\cdot n dA_{in} + \frac{1}{A_{out}}\int _{out}P_{t,out}U_{i}\cdot n dA_{out}\right] , \end{aligned}$$
(4)

where \(\cdot n\) indicates the component normal to the corresponding boundary. An algebraic relation between \(\zeta _{2}\) and \(C_{p}\) is not easily derived. However, another form of the pressure recovery coefficient is known to display a similar relation with \(\zeta _{2}\) (Andersson et al. 2004).

Despite the apparent connection between the various forms of \(C_{p}\) and \(\zeta\), a number of publications have carried out a multi-objective design optimisation of the draft tube using \(C_{p}\) and \(\zeta\) as objective functions. A mentioned above, Marjavaara and Lundström (2006) performed a multi-objective optimisation focusing on the shape of the outer-heel, minimising \(\zeta _{2}\) as the second objective. Regardless of the additional \(U_{i}\cdot n\) components, with increasing heel radius, a trend between \(C_{p}\) and \(\zeta _{2}\) is evident in their results, but is not addressed in their discussion; in an earlier publication (Marjavaara and Lundström 2003), the same design optimisation was carried out using a single objective consisting of a weighted-sum of the two performance quantities. Demirel et al. (2017) performed a multi-objective shape optimisation of another elbow-type draft tube, minimising a slight variation of \(\zeta _{1}\) (head-loss) as the second objective. The authors note that a close relation between pressure recovery and head-loss was evident, and subsequently analysed their design based on the pressure recovery. Over two papers, Galván et al. (2013a, b) optimised the inflow velocity profiles of the draft tube, minimising \(\zeta _{1}\) as a single objective. Other performance quantities (e.g., \(C_{p}\)) were calculated from the evaluated designs. Again, a connection between \(C_{p}\) and \(\zeta _{1}\) can be observed in their work, but is not described in their discussion. To justify the application of a multi-objective optimisation procedure, this work will adopt a different expression for frictional losses as the second design objective, and will demonstrate the close relation between \(\zeta _{1}\) and \(\zeta _{2}\) with \(C_{p}\) by deriving the values of \(\zeta\) from the evaluated designs.

In general, frictional losses are a significant economic concern for engineering designs. In order to recover the residual kinetic energy, that would otherwise be converted into static head, the draft tube design must minimise frictional losses while simultaneously reducing the fluid’s velocity. One way this can be achieved is by gradually enlarging the cross-sectional areas along the draft tube. The draft tube shapes generated during an optimisation run may include rapid cross-sectional expansions, leading to large flow separation and other detrimental flow phenomena, which may impair the draft tube’s efficiency. On the other hand, if the expansion along the draft tube is too gentle, the fluid will be exposed to an excessive wall surface area, leading to large frictional losses, flow separation and a more expensive construction. As reported in the literature, frictional losses, and other detrimental flow phenomena in the draft tube (e.g. pressure surging and cavitation) can be mitigated using actively controlled methods, such as the use of air bubble jets at the vicinity of the inflow wall (Chirkov et al. 2017), or by passive methods, such as implementation and modification of ribs (Štefan et al. 2012). To the best of the authors’ knowledge, consideration for the effects of the draft tube shape on the wall-frictional losses in the design optimisation has yet to be considered in the literature, and will therefore be conducted in this work.

In this multi-objective shape optimisation, the energy losses along the draft tube are quantified by the frictional forces along the draft tube walls. In one form, the pressure losses at the walls can be quantified in terms of the nondimensional wall pressure, \(C_{p,wall}\) (see Daniels et al. 2020), as shown in Galván et al. (2013a) or Galván et al. (2013b). Alternatively, for this work, the frictional losses are quantified by the area-weighted average of skin friction across the draft tube walls:

$$\begin{aligned} \overline{C_{f}} = \frac{1}{A_{wall}}\int _{wall}\frac{(\nu + \nu _{t}) \big |\frac{\partial U}{\partial y_{1}}\big |}{\frac{1}{2} \left( \frac{Q}{A_{in}} \right) ^{2}} dA_{wall}, \end{aligned}$$
(5)

where \(A_{wall}\) is the surface area of the draft tube, \(y_{1}\) corresponds to the height of the first wall-normal cell to the draft tube wall, and U is the velocity magnitude. \(\nu\) and \(\nu _{t}\) are the kinematic and turbulent viscosity of the fluid. It is expected that this objective is equivalent to maximising the wall pressure as a quality function (Galván et al. 2013a, b), as \(C_{p,wall}\) and \(C_{f}\) are inherently related (Cervantes and Gustavsson 2006). Thus, in this work, the second objective for optimisation is to minimise the average skin friction along the draft tube walls, which, while expected to display an anti-correlated relation, has a greater independence to \(C_{p}\) than \(\zeta\).

2.2 CFD evaluation methodology

The reliability of the optimization process depends on the underlying CFD results; it is thus critically important to make sure that the CFD methodology is reliable. This has been the subject of an earlier paper (Daniels et al. 2020) which assesses the numerical methodology for evaluating the performance of a draft tube design, the results of which are summarised here. The CFD simulations in this work were performed using the open-source code OpenFOAM-4.x. The numerical modelling follows the specified setup by the organisers of the 2nd Turbine-99 Workshop (Engström et al. 2001). As part of this setup, the velocity components from the Kaplan turbine were linearly interpolated onto the inflow boundary of the draft tube. For the optimisation run in this work, the operating conditions were kept constant at the ‘T(n)’ mode (Cervantes et al. 2005) detailed in Table 1. The flow through each draft tube design was simulated using the Reynolds-Averaged Navier–Stokes equations, using the standard \(k-\epsilon\) model for the turbulent viscosity. The choice of the \(k-\epsilon\) model was primarily based on existing studies in the literature (e.g. Abbas and Kumar 2017; Cervantes et al. 2005) although in our previous paper (Daniels et al. 2020) we also investigated the use of the \(k-\omega -SST-SAS\) model.

Table 1 Kaplan turbine operating modes ‘T(n)’ and ‘R(n)’ Cervantes et al. (2005)

The sharp-heeled Hölleforsen–Kaplan draft tube was used as the base geometry. Steady-state and time-averaged transient simulations were carried out to assess the performance measures of the draft tube (see Sect. 2.1) with near-identical results as shown in our earlier paper (Daniels et al. 2020). Sensitivity of the results to the inlet and outlet conditions was also investigated in this earlier work. The inlet conditions, which are also the exit conditions from the turbine, may have a significant impact on the flow (and in particular flow separation) in the draft tube and hence the actual optimal solution. However since the aim of the current paper is to investigate the optimisation methodology rather than any precise installation, this is possibly less important than it might be. Similarly the outflow was extended by \(2\;\mathrm{m}\) as recommended in Engström et al. (2001) to avoid backflow and ensure convergence of the results. The simulations were also extensively validated against experimental and numerical data by various authors.

In the earlier paper, a number of proposed draft tube design features from the literature and the present optimisation study were also evaluated using this CFD methodology. Design variants investigated included changes to the inflow cone section (overall diameter, convex and concave designs), and the elbow section (curved-, expanded- and chamfered-heel designs). The objective of this work was to support and explore the CFD methodology across the range of possible designs that might be generated during the optimisation process. A predominately uniform hexahedral CFD grid was generated within each draft tube geometry using cfMesh (Juretić 2017). The appealing aspect of this software is its ability to generate a quality grid within an enclosed geometry regardless of its cumbersome shape. The geometries evaluated during the optimisation run in this article were assessed using the user-defined parameters for cfMesh outlined in Table 2. The Grid Convergence Index (Celik et al. 2008) (Richardson extrapolation) was applied to a number of draft tube designs with various geometrical features. The uncertainty of the pressure recovery relating the grid topology/resolution and extreme design features considered in this work are shown in Table 3.

Table 2 User-defined parameters used in cfMesh and resulting total number of cells for the base design
Table 3 Uncertainty of \(C_{p}\) (Eq. 1) determined through Grid Convergence Index (Celik et al. 2008) performed in Daniels et al. (2020) on the base design and with varying geometrical features

3 Multi-objective shape optimisation

The primary goals of this work are to simultaneously maximise the pressure recovery factor and minimise the average skin friction coefficient along the walls of the draft tube. Any particular geometry can be represented as an n-dimensional decision vector \({\mathbf {d}}\in {\mathcal {D}}\), where \({\mathcal {D}}\subseteq \mathbb {R}^n\) consists of all feasible shapes. From this the objectives may be written as:

$$\begin{aligned} \min _{{\mathbf {d}}\in {\mathcal {D}}} f_{1}({\mathbf {d}})&= -C_{p}, \end{aligned}$$
(6)
$$\begin{aligned} \min _{{\mathbf {d}}\in {\mathcal {D}}} f_{2}({\mathbf {d}})&= \overline{C_{f}}. \end{aligned}$$
(7)

These two objectives potentially conflict and are to a certain degree anti-correlated (see Sect. 2.1): a design that improves one objective results in an almost proportional decrease in performance in the other objective. This is a common situation in multi-objective optimisation; the result is that there is no unique solution to this problem, but a set of optimal designs referred to as the Pareto set which represents the optimal trade-off between the different objectives. This trade-off between objectives is described using the concept of dominance (Coello Coello et al. 2007). A shape \({\mathbf {d}}\) is said to dominate another shape \({\mathbf {d}}'\), denoted as \({\mathbf {d}}\prec {\mathbf {d}}'\), if,

$$\begin{aligned}&f_1({\mathbf {d}})< f_1({\mathbf {d}}') \text { and } f_2({\mathbf {d}}) \le f_2({\mathbf {d}}')\nonumber \\ \text {or }&f_1({\mathbf {d}}) \le f_1({\mathbf {d}}') \text { and } f_2({\mathbf {d}}) < f_2({\mathbf {d}}'). \end{aligned}$$
(8)

Using this concept, the Pareto set is the set of non-dominated solutions :

$$\begin{aligned} {\mathcal {P}}= \{{\mathbf {d}}~ | ~{\mathbf {d}}' \nprec {\mathbf {d}}~ \forall {\mathbf {d}},{\mathbf {d}}' \in {\mathcal {D}}\wedge {\mathbf {d}}\ne {\mathbf {d}}'\}. \end{aligned}$$
(9)

A solution in the Pareto set cannot be improved in any objective without worsening in the others. The image of the Pareto set \({\mathcal {P}}\) in the objective space is known as the Pareto front \({\mathcal {F}}\). It may be impossible to identify the exact Pareto set within a reasonable time period, even if the objective functions are cheap to evaluate. Thus the goal of any practical optimisation algorithm is to find an approximation of the Pareto set (which we denote as \({\mathcal {P}}^*\subseteq {\mathcal {D}}\)) and associated estimated Pareto front (denoted as \({\mathcal {F}}^*\)).

The optimisation methodology that we use to find this Pareto front is that of Bayesian Optimisation (Shahriari et al. 2016). In Bayesian Optimisation, an approximation to the system cost function surface is iteratively learned, using an acquisition function to identify the best location at each iteration to evaluate the (expensive) CFD simulation. To carry out this process, we must first parameterise the geometry, details of this being discussed in Sect. 3.1. Section 3.2 presents the theory of Multi-objective Bayesian Optimisation, whilst we set out the details of our implementation in Sect. 3.3.

3.1 Draft tube geometry representation

The focus of this paper is to manipulate the inflow cone, outer-heel, and diffuser sections in order to locate a set of geometries \({\mathcal {P}}^*\) that closely approximates the optimal trade-off between pressure recovery factor and wall-frictional losses. Generally, an effective optimisation method takes a decision vector and alters parts of the vector intelligently to locate designs that potentially improve the current approximation \({\mathcal {P}}^*\). However, to evaluate the performance of a design, a Computer-Aided Design (CAD) for CFD is required. Therefore, a useful representation for this problem should be easily convertible from vectorial form to a CAD model. It should be noted that the computational cost of an optimisation problem is highly influenced by the number of design variables. Hence, there is much to gain in finding as minimal a number of relevant design variables as possible, as well as design space constraints. In addition, properly chosen design variables in CFD applications will alter the geometry and CFD grid in a proper and sound manner. Taking this into consideration, Sect. 3.1.1 will outline the design variables for the optimisation process used in the present work.

3.1.1 Vectorial representation

In this section, the construction of the design vector \({\mathbf {d}}= [d_1, \dots , d_n]^\top \in {\mathcal {D}}\) is discussed.

Inflow cone section

Fundamentally, the inflow cone is a surface in three dimensions. Considering a cross-section in parallel to the xz-plane, the perimeter of this surface may be concisely represented with a radius r. This is considered as the first element of a decision vector. Thus, \(d_1 = r\) where \({\mathbf {d}}= [d_1, \dots , d_n]^\top\) is an arbitrary shape. The side of the inflow cone was represented by a single Catmull–Rom spline (Catmull and Rom 1974), possessing \(C^{1}\) parametric continuity. The spline implementation is indicated in Fig. 2a.

Fig. 2
figure 2

A demonstration of the inflow cone design generation considered for optimisation; a a schematic of the inflow cone with the bounds and location for the control point, and Catmull–Rom spline implementation; b the base design. All dimensions are normalised by the inflow cone diameter

Elbow section

For the outer-heel, a single control point \({\mathbf {h}}= [h_x, h_z]^\top \in \mathbb {R}^2\) and two end points of the Catmull–Rom spline is used within a two-dimensional bounding box defined on the cross-section parallel to the xz-center-plane to represent a curve on the heel. The implementation of this spline is outlined in Daniels et al. (2020). \({\mathbf {h}}\) is included in a decision vector \({\mathbf {d}}\): \(d_2 = h_x\) and \(d_3 = h_z\). Figure 3b shows the resulting Catmull–Rom spline implementation, control point, and bounding box for the heel design. This representation is capable of reproducing the base design in the optimisatio run.

Fig. 3
figure 3

A demonstration of the heel design generation considered for optimisation; a base heel construction using proposed heel representation; b schematic of the Catmull–Rom spline implementation, control point, and bounding box; c a demonstration of the deformed heel using the spline formation in (b). All dimensions are normalised by the inflow cone diameter

Diffuser section

Catmull–Clark subdivision curves (Stam 1998) were utilised to represent the profile of the upper and lower diffuser walls. The three-dimensional surfaces were generated by extruding the two-dimensional curves along the remaining axis. Since the curves are completely defined by a finite set of control points, optimisation comprises the process of manipulating the positions of the control points to generate the optimal shape. The design space for the k control points can be represented as a vector of control point coordinates, where the ith control point \({\mathbf {c}}^i = [c_x^i, c_z^i]^\top \in \mathbb {R}^{2}\). Horizontal (x) and vertical (z) components of the control points are combined into a decision vector thus: \(d_{2(i-1)+4} = c_x^i\) and \(d_{2(i-1)+5} = c_z^i\), where \(i\in [1,k]\). In addition, the optimisation process is constrained by the requirements that all control points must reside within a predefined volume of space, and that the resulting curves are simple (non-intersecting) polygons.

To demonstrate the Catmull–Clark subdivision curve methodology, Fig. 4 (top-right) shows the use of one control point to generate the top and bottom walls of the diffuser section. The control points (red triangles) are randomly generated within predefined bounding-boxes (circle markers, dotted lines). The bounding-box for the lower wall was defined to ensure that the lowest point does not fall below the base of the heel, while the bounding-box for the top surface is constricted to be no higher that the outflow cross-section. To ensure a smooth transition between the diffuser and neighbouring sections, fixed points are prescribed at either end of the top/bottom subdivision curves (red/blue squares). For the optimisation run, \(k = 3\) was chosen for the subdivision curves; this allows for the base geometry to be considered in the design space. From a practical perspective, usually only a few iterations of subdivision are needed to produce a visually smooth curve; the iteration limit was set to 5 here.

Fig. 4
figure 4

Schematic of the Catmull–Clark subdivision curve setup, and a randomly generated subdivision curve using one control point for the top and bottom walls. All dimensions are normalised by the inflow cone diameter

Given the representations of the three geometries outlined in this section, optimisation strategies such as EAs, for instance NSGA-II Deb et al. (2002), would be able to identify a good approximate Pareto set \({\mathcal {P}}^*\). However, EAs generally require thousands of function evaluations, which renders them impractical with computationally expensive evaluations. An alternative is to use multi-objective Bayesian optimisation, which has proved to be an effective approach with limited number of function evaluations (Rahat et al. 2017); the methodology for this will be outlined in the proceeding section.

3.2 Multi-objective Bayesian optimisation

Multi-objective Bayesian optimisation (MBO) is a surrogate model-based global search strategy that sequentially samples the design space at likely locations of shapes which may improve the current approximation of the Pareto set. One approach towards MBO is to reduce (or scalarise) the multi-objective problem into a single objective problem (Rahat et al. 2017). This enables the use of standard single objective Bayesian optimisation method in which a data-driven stochastic model (usually a Gaussian process model) is constructed, and then the prediction from the model is used to locate promising solutions (Shahriari et al. 2016).

To ensure that sequential sampling improves the current approximation of the Pareto set \({\mathcal {P}}^*\), the formulation of scalarisation function may be relative to the quality of the \({\mathcal {P}}^*\) (Rahat et al. 2017). In this paper, the hypervolume (Zitzler 1999)—a set based indicator that measures the volume of objective space covered between a non-dominated set and a predefined reference vector—is used as a quality measure for the \({\mathcal {P}}^*\). It is an exceptional indicator as maximising it is equivalent to locating the optimal Pareto set (Fleischer 2003). The details of the MBO approach with this scalarisation technique are given in Rahat et al. (2017), but for completeness a brief description of the method is provided below and summarised in Algorithm 1.

figure a

MBO starts (line 1) with a space filling design (e.g. Latin Hypercube Sampling (McKay et al. 2000)) of the feasible design space. The initial set of M feasible shapes \(D = \{{\mathbf {d}}^1, \dots , {\mathbf {d}}^M\}\) are then expensively evaluated with appropriate CFD simulations, and this provides a data set \(S = \{({\mathbf {d}}^t, f_1({\mathbf {d}}^t), f_2({\mathbf {d}}^t))\}_{t=1}^M\) (lines 3 –6). On each subsequent step (lines 7–18), an approximation of the Pareto set may be made as follows:

$$\begin{aligned} {\mathcal {P}}^*= \texttt {nondom}(S) = \{{\mathbf {d}}~|~ {\mathbf {d}}' \nprec {\mathbf {d}}, {\mathbf {d}}' \in S \wedge {\mathbf {d}}\ne {\mathbf {d}}'\}. \end{aligned}$$
(10)

The hypervolume improvement scalarisation \(g({\mathbf {d}}, {\mathcal {P}}^*)\), relative to the approximate Pareto set, for all the evaluated designs is then calculated (line 10), which allows a Gaussian process surrogate model to be trained capturing the relationship between design space and the hypervolume improvement (line 12).

The predictive density from the \(\mathcal {GP}\) model for a shape \({\mathbf {d}}\) is denoted \(p({\hat{g}}({\mathbf {d}})| S)\). This enables the calculation of a utility function to be obtained by querying the surrogate model for any desired design. Here the expected improvement in hypervolume improvement (with respect to the best value observed so far) is used as the utility. Given the best evaluated shape \({\mathbf {d}}^* = {{\,\mathrm{arg\,max}\,}}_{{\mathbf {d}}\in D}g({\mathbf {d}}, {\mathcal {P}}^*)\) (line 13), the expected improvement of an arbitrary feasible shape \({\mathbf {d}}'\) is defined as:

$$\begin{aligned} \alpha ({\mathbf {d}}', {\mathbf {d}}^*) = \int _{-\infty }^{\infty } \max ({\hat{g}}({\mathbf {d}}') - g({\mathbf {d}}^*, {\mathcal {P}}^*), 0) ~p({\hat{g}}({\mathbf {d}})| S)~ d{\hat{g}}({\mathbf {d}}). \end{aligned}$$
(11)

As the predictive distribution is Gaussian, this integral can be calculated in closed form. Thus, selecting the next shape to evaluate is decided by finding the optimum of the infill criterion, namely the solution of the following sub-problem: \({\mathbf {d}}^{t+1} = {{\,\mathrm{arg\,max}\,}}_{{\mathbf {d}}\in {\mathcal {D}}} \alpha ({\mathbf {d}}, {\mathbf {d}}^*)\) (line 14). Bi-POP-CMA-ES was used to locate a good approximation of the optimum in the sub-problem (Hansen 2009). The candidate solution \({\mathbf {d}}^{t+1}\) is expensively evaluated (15) and the data sets D and S augmented with the newly evaluated shape (lines 16 and 17). This process continues until the computational budget is exhausted, at which point the current approximation of the Pareto set \({\mathcal {P}}^*\) is returned. A single design from \({\mathcal {P}}^*\) may then be selected by the decision maker based on the desired performance and a knowledge of the optimal set of possible trade-offs between the objectives.

3.3 Computational details

The MBO method as proposed in Rahat et al. (2017) was used to generate the results in this paper.Footnote 1 The optimiser requires a few parameters to be set a priori, which are discussed here.

At the outset it is essential to set the number of initial evaluations so that an initial surrogate model may be constructed. A conventional rule of thumb for this is to set \(M = 11n - 1\), where n is the number of dimensions in the decision vector \({\mathbf {d}}\) (Jones et al. 1998). Here \(n=11\) and so the number of initial samples was set to 120. The boundary of the decision space is specified by lower and upper limits corresponding to each dimension, and Latin Hypercube Sampling used to generate an initial set of decision vectors \(D = \{{\mathbf {d}}^1, \dots , {\mathbf {d}}^{M}\}\). Additional constraints were imposed to ensure that no control points resided beyond the predefined bounding box as shown in Figs. 3 and 4. Any Latin Hypercube samples violating these constraints were rejected. As a consequence, in one instance, only 22 out of the 120 initial samples were deemed feasible (thus \(M=22\)). Note that the precise number of initial feasible solutions may vary as repeated Latin Hypercube samples do not have to generate the same pattern of cases. The other settings used in the MBO framework are as follows: the maximum number of surrogate evaluations to maximise the expected improvement using Bi-POP-CMA-ES was set to 12000n, the reference point used to compute the hypervolume \({\mathbf {l}}= [l_1, l_2]^\top\) was specified dynamically as \(l_i = \max _{{\mathbf {d}}\in D} f_i({\mathbf {d}}) + 0.1 (\max _{{\mathbf {d}}\in D} f_i({\mathbf {d}}) - \min _{{\mathbf {d}}\in D} f_i({\mathbf {d}}))\). Additionally, while locating \({\mathbf {d}}^{t+1}\) using Bi-POP-CMA-ES, each candidate solution is checked for feasibility (which is computationally cheap to determine). An infeasible solution \({\mathbf {d}}'\) is deemed to have zero expected improvement, i.e. \(\alpha ({\mathbf {d}}', {\mathbf {d}}^*) = 0\). This discounts infeasible solutions, and restricts the CFD calculations to only evaluate feasible solutions. For each CFD calculation, a steady-state simulation was performed using the numerical settings specified in Sect. 2.2 (see Daniels et al. (2020)). Each CFD simulation was deemed to be converged when all residual tolerances had reduced to less than \(10^{-6}\), a cutoff value which is deemed sufficient for most engineering problems (Michael and Torsten 2000). The iterative convergence for the evaluation of \(C_{p}\) was also included in this criterion. Simulations were carried out in parallel on one node of a dedicated CFD computer, comprising 2x Intel Haswell E5-2640v3 2.6GHz processors (16 cores total). The wall-clock time for individual evaluations was approximately 2 hours. Due to such exorbitant computational time for each evaluation, the overall number of CFD evaluations were limited to 200 evaluations. This means a complete optimisation cycle required approximately 2.5 weeks of total computation time.

Fig. 5
figure 5

Performance of the multi-objective Bayesian optimiser over ten independent runs depicted in terms of summary attainment surfaces. The blue, red and green lines show the attainment of \(20\%, 50\%\) and \(80\%\) of the fronts. The small variations between 20th and 80th percentile surfaces indicate the desirable repeatability and convergence properties of the algorithm. The green star and the associated dashed green lines depict the performance of the base geometry

Fig. 6
figure 6

Hypervolume progression for the last 150 evaluations of the MBO over 10 repetitions. The narrowing width in addition to the little progress in the median hypervolume achievement in the last 25 evaluations indicates that further evaluations are unlikely to deliver major improvements in the solution

The performance of the multi-objective Bayesian methodology can be described with a summary-attainment-surface (see either Fonseca and Fleming 1996 or Knowles 2005). To achieve this, the optimisation run was repeated 10 times, the resulting summary-attainment-surface of the estimated Pareto fronts is shown in Fig. 5. The small variations between the 20% and 80% percentile surfaces indicate the desirable repeatability and convergence properties of the optimiser. In addition, the hypervolume progression in Fig. 6 indicate that major improvements may not be achieved with further evaluations. It can be seen that with only 200 evaluations, all solutions in \({\mathcal {P}}^*\)s (found by the MBO in independent runs) outperform the base shape in \(\overline{C_f}\), while most solutions in \({\mathcal {P}}^*\)s also provide an improvement in \(C_p\). The MBO locates a range of solutions in each run trading off \(C_p\) and \(\overline{C_f}\), and on average the solutions in \({\mathcal {P}}^*\)s provide a change of approximately \(6\%\) in \(C_p\) and \(29\%\) in \(\overline{C_f}\) from the best to the worst performing solution in respective criteria across all runs. Thus, the proceeding section will focus its analysis on an arbitrarily chosen optimisation run.

4 Optimisation results and discussion

To develop the analysis of the results, an instance from the 10 independent optimisation runs must be chosen. Here, the run with the hypervolume closest to the median hypervolume over the repetitions (see Fig. 6) is selected. The close-up of the estimated Pareto front (\({\mathcal {F}}^*\)) from this optimisation run is shown in Fig. 7, and the overall objective space is shown in Fig. 8. A few draft tube designs evaluated during the optimisation run are indicated in these diagrams for discussion. The data points in these diagrams are grey-scaled according to the evaluation number in the optimisation run. The base design is indicated by a green star.

Fig. 7
figure 7

The estimated Pareto front (\({\mathcal {F}}^*\), red squares) of the two cost functions (Eqs. 1, 5) for the median dominated hypervolume optimisation run, annotated with associated draft tube designs. Design (vii), which has the lowest value of \(\zeta\), is not shown here as it is shown in Fig. 12 and discussed in more detail in the text. (Color figure online)

Fig. 8
figure 8

Objective space of the two cost functions (Eqs. 1, 5) for the median dominated hypervolume optimisation run, annotated with the estimated Pareto front (\({\mathcal {F}}^*;\) red squares) and associated draft tube designs. (Color figure online)

The objective space (Fig. 8) shows a wide range of conventional and unusual designs for the draft tube. When considering the first objective, \(\min -C_{p}\), it can be seen that the designs that perform poorly relative to the base design (i.e. right-hand side of the base geometry) are also unconventional in design. It can be seen that these designs largely possess sharp increases of curvature in the inflow cone and diffuser sections, resulting in aggressive changes in cross-sectional area along the geometry. The incurred flow separation has a detrimental impact on the pressure recovery, as discussed in our earlier paper (Daniels et al. 2020). However, it is not necessarily the case that sharp increases of curvature along the diffuser has a negative effect on efficiency. For example, it can be seen in Fig. 7 that the diffuser shape for ‘v’ forms a vortex-cavity along the upper-wall, which is known to reduce wall-frictional losses and improve pressure recovery (Mariotti et al. 2013), compare this to ‘vi’, whose design incurs a large separation region along the diffuser and poorer pressure recovery. It should be noted that under the current representation, increasing the number of control points will expand the design space, and potentially discover better performing diffuser designs with increasing undulations. However, this simultaneously increases the difficulty to manufacture the design. Thus, the engineer has to consider the trade-off between potentially improving the draft tube performance against computational cost with increasing number of control points, as well as the manufacturability of the resulting design. Finally, for the more conventional designs (those whose features are currently seen on existing draft tubes), to the left-hand-side of the base geometry, can be seen closer to (and on) the estimated Pareto front (\({\mathcal {F}}^*\))—indicated by red squares.

Figure 7 shows a close-up of \({\mathcal {F}}^*\) from the optimisation run. As is the nature of multi-objective optimisation, the estimated optimal trade-off front leaves the final choice of optimal design to the engineer. In the present work, a few designs selected along \({\mathcal {F}}^*\) are labelled ‘\(i-vii\)’ for further examination. It can be seen that these designs show three distinct trends:

  1. 1.

    The radius of the inflow cone is consistently wider than that of the base design for all designs along \({\mathcal {F}}^*\);

  2. 2.

    The designs of the heel along \({\mathcal {F}}^*\) (from top-left to bottom-right) progresses from a chamfered \(\rightarrow\) curved \(\rightarrow\) sharp \(\rightarrow\) expanded shape, i.e. the cross-sectional area in the heel increases along \({\mathcal {F}}^*\);

  3. 3.

    The curvature and divergence (cross-sectional area) of the diffuser section increases along \({\mathcal {F}}^*\) (top-left to bottom-right).

For a deeper insight into the performance of these designs, a series of sample planes are placed along each geometry as indicated in Fig. 9 (top). The pressure recovery was calculated on these planes using Eq. 1, where out is synonymous with the position of the plane (e.g., \(p_{out} = p_{A}\) at position A). The heel design impacts the pressure recovery along the diffuser section, as can be seen in Fig. 9 for ‘i’ and ‘vii’, which possess an inflected and chamfered (slightly larger) heel design respectively. It can be seen that the pressure recovery for the inflected heel design (‘i’) increases the pressure recovery along the diffuser compared to the equivalent draft tube design with a chamfered heel (‘vii’). Therefore reducing the cross-sectional area of the heel increases the pressure recovery along the draft tube.

Fig. 9
figure 9

Pressure recovery evaluated along the draft tube cross-sections for various designs identified along the Pareto front \({\mathcal {F}}^*\) (Fig. 7)

The optimisation run has also identified a region in which the estimated optimum inflow cone radius was just slightly larger than the base geometry (See Fig. 10), with the optimum radius (those along \({\mathcal {F}}^*\)) was approximately 0.02 m wider. The equivalent divergence angle of the inflow cone has increased from \(2\theta _{cone} = 21.96^{\circ }\) to \(37.02^{\circ }\); this new angle is within the design guidelines Gubin (1973). Figure 10 shows the relation of pressure recovery to inflow cone radius around the hub for the evaluations of the optimisation run. It can be seen that for the smaller inflow cone radius (concave), the pressure recovery reduces as the axial velocity increases through this region. As the radius of the inflow cone increases (convex), the pressure recovery develops an indifference to the inflow cone radius above 0.375 m.

Fig. 10
figure 10

Pressure recovery versus inflow cone radius from the evaluated designs during the optimisation run. The x-axis (radius) is the normalised by the diameter of the inflow cone, \(D_{0} = 0.5\) m (Andersson et al. 2004)

Overall, from the optimisation run, the maximum pressure recovery attained was 0.992 (‘i’, Fig. 7), an increase of 3.7% on the base design. The minimum average skin friction attained, while maintaining an improvement on pressure recovery (‘v’), decreased by 22.3% to the base design. It can be seen in Fig. 8 that to a certain degree, an anti-correlated relation exists between pressure recovery and wall-fictional losses (Sect. 2.1). An alternative quantification of frictional losses, and typical second objective: \(\min \zeta _{i}\) (Eqs. 2, 4), can be quantified from the evaluated designs.

Fig. 11
figure 11

Objective space showing \(-C_{p}\) vs \(\zeta _{i}\) for designs from the median dominated hypervolume optimisation run. In a the blue line shows the relation Eq. 3, whilst the red line is a best-fit line through the calculated results

Table 4 Correlation coefficients for \(C_{p}\) and \(\overline{C_{f}}\) with \(\zeta _{i}\)

Figure 11 shows the objective space relation between \(\zeta\) and \(C_{p}\) for designs evaluated from the optimisation runs, together with the expression Eq. 3 for \(\zeta _1\). Clearly, in both diagrams, it can be seen that there is a distinct anti-correlation between \(C_{p}\) and \(\zeta _{i}\). The corresponding Pearson correlation coefficients between \(C_{p}\) and \(\zeta _{i}\) have been quantified in Table 4, confirming the highly anti-correlated relationship discussed in Sect. 2.1. The correlation coefficients between \(\overline{C_{f}}\) and \(\zeta _{i}\) have also been quantified in Table 4, demonstrating the inherent relation between wall-frictional losses and \(\zeta\) to describe the frictional losses within the draft tube. To recall that a linear relationship (Eq. 3) between \(C_{p}\) and \(\zeta _{1}\) can be derived by algebraic manipulation of Eqs. 1 and 2—the blue-dotted line in Fig. 11a represents this expression. A curve-fitted (red-dotted) line shows a similar gradient and axis interception to Eq. 2. However, due to the assumptions outlined in Sect. 2.1, discrepancies are seen between the two lines. Nevertheless, it can be seen that Eq. 2 makes a good approximation to the relation between \(C_{p}\) and \(\zeta _{1}\), provided that the cross-sectional area of the outflow is constant throughout the optimisation run. The analysis of these cost functions demonstrate that the use of \(C_{p}\) and \(\zeta\) should be used as a single objective optimisation case, and engineers should consider alternative performance quantities for multi-objective optimisation.

Fig. 12
figure 12

Estimated design of the draft tube for minimal \(\zeta _{i}\) and \(-C_{p}\). This design is labelled ‘vii’ in Fig. 7

The design which has the lowest value of \(\zeta _{i}\) is shown in Fig. 12. This design is located on the Pareto front in Fig. 7 (‘vii’). It can be seen in Fig. 12 that the design has a slightly wider inflow cone than the base design, chamfered heel, and a straight-tapered diffuser section. This diffuser design is typically seen in modern elbow-type draft tubes, such as the U9 (Mulu et al. 2012). The use of a slightly convex inflow cone and chamfered heel, on the other hand, is a new design characteristic determined in this work. However, a reduction of heel size is another design characteristic seen in modern draft tubes, which invariably possess a curved heel.

4.1 Efficiency curves of the designs

It is important to analyse the performance of the identified designs across a range of turbine operating conditions. The flow entering the draft tube has both a tangential and axial velocity component, i.e. swirl. The magnitude of the swirl is, for fixed blade runners (Francis), usually a function of the discharge (Q). A Kaplan turbine is doubly-regulated (runner blade-guide vanes) which can adjust the level of swirl at each operating point. In this work, a constant associated blade angle is used for all loads, producing conditions similar to those of Francis turbines. Altering the discharge through the draft tube alters the swirl characteristics through the draft tube. For example, at part-load conditions (e.g. T(n)), the swirl is rather large, which potentially gives rise to vortex breakdown below the runner hub (Mulu et al. 2012); In this case, local negative axial velocity is created, impeding the overall performance of the draft tube (Čarija et al. 2006).

The optimisation run in this work used a constant turbine operating condition T(n) (see Table 1). However, it is not necessarily the case that the new designs are operating at their peak efficiency at this condition. From a range of unit discharges at the inflow, an optimal working condition of the draft tube can be identified from the resulting efficiency curve. The unit discharge is defined:

$$\begin{aligned} Q_{11} = \frac{4Q}{\pi D_{0}^{2} \sqrt{2gH}}, \end{aligned}$$
(12)

where \(D_{0}\) is the diameter of the inflow cone, g is the gravitational acceleration, and H is the pressure head of the turbine. \(D_{0}\) and H were taken as 0.5 m and 4.5 m respectively based on the experimental setup by Andersson (Andersson et al. 2004).

For the 1st Turbine-99 Workshop (Gebart et al. 2000), two operating points on the efficiency curve were chosen for the benchmark study for the base geometry. These points were originally denoted T(r) and R(r)—‘T’ for top, and ‘R’ for the right leg on the turbine efficiency curve. However, during measurements, a mechanical breakdown occurred on the benchmark experimental setup (Andersson et al. 2004) resulting in a new set of operating conditions. The two new points were specified—T(n) and R(n) (Table 1), both of which had a lower discharge than their original counterparts.

Fig. 13
figure 13

Pressure recovery versus unit discharge curves for various heel configurations

The efficiency curves of the draft tube with various heel designs are shown in Fig. 13. Based on their wall-pressure results, Lövgren et al. (2007) determined that the positions of the T(n) and R(n) conditions are located on the left-leg and top of the efficiency curve. In the present work, the operating conditions T(n) and R(n) are located on the left-leg and top of the propeller curve respectively, in agreement with their deduction. Based on this, the measured peak pressure recovery for the base design is 0.975 at 0.46 unit-discharge—with an increase of approximately 1.94% from the T(n) operating condition. For the curved heel, proposed by Dahlbäck (1996), it can be seen that the pressure recovery is consistently larger than the base design along the left-leg of the base efficiency curve, and is characterised by a shallower curve. The peak efficiency is located at a lower unit-discharge to the base design, with an increase of 0.959% at peak performance. For the chamfered heel design, it can be seen that this has a shallower curve than the curved heel design, i.e., less sensitive to the volumetric discharge. The peak performance for this is located at a lower discharge than the curved and sharp heel designs; with a potential increase of 1% to the base design at the R(n) operating condition. Finally, for the expanded heel (vortex-generator design (Duarte et al. 2016)), the peak performance is located at a greater discharge than the base design, with an efficiency increase of approximately 1% at the same discharge as the base design. Overall, the chamfered heel has shown to be more robust than the other designs for off-design working conditions. The designs other than the sharp-heel show a marked improvement in pressure recovery at the left-leg of the efficiency curve, while no improvement can be seen around the best operating mode of the base design.

Fig. 14
figure 14

Pressure recovery versus unit discharge curves for various designs identified on the estimated Pareto front \({\mathcal {F}}^*\) (Fig. 7)

Figure 14 shows the efficiency curves for the designs identified along \({\mathcal {F}}^*\) (see Fig. 7). For designs ‘i’ (concave heel) and ‘vii’ (chamfered heel), located at the top-left of the Pareto front, it can be seen that reducing the size of the heel section increases the peak pressure recovery and moves its location to a lower discharge. On the other hand, the designs located towards the bottom-right of \({\mathcal {F}}^*\), comparing ‘v’ and ‘vi’, or ‘ii’ and ‘vii’, suggests that the high expansion of diffuser sections produces a steeper curve—reducing the robustness of the draft tube. Flatter curves will operate across a wider range of operating conditions, which may be valuable as Kaplan turbines operate over a wide range of discharges; this is especially appealing for lower-head sites such as lowland rivers. Consistently, the location of the peak efficiency of the draft tube for all designs along \({\mathcal {F}}^*\) are located at a lower discharge than the base design. It may however, be the case that the optimum design of draft tube is obtained at the same discharge as the base design. To overcome this, a possible design objective to consider in a future optimisation run is to maximise the pressure recovery at multiple discharges.

5 Conclusions and future work

5.1 Conclusions

Automated shape optimisation of a sharp-heeled Hölleforsen–Kaplan draft tube was performed based on the combined use of Computational Fluid Dynamics and a multi-objective Bayesian methodology for the first time. Geometric variables that specify the diameter of the inflow cone, shape of the outer-heel and diffuser were chosen. The total number of parameters used was 17, requiring approximately 200 evaluations to establish the Pareto front. Analysis of the performance quantities showed the typically used energy-loss factor and pressure recovery were highly correlated in cases of constant outflow cross-sections, and therefore unsuitable for use in a multi-objective optimisation. Thus, the pressure recovery factor and the wall-frictional losses were selected as objective functions. The objective space included a wide range of design features, some conventional, some more unusual. Sharp changes in curvature resulting in flow separation have a detrimental effect on the pressure recovery but not necessarily on the overall efficiency, eg through reduction in wall friction. The design with the lowest \(\zeta\) values was found to have a slightly convex inflow cone, chamfered heel, and straight-tapered diffuser. For this design, the pressure recovery factor was increased by 3.7% and wall-frictional losses were reduced by 22.3% compared to the base design. Finally, designs identified along the estimated Pareto front were analysed over a range of turbine discharges. It was found that reducing the size of the heel increased the draft tube performance over a wide range of turbine discharges. Overall, this work demonstrates the potential to optimise new or existing parts of the hydraulic turbine using a Bayesian methodology.

5.2 Challenges and future work

The comment “all models are wrong [but some are useful]” is often attributed to the statistician Box (1976). In this context it is worth examining the limitations of the study presented here. We have chosen a CFD methodology which has previously been shown to be robust at estimating the performance quantities used in the optimisation (as discussed in the earlier paper (Daniels et al. 2020)), however it is a steady-state approximation to a probably transient flow. Hence the effects of pressure fluctuation such as vibration and cavitation are not examined in the design optimisation. Bayesian optimisation itself is suitable for relatively small numbers of design parameters, in fact the 17 considered here is probably towards the limit of what can be achieved, and increasing the number of parameters will significantly increase the number of iterations needed. Additionally, this is a sequential optimisation, whilst in engineering contexts, parallelisation of the algorithm to run on multiprocessor HPC systems would certainly be desirable.

This work naturally leads to the following topics of investigation involving optimisation:

  • Additional design variables including the inflow velocity profiles, with an inference of optimising the Kaplan turbine design;

  • Optimising the runner hub geometry—providing a greater potential for pressure recovery and geometric flexibility than the inflow cone;

  • Optimising the draft tube while considering the robustness. To achieve this, evaluating the design at a range of discharges should be conducted, with the aim of maximising the pressure recovery for all cases. However, this would require a CFD simulation to be run multiple times for each evaluation—effectively increasing the computational cost.