Optimization of hydraulic turbine diffuser

Hydraulic turbine diffuser recovers pressure energy from residual kinetic energy on turbine runner outlet. Efficiency of this process is especially important for high specific speed turbines, where almost 50% of available head is utilized within diffuser. Magnitude of the coefficient of pressure recovery can be significantly influenced by designing its proper shape. Present paper focuses on mathematical shape optimization method coupled with CFD. First method is based on direct search Nelder-Mead algorithm, while the second method employs adjoint solver and morphing. Results obtained with both methods are discussed and their advantages/disadvantages summarized.


Introduction
Swirl turbine [1] is suited for low heads and relatively large discharges.Almost a half of the available head should be exploited in the draft tube by recovery of kinetic energy into pressure energy.Velocity field in the draft tube is characterized by swirl component on its inlet.This imposes stress on the draft tube performance [2].Liquid, which flows through draft tube (diffuser), is affected by hydraulic losses.These losses are undesirable and there is an effort to minimize them.The term minimization means process of finding the minimum of examined function.In our case the examined function is coefficient of pressure recovery c p , defined: where p 1 is static pressure at inlet, p 2 is static pressure at outlet; ȡ is density of liquid, v 1 is stream velocity at inlet.

Nelder-Mead method
Nelder-mead method, known as a simplex method, is a very popular method that does not use derive of the examined functions.
A method is described for the minimization of a function of n variables, which depends on the comparison of function values at the (n+1) vertices of a general simplex, followed by the replacement of the vertex with the highest value by another point.The simplex adapts itself to the local landscape, and contracts on to the final minimum.The method is shown to be effective and computationally compact [3].
Example of Nelder-Mead method: finding a global minimum of function: ).The path of simplex is shown in Figure 2.

Choice of initial simplex and parameters
The first step of the algorithm is a creation of regular simplex with side length a, using following definition: where x 1 is specified vertex, e k are unitary basis vectors and where n is a dimension.Important part of algorithm is also an ending tolerance criterion: where İ is required precision (tolerance).Algorithm consists of four basic operations that depend on chosen parameters: Parameter values are often chosen as follows:

Vertices ordering
The next step involves ordering according to the values of objective function at vertices: after that comes control of tolerance criterion, if epsilon is reached, the end of algorithm will occur.First, we must calculate ‫ݔ‬ҧ -the centre of gravity of the centroid (Figure 3):

Reflection
second, coordinates of reflected vertex x r are computed (Figure 3): finally, if ݂ ଵ ݂ ݂ , then vertex x r is accepted as replacement of x n+1 .Also function value in x r is calculated ݂ ൌ ݂ሺ‫ݔ‬ ሻ.Expansion only takes place when ݂ ݂ ଵ .Then coordinates of expanded vertex are computed:

Expansion
finally, if ݂ ൏ ݂ then vertex x e is accepted as replacement of x n+1 (else x r is).Also function value in x e is calculated ݂ ൌ ݂ሺ‫ݔ‬ ሻ.Contraction occurs when a given condition ݂ ݂ is satisfied.Depending on which vertex is better suited (x r or x n+1 ), comes external or internal contraction.

Contraction
Coordinates of internally contracted vertex are computed (if ݂ ݂ ାଵ ): then if ݂ ൏ ݂ ାଵ , computed vertex x ci is accepted as replacement of x n+1 , else coordinates of externally contracted vertex (if ݂ ൏ ݂ ାଵ ) are computed: EFM 2015 In this step of Nelder-Mead algorithm is an edge reduction of initial simplex:

Reduction
where v i , …, v n+1 (v 1 = x 1 ) are vertices of a new reduced simplex (Figure 6.).Notice: After each step of Nelder-Mead algorithm the control of tolerance criterion must be checked.If is not satisfied, algorithm returns to step called vertices ordering.Adjoint solver is a specialized tool that extends the scope of the analysis provided by a conventional flow solver by providing detailed sensitivity data for the performance of a fluid system [7].

Adjoint solver
Optimization via adjoint solver consists of several main steps -Figure 7.
There are two types of adjoint approachescontinuous and discrete.Ansys Fluent uses only the discrete one.

Discrete adjoint approach theory
A detailed description of discrete adjoint approach may be found in [6].
The adjoint solver needs a fully converged flow simulation, from which the gradients of objective function are derived.
The gradient of objective function ‫ܬ‬ ൌ ‫ܬ‬ሾܺ ǡ ሺܺ ሻሿ is [6]: where ܺ is vector of design variables (in [6] is reduced to scalar) and fraction is referred as flow sensitivities (Q is vector of conserved flow variables).This gradient could be rewritten in form [6]: where ܴ ൌ ܴሾܺ ǡ ܳሺܺ ሻሿ is the discretized RANS residual vector, ɗ represents vector of adjoint variables and is obtained by solving the linear system of equations (known as adjoint equation) [6]:

Overview
Initial diffuser design which was used in mathematical optimization methods (Nelder-Mead, adjoint solver) is shown in Figure 8. Table 1.displays types of boundary conditions, their locations, magnitudes and turbulence intensity values.Notice: Value of velocity at inlet is constant and does not match with real velocity profile behind the runner of swirl turbine.
For initial (Figure 8.) design was created computational mesh (Figure.9).Value of pressure recovery coefficient for this initial shape diffuser design was 0.736824.Shape optimization using Nelder-Mead algorithm consists of several steps (Figure 10.).Very important part is the first step, in which we define the criterial function (c p in our case), geometry (Figure 8.) and boundary conditions (Table 1

.).
Preprocessing involves creation of geometry and computational mesh.Everything is handled in Gambit software using programmed script.
The next two steps -CFD solver and postprocessortake place in Ansys Fluent, which also uses programmed script.Results (Fluent outcome) are later saved in text files due to easier access to data of examined variables.
The rest of the loop is realized in software called Matlab (including script calling from Gambit and Fluent).
Notice: The geometry change occurs only between vertices E and D.

Optimization using Bezier curves
In the next step of optimization, Bezier curves were used as restriction of wall E-D.

First order Bezier curve
First order Bezier curve is characterized only by one point of control polygon -P 1 .
Number of iterations in this approach was 52.Value of coefficient of pressure recovery c p reached up to 0.810309.Coordinates of control point P 1 were computed during Nelder-Mead algorithm as follows:

Adjoint solver
Adjoint solver is incompatible with an axis boundary condition in Fluent.Due to this fact, the working computational domain was created entirely in 3D (Figure 15.).After optimization in 3D, shape of diffuser was transformed back into 2D.In [10] author recommends polyhedral type of mesh (Figure 16.) to ensure better mesh behaviour during mesh morphing.Types of boundary conditions, their locations, magnitudes and turbulence intensity values used in Fluent software during optimization via adjoint solver are the same as in Table 1.
Simulation objective of adjoint solver was set to minimize static pressure at inlet, which is consistent with minimization of c p .
The next step of adjoint solver optimization approach (Figure 7.) involves shape modification -mesh morphing.Morphing takes place, when sensitivity field is computed (through adjoint equations).This morphing procedure has a double role.Firstly, it smooths the surface sensitivity field and secondly provides smooth distortion in boundary and interior mesh [10].In adjoint solver rectangular volume encloses boundaries which will be morphed.An array of control points is then distributed in this volume.These points in combination with the local mesh coordinates will define the proper movement both on the boundary and in the interior mesh.The local coordinate system relies on the properties of Bernstein polynomials [10].Bernstein polynomials are the basis of Bezier curves.
The most important information, which adjoint solver offers are optimal displacements (Figure 17.).Vectors of optimal displacement show location and intensity of diffuser shape change.

Conclusions
The values of c p obtained in both optimization methods -Nelder-Mead and adjoint solver -are summarized in Table 2.
Optimized shape of diffuser via Nelder-Mead method has 10.43 % higher pressure recovery than the original (initial) shape.This value was achieved through the use of 3. order Bezier curve as a wall constraint.The benefits of this method include: achieved value of c p and ability to control changed boundaries.Disadvantages include: the need to use additional software (Matlab) and the higher number of loops to achieve the final value of c p .
Optimized shape of diffuser via adjoint solver in Fluent has 10.25 % higher pressure recovery than the original (initial) shape, but this shape had extensively corrugated boundary.Problem was solved by polynomial approximation of vertices that form the diffuser.Value of c p increased by 0.04 %.The main advantages of this method are: the implementation in Ansys Fluent and shape modification is also present in the solver.Disadvantages include mainly the absence of better restriction of changed boundaries and a large number of settings, which adjoint solver uses.
Comparison of execution time for both methods has not been done due to use of two different computers (clusters).
It should be also mentioned that the value of velocity at inlet was constant and does not correspond with the real velocity profile behind the runner of swirl turbine.
Coefficient of pressure recovery c p diffuser increases as a result of weaker velocity gradients and shear layers (it was also confirmed by theoretical analysis [11]).

Figure 9 .
Figure 9.Initial 2D domain mesh Complete setting of the chosen discretization schemes is listed in the following chart.Pressure: Standard Momentum: Quick Turbulent diss.rate: Second order upwind Turbulent kinetic energy: Second order upwind

Figure 11 .
Figure 11.shows the change of diffuser geometry (red line).

Figure 17 .
Figure 17.Optimal displacements -vectors After 30 loops optimization was stopped, 3D geometry was transformed back to 2D and 2D geometry was recalculated using the same CFD settings as in Nelder-Mead algorithm.Coefficient of pressure recovery in this case reached up to value 0.812354 (22.loop -Figure 18.).

Figure 18 .
Figure 18.Adjoint solver approachDue to the absence of better restriction of shape (for example: using only the first order Bezier curve as a wall restriction) vertices of changed wall were later approximated in Excel using the 4. order polynomial (Figure19.).Coefficient of pressure recovery changed to 0.812614.