A CFD Based Application of Support Vector Regression to Determine the Optimum Smooth Twist for Wind Turbine Blades

Abstract

Computational fluid dynamics (CFD) is a powerful tool to estimate accurately the aerodynamic loads on wind turbine blades at the expense of high requirements like the duration of computation. Such requirements grow in the case of blade shape optimization in which several analyses are needed. A fast and reliable way to mimic the CFD solutions is to use surrogate models. In this study, a machine learning technique, the support vector regression (SVR) method based on a set of CFD solutions, is used as the surrogate model. CFD solutions are calculated by solving the Reynolds-averaged Navier–Stokes equation with the k-epsilon turbulence model using a commercial solver. The support vector regression model is then trained to give a functional relationship between the spanwise twist distribution and the generated torque. The smooth twist distribution is defined using a three-node cubic spline with four parameters in total. The optimum twist is determined for two baseline blade cases: the National Renewable Energy Laboratory (NREL) Phase II and Phase VI rotor blades. In the optimization process, extremum points that give the maximum torque are easily determined since the SVR gives an analytical model. Results show that it is possible to increase the torque generated by the NREL VI blade more than 10% just by redistributing the spanwise twist without carrying out a full geometry optimization of the blade shape with many shape-defining parameters. The increase in torque for the NREL II case is much higher.

1. Introduction

Power extraction from wind energy is a sustainable and popular option in the renewable energy sector [1,2,3]. Its growing popularity is mainly due to the fact that fuel cost is getting more expensive in addition to the global warming problems based on the usage of fossil fuels [4]. Moreover, the cost of energy harvesting from wind has been significantly reduced during the last years [5]. Therefore, the budget for other costs such as the reinforcement of the blade structure may be relaxed so that less restrictive bending moment constraints are considered in aerodynamic shape optimization of the turbine blades compared to the designs made in the past.

The full design of a wind turbine is a multidisciplinary optimization problem that incorporates the aerodynamic loads, structural requirements, and the investment/operation costs [6,7]. However, such an optimization includes many constraints that make the problem highly difficult and complex. Studies by Yu and Kwon [8] and Imiela et al. [9] are a few examples in the literature for the multidisciplinary analysis of wind turbines using high fidelity methods. A common approach to simplify the complexity of such a design is to use low fidelity methods on one or more branches of the multidisciplinary research. Employing the blade element momentum (BEM) method instead of computational fluid dynamics (CFD) to predict the aerodynamic performance of the turbine blade or making modal analysis of the blade structure instead of computing the loads using a solver based on the finite element method (FEM) is commonly preferred. Low fidelity models are easily implemented using low computational resources. For example, the blade design study conducted by Giguère et al. [10] rook into account both aerodynamic and structural considerations based on the BEM method constrained with the analytical relations between the mass, moment of inertia, and strength. However, the turbine performance calculations based on the low fidelity methods generally lead to less accurate results [11].

Recent works on the aerodynamic performance optimization of wind turbines by solving the Reynolds-averaged Navier–Stokes (RANS) equations to compute the three-dimensional flowfields around the blades include the studies of Dhert et al. [11], Vorspel et al. [12], Economon et al. [13], Elfarra et al. [14], and Shrestha [15]. Dhert et al., Vorspel et al., and Economon et al. implemented the adjoint method to determine the optimum blade geometry for maximum torque by choosing the NREL VI rotor blade as the baseline. Dhert et al. applied a discrete adjoint method with about 250 design variables, while Vorspel et al. and Economon et al. employed a continuous adjoint method with nine and 50 shape variables, respectively. Elfarra et al. optimized the winglet shape to be mounted on the NREL VI rotor blade to maximize the power generation using a genetic algorithm. Shrestha also optimized the NREL VI rotor blade for maximum torque. In that work of Shrestha, the blade shape was characterized by 14 design variables for the discrete chord length and twist values along the span and optimized employing the adaptive single objective algorithm.

The researchers cited above focused on demonstrating CFD-based wind turbine optimization after validating their CFD model against the results of the NREL experiments. However, the validation of the turbulence models to understand their suitability for computing the flowfields around the wind turbine blades is an ongoing study in the wind energy community [11]. Currently, the common approach is to use the k-$\epsilon$ turbulence model for the attached flow conditions in wind energy applications [11]. However, some researchers have shown that the computed results based on the k-$\epsilon$ turbulence model do not agree well with the experimental results in the separated flow conditions [16,17]. Another important work is by Rogowski et al. [18] who presented CFD solutions based on the detached eddy simulation (DES) technique. DES is a hybrid approach that combines the best aspects of the Reynolds-averaged Navier–Stokes and large eddy simulation methods. The numerical results obtained by Rogowski et al. showed that DES captures well the aerodynamic characteristics in the case of separated flows at high Reynolds numbers. Fortunately, since the attached flow conditions are considered in this study, the k-$\epsilon$ turbulence model is selected for the compatibility with the other researchers’ CFD validation works against the NREL experiments [14,19,20,21].

In contrast to a few CFD-based wind turbine optimization studies, there are many applications of the BEM method to design wind turbines for an objective function [22,23,24,25,26,27,28,29,30]. Ma et al. [22] optimized the offshore wind turbine blades using the particle swarm optimization technique. They determined optimum chord and twist distributions in terms of a number of discrete values along the span. Their results showed about a 4% increase in the generated power. Another discrete distribution of chord and twist was optimized by Chaudhary and Prakash [23] with an increase of around 8% in the power coefficient. Erturk [24] modified the baseline NREL VI rotor blade by defining piecewise constant chord and twist distributions along the span. The optimized blades were observed to produce almost the same power as the baseline blade. However, the author reported that they were less expensive and easier to manufacture. Using an analytical approach for the optimum chord length and twist angle values at the sections in the spanwise direction, Tenghiri et al. [25] designed a small wind turbine blade that was expected to generate the maximum power. The ant colony optimization algorithm implemented by Tahani et al. [26] successfully increased the power output from a wind turbine by about 14%. In their optimization process, six chord distribution functions (polynomial or logarithmic), ten twist distribution functions (polynomial or exponential), and also, 12 airfoils were considered. Capellaro and Cheng [27] obtained the optimum discrete twist angle values such that the equilibrium torsional deflection led to the maximum power generation. The study by Liu et al. [28] presented a heuristic approach for optimizing the blade shape of fixed-pitch fixed-speed small wind turbines by linearizing the chord and twist angle distributions along the span. They achieved about 3% higher annual energy production. Polat and Tuncer [29] optimized the shape geometry of a wind turbine blade using the genetic algorithm for maximum power. They took the sectional chord length, the sectional twist, and the blade profiles at the root, mid, and tip regions of the blade as design variables. Their results showed an increase of around 10% in the power production. Another study that used the genetic algorithm to optimize the blade shape was performed by Selig and Coverstone-Carroll [30]. They maximized the annual energy production by determining the optimum chord and twist distributions based on movable spline supports and the blade pitch. Then, they examined the sensitivity of annual energy production to changes in the rotor blade length and peak rotor power. Their conclusion stated that if the blade radius and peak power were fixed constraints, wind turbine blades designed for high-speed wind sites were likely to be equally well designed for lower speed sites.

In order to account for the advantages of both the CFD and BEM methods, some researchers conducted wind turbine optimization by combining these two methods. In this way, the design accuracy and flexibility were improved while maintaining the low computational cost [31,32]. Another approach to reduce the computational cost of a CFD solution to a level as low as that of a BEM solution is to use a surrogate model such as the response surface method. For example, Zahle et al. [33] optimized the blade tip extension for maximum energy production using the response surface method based on the incompressible flow solutions by a Navier–Stokes solver.

Similarly, other surrogate models like machine learning tools may also be implemented. Once a set of CFD-based flow solutions has been obtained in terms of the given values of a number of design variables, a machine learning model was trained to obtain the flow solutions for the other values of the design variables. In this study, the support vector regression method (SVR), which was firstly introduced by Vapnik et al. [34], is implemented as the machine learning model. This method has been used by many researchers [35,36,37,38,39,40,41,42,43] in various engineering problems in the last two decades. As other machine learning tools, it is generally preferred as an approximation tool for the solutions of the problems requiring high computational resources (CPU time, floating operations, memory, etc.) [35,36,37,38,39]. A well-defined support vector regression model may provide a highly-accurate functional relationship between some arbitrary input-output pair by minimizing the approximation error bound [42,43]. Balabin and Lomakina [44] observed that SVR was superior in both extrapolating and interpolating problems over the artificial neural network method in terms of the stability, generality, and robustness of the final model. They also stated that the SVR model did not require a wider dataset for accurate prediction results. For this reason, the SVR method is selected in this study to build the regression models for the twist distributions.

SVR applications cited in the wind energy literature are mainly about the fault detection and condition monitoring of wind turbines [45,46,47,48,49,50,51]. The studies by Shamshirband et al. [52], Erfort et al. [53], and Mohandes et al. [54] were more directly related to wind turbine performance. Sharmshirband et al., using SVR, obtained the optimum combination of the wind speed and the rotor speed for the maximum torque generated by a wind turbine. Their results showed that SVR can be considered as a promising alternative for other surrogate models. Erfort et al. applied SVR to determine the aerodynamic loads quickly as a function of the shape geometry of wind turbine profiles. It was indicated in the comparative study by Mohendes et al. that SVR is better in the predicting wind speed compared to the multilayer perceptron neural networks.

In this study, continuous spanwise twist distributions are optimized for maximum torque. The generated torque values are estimated using the support vector regression based on a set of CFD solutions. These CFD-based flow solutions are obtained employing a three-dimensional Reynolds-averaged Navier–Stokes solver. The continuous and smooth variation of the twist along the blade span is defined using a cubic spline. The parameters defining the cubic spline are input into the SVR model, which outputs the torque. Since the SVR provides an analytical function for the generated torque, the optimum spline parameters are easily determined. The baseline blades to be optimized are selected as the National Renewable Energy Laboratory (NREL) Phase II [55,56] and Phase VI [57,58] rotor blades. The taper distribution of the blades is kept fixed while optimizing the twist. All CFD solutions are obtained assuming the blade is rigid. The influence of the pitching moment on the blade twist is not considered.

2. Methodology

2.1. Flow Solver

The commercial CFD software FINE/Turbo by NUMECA International was employed to compute the flowfields around the wind turbine blades. The FINE/Turbo is a three-dimensional, compressible, structured, multi-block finite volume solver, which is specially developed to simulate the internal or external flows for turbomachinery, including wind turbines [59]. It can simulate the incompressible flows using low Mach number preconditioning.

The solver FINE/Turbo solves the unsteady Reynolds-averaged Navier–Stokes (URANS) equations formulated in the rotating frame. The equations are solved for the absolute velocities so that the freestream flow is taken as uniform. The 3-dimensional URANS equation according to this approach is given below:

$$\frac{\partial }{\partial t}{\int }_{V}QdV+{\oint }_S(\mathbf{F}·\mathbf{n})dS-{\oint }_S({\mathbf{F}}_{\mathbf{v}}·\mathbf{n})dS={\int }_V{s}_{T}dV$$

(1)

where Q is the vector of conservative variables Favre-averaged over some time interval, T. T is sufficiently large compared to the time scales of the turbulent fluctuations, but also less than all other time-dependent effects. The difference between the averaged and the instantaneous values of a variable is called the turbulent fluctuation of this variable. Favre-averaging of a variable, $\varphi$, is defined as $\tilde{\varphi }=\overline{\rho \varphi }/\overline{\rho }$ where:

$$\begin{array}{cccc}\hfill \overline{\rho }& =\frac{1}{T}{\int }_t^{t+T}\rho (\overrightarrow{x},t)dt\hfill & \hfill \overline{\rho \varphi }& =\frac{1}{T}{\int }_t^{t+T}\rho (\overrightarrow{x},t)\varphi (\overrightarrow{x},t)dt\hfill \end{array}$$

(2)

In the above expressions, V is the control volume inclosing the flowfield and S is the surface, which surrounds the control volume, V. $\mathbf{n}$ is the normal direction of S, and $\overrightarrow{x}$ is the 3-dimensional coordinates of a point within V.

The conservative variables, Q, the $j^{\mathrm{th}}$ component of the inviscid fluxes, $\mathbf{F}$, the $j^{\mathrm{th}}$ component of the viscous fluxes, ${\mathbf{F}}_{\mathbf{v}}$, and the source term, $s_T$ are defined as follows:

$$\begin{array}{cc}\hfill Q& =\left[\begin{array}{c}\overline{\rho }\\ \overline{\rho }{\tilde{u}}^1\\ \overline{\rho }{\tilde{u}}^2\\ \overline{\rho }{\tilde{u}}^3\\ \overline{\rho }\widetilde{e_0}+k\end{array}\right]\hfill \end{array}$$

(3)

$$\begin{array}{cccc}\hfill F^j& =\left[\begin{array}{c}\overline{\rho }{\tilde{w}}^j\\ \overline{\rho }{\tilde{w}}^1{\tilde{w}}^j+\overline{p}{\delta }_{1j}\\ \overline{\rho }{\tilde{w}}^2{\tilde{w}}^j+\overline{p}{\delta }_{2j}\\ \overline{\rho }{\tilde{w}}^3{\tilde{w}}^j+\overline{p}{\delta }_{3j}\\ \overline{\rho }\widetilde{h_0}{\tilde{w}}^j+k{\tilde{w}}^j\end{array}\right]\hfill & \hfill F_v^j& =\left[\begin{array}{c}0\\ {\tilde{\tau }}^{1j}-{\tau }_T^{1j}\\ {\tilde{\tau }}^{2j}-{\tau }_T^{2j}\\ {\tilde{\tau }}^{2j}-{\tau }_T^{2j}\\ {\tilde{u}}^i{\tilde{\tau }}^{ij}-{\tilde{q}}^j+{\mathsf{\Theta }}_T^j\end{array}\right]\hfill \end{array}$$

(4)

where $i,j=1,2,3$ refers to the component in each coordinate. $\tilde{u},\tilde{w},\widetilde{e_0},\overline{p}$, and $\widetilde{h_0}$ are absolute velocity, relative velocity, total energy, pressure, and total enthalpy, respectively. k is the kinetic energy of turbulent fluctuations. ${\delta }_{ij}$ is the Kronecker delta, which is equal to one if $i=j$ or equal to zero $i\ne j$. ${\tilde{\tau }}^{ij}$ is the viscous stress tensor, while ${\tau }_T^{ij}$ is the Reynolds stress tensor. ${\tilde{q}}^j$ is the $j^{\mathrm{th}}$ component of the heat flux vector. ${\mathsf{\Theta }}_T$ consists of the turbulent heat flux and all other turbulent terms evolving from density-velocity correlations and triple velocity correlations of the turbulent fluctuations.

The viscous stress tensor was calculated based on Stokes’ hypothesis and assuming the air is a Newtonian fluid. Sutherland’s law was used to determine the value of the viscosity term in this tensor. The CFD software FINE/Turbo is capable of computing the Reynolds stress tensor using a number of different algebraic and differential turbulence models. In this study, a two-equation model, the k-$\epsilon$ turbulence model, was selected to simulate the turbulent flow around the wind turbine blades. The motivation for selecting this turbulence model is explained in Section 1.

$\overrightarrow{\omega }$ is the angular velocity of the relative frame of reference. The source term due to the rotation of the blades is given as:

$$\begin{array}{cc}\hfill s_T& =\left[\begin{array}{c}0\\ -\overline{\rho }{(\overrightarrow{\omega }\times \overrightarrow{\tilde{u}})}^1\\ -\overline{\rho }{(\overrightarrow{\omega }\times \overrightarrow{\tilde{u}})}^2\\ -\overline{\rho }{(\overrightarrow{\omega }\times \overrightarrow{\tilde{u}})}^3\\ 0\end{array}\right]\hfill \end{array}$$

(5)

Finally, the URANS equations were closed under the perfect gas assumption, which relates the pressure and the internal energy.

The computational domain was discretized using a multi-block structured grid topology, which is described in the next section. Then, the solution was computed based on the finite volume method.

Grid Topology

The mesh around the blades was generated using the so-called O4Hgrid topology. As an example of the O4H grid topology, the schematic of a 5-block mesh is given in Figure 1. The description of the figure is as follows:

• An O block around the blade
• An H block upstream the leading edge of the blade
• An H block downstream the trailing edge
• An H block up to the blade section
• An H block down to the blade section

2.2. Boundary Conditions

The blade, hub, and shroud components of the wind turbine were modeled as solid boundaries on which adiabatic and nonslip boundary conditions were applied. The nonslip condition was satisfied by setting the fluid velocity on the solid surface equal to the blade velocity. The flow variables on the far-field boundary were estimated using the Riemann invariants based on the freestream values taken from the experimental data [55,57]. The far-field value of the turbulent kinetic energy, k, was calculated assuming a 1% freestream turbulent intensity. A 1/1 eddy viscosity ratio (the ratio of the turbulent viscosity to the molecular dynamic viscosity) at the far-field was used to predict the far-field value of the turbulent dissipation, $\epsilon$. This 1/1 eddy viscosity ratio was in the value range suggested by the developers of the solver FINE/Turbo [59]. A schematic for the flow domain is shown in Figure 2.

2.3. Unsteady Aerodynamics Experiments by the National Renewable Energy Laboratory

The National Renewable Energy Laboratory (NREL) conducted the unsteady aerodynamics experiments for quantifying the full-scale 3D aerodynamic behavior of horizontal axis wind turbines [55,57]. The experimental setup was installed to characterize rotating blade aerodynamic performance, structural responses, and atmospheric inflow conditions. Many blade configurations were tested and reported under different phase and sequence names.

The experimental results obtained by NREL have been used by many researchers to validate their numerical solutions for the aerodynamic and structural analyses. Especially, the NREL VI blade has been taken into consideration as the baseline blade by most of the studies on the blade shape optimization using the three-dimensional CFD solutions [11,12,13,14,15]. There are also recent works by some researchers who have validated their CFD solutions against the experimental results of the NREL II blade, which has a simpler geometrical shape [60,61]. These two popular blades are also selected in this study. The details of the NREL Phase II and NREL Phase IV Sequence H experimental configurations are given in the next section.

2.3.1. National Renewable Energy Laboratory II Rotor Blade

These experiments [55] were carried out at NREL’s National Wind Technology Center located in Colorado. The wind turbine of this phase was a downwind, three-bladed, horizontal axis wind turbine. The blades were untapered and untwisted. This wind turbine had the S809 airfoil as the blade profile from root to tip with a span length of 5.029 m and a constant chord length of 0.4572 m. A pitch of 12 deg was assigned to the blade at the $30\%$ chord location from the leading edge. The blades were attached to the hub through a circular section. The airfoil thickness was $43.0\%$ chord length at $14.4\%$ blade span in the radial direction, where the transition from the circular section ended. The thickness decreased linearly from $43.0\%$–$20.95\%$ until the $30\%$ span location. Outboard of 30% span, the thickness was maintained constant at a $20.95\%$ chord length value. The rotational speed of the NREL II rotor blade was 72 rpm. An illustration of a single blade of the NREL II rotor is given in Figure 3.

2.3.2. National Renewable Energy Laboratory VI Rotor Blade

These experiments [57] were carried out at the NASA Ames wind tunnel facility in California. The Phase VI experiments included several tests and sequences. In this study, Sequence H was selected as the baseline configuration. The wind turbine of this sequence was an upwind, two-bladed, horizontal axis wind turbine. The blades were tapered and twisted. This wind turbine had the S809 airfoil as the blade profile from root to tip with a span length of 5.029 m. A pitch of 4.815 deg was assigned to the blade at the $30\%$ chord location from the leading edge. The blades were attached to the hub through a circular section. The airfoil transition started at $17.6\%$ of the blade span in the radial direction. The taper and twist axes were located at the $30\%$ chord location from the leading edge. The taper and twist distributions of this blade were optimized based on a trade-off design [58]. The rotational speed of the NREL VI Sequence H rotor blade was 72 rpm. Figure 3 also shows a single blade of the NREL VI rotor.

2.4. Support Vector Regression

The general form of SVR to approximate the function $y=y(\overrightarrow{x})$ is expressed as follows:

$$y^{*}(\overrightarrow{x})=\langle \overrightarrow{w},\overrightarrow{\varphi }(\overrightarrow{x})\rangle +b$$

(6)

where $\langle .,.\rangle$ denotes the dot product, $\overrightarrow{x}$ is the vector of input variables, $y^{*}$ is an approximation function to the target function y, $\overrightarrow{w}$ is the weight vector, $\overrightarrow{\varphi }$ is a vector-valued function of $\overrightarrow{x}$, and b is a constant. In the literature, $\overrightarrow{\varphi }$ and b are respectively called the (non-linear) feature mapping function and the bias.

There were two aims while building the SVR model. The first aim was to determine the approximation function, $y^{*}(\overrightarrow{x})$, which had at most $\epsilon$ deviation from the actual target, $y(\overrightarrow{x})$. The second aim was to make $y^{*}(\overrightarrow{x})$ as flat as possible. Therefore, the following optimization problem is solved:

$$\begin{array}{cc}\underset{\overrightarrow{w},{\xi }_i^{+},{\xi }_i^{-}}{\mathrm{minimize}}\hfill & \hfill \frac{1}{2}{\parallel \overrightarrow{w}\parallel }^2+C\sum _{i=1}^m\left({\xi }_i^{+}+{\xi }_i^{-}\right)\\ \mathrm{subject}\mathrm{to}\hfill & \hfill \langle \overrightarrow{w},\overrightarrow{\varphi }(\overrightarrow{x_i})\rangle +b-y_i\le \epsilon +{\xi }_i^{+}\\ & \hfill y_i-\langle \overrightarrow{w},\overrightarrow{\varphi }(\overrightarrow{x_i})\rangle -b\le \epsilon +{\xi }_i^{-}\\ & \hfill {\xi }_i^{+},{\xi }_i^{-}\ge 0,i=1,\dots ,m\end{array}$$

(7)

where $\overrightarrow{x_i}$ and $y_i$ denote the $i^{\mathrm{th}}$ input-output pair in the training dataset. m is the number of data pairs in the entire set or in a subset of the entire set, depending on the training algorithm. $\epsilon$ is called the loss function and is a model parameter that must be supplied before solving the minimization problem in Equation (7). The penalty parameter, $C>0$, is also a model parameter and must be supplied a priori as well. It determines the trade-off between the flatness of $y^{*}(\overrightarrow{x})$ and the amount up to which deviations larger than $\epsilon$ are tolerated. Finally, ${\xi }_i^{+}$, ${\xi }_i^{-}$ are the slack variables to provide a solution to Equation (7) in the feasible domain by coping with the $\epsilon$ constraint. The explanations and comments about $\overrightarrow{\varphi }(\overrightarrow{x})$ and b will be given at the end of this section.

Since the Lagrange dual form is simpler (in terms of constraint handling), it is solved rather than Equation (7) itself:

$$\begin{array}{cc}\underset{{\alpha }_i^{+},{\alpha }_i^{-}}{\mathrm{maximize}}\hfill & \hfill -\frac{1}{2}\sum _{i,j=1}^m({\alpha }_i^{+}-{\alpha }_i^{-})({\alpha }_j^{+}-{\alpha }_j^{-})k(\overrightarrow{x_i},\overrightarrow{x_j})+\sum _{i=1}^m({\alpha }_i^{+}-{\alpha }_i^{-})y_i-\sum _{i=1}^m({\alpha }_i^{+}+{\alpha }_i^{-})\epsilon \\ \mathrm{subject}\mathrm{to}\hfill & \hfill \sum _{i=1}^m({\alpha }_i^{+}-{\alpha }_i^{-})=0,i=1,\dots ,m\\ & \hfill 0\le {\alpha }_i^{+},{\alpha }_i^{-}\le C,i=1,\dots ,m\end{array}$$

(8)

where ${\alpha }_i^{+}$ and ${\alpha }_i^{-}$ are the Lagrange multipliers. $k(\overrightarrow{x_i},\overrightarrow{x_j})=\langle \overrightarrow{\varphi }(\overrightarrow{x_i}),\overrightarrow{\varphi }(\overrightarrow{x_j})\rangle$ is called the kernel function. As seen from its definition, it is symmetric and positive definite. The Lagrange multipliers are determined by satisfying the Karush–Kuhn–Tucker conditions assuming that the kernel function, $k(\overrightarrow{x},\overrightarrow{x_i})$, and C are both supplied. Finally, SVR is obtained as follows:

$$y^{*}(\overrightarrow{x})=\sum _{i=1}^m({\alpha }_i^{+}-{\alpha }_i^{-})k(\overrightarrow{x},\overrightarrow{x_i})+b$$

(9)

It is seen from Equation (9) that the weight vector $\overrightarrow{w}$ and the feature mapping function $\overrightarrow{\varphi }(\overrightarrow{x})$ were eliminated from the original definition of SVR in Equation (6). The most popular kernel function in the literature is the Gaussian kernel [62]. It has only one model parameter: $\gamma$. This kernel function is defined as follows:

$$k(\overrightarrow{x},\overrightarrow{x_i})=e^{-\gamma \parallel \overrightarrow{x}-\overrightarrow{x_i}{\parallel }^2}$$

(10)

In this study, the Gaussian kernel was chosen as the kernel function when building the SVR models due to its popularity and simplicity. The last parameter to be determined is the bias, b. The detailed information about choosing b was given by Keerthi et al. [63]. As the simplest way, it may be chosen arbitrarily provided that the following condition is satisfied:

$$\begin{array}{cc}\hfill max\{y_i-\langle \overrightarrow{w},\overrightarrow{\varphi }(\overrightarrow{x_i})\rangle -\epsilon |{\alpha }_i^{+}<C\mathrm{or}{\alpha }_i^{-}>0,i=1,\dots ,m\}& \le b\le \\ \hfill min\{y_i-\langle \overrightarrow{w},\overrightarrow{\varphi }(\overrightarrow{x_i})\rangle -\epsilon |{\alpha }_i^{+}>0\mathrm{or}{\alpha }_i^{-}<C,i=1,\dots ,m\}& \end{array}$$

(11)

In this work, the open source library LIBSVM by Chang and Lin [64] was employed to build the SVR models. The parameters $\epsilon$, C, and $\gamma$ that must be supplied a priori for building the model were selected such that the mean squared error between $y(\overrightarrow{x})$ and $y^{*}(\overrightarrow{x})$ was minimum based on a cross-validation study. For this purpose, the strategy stated by Yan et al. [62] was used. A brief explanation of the cross-validation is also given in Section 4.

3. Optimization

The SVR model was trained for the twist parameters explained in Section 4 as input variables and the torque values calculated based on the CFD solutions as output variables. Once the model was built, we had an analytically-defined fitting function (Equation (9)) between the input and the output. The extremum of this known analytical function, that is the maximum generated torque value, can be easily found by making its first derivative vanish since there is no constraint on the twist parameters:

$$\frac{d{y}^{*}(\overrightarrow{x})}{d\overrightarrow{x}}=\sum _{i=1}^m({\alpha }_i^{+}-{\alpha }_i^{-})\frac{dk(\overrightarrow{x},\overrightarrow{x_i})}{d\overrightarrow{x}}=0$$

(12)

or:

$$\begin{array}{ccc}\hfill \sum _{i=1}^m({\alpha }_i^{+}-{\alpha }_i^{-})\frac{dk(\overrightarrow{x},\overrightarrow{x_i})}{d{x}^{(1)}}& =& 0\hfill \\ \hfill ⋮\\ \hfill \sum _{i=1}^m({\alpha }_i^{+}-{\alpha }_i^{-})\frac{dk(\overrightarrow{x},\overrightarrow{x_i})}{d{x}^{(4)}}& =& 0\hfill \end{array}$$

(13)

where $x^{(1)},\dots ,x^{(4)}$ are the components of the input vector $\overrightarrow{x}$, which are the parameters of the cubic spline to define the spanwise twist distribution (see Section 4 for details):

$$\begin{array}{cc}\hfill \overrightarrow{x}& =\left[\begin{array}{c}{x}^{(1)}\\ x^{(2)}\\ x^{(3)}\\ x^{(4)}\end{array}\right]=\left[\begin{array}{c}{\theta }_{root}\\ {\theta }_{mid}\\ {\theta }_{tip}\\ \frac{d\theta }{dr}{|}_{root}\end{array}\right]=\left[\begin{array}{c}\mathrm{Root}\mathrm{twist}\mathrm{angle}\\ \mathrm{Mid}\mathrm{twist}\mathrm{angle}\\ \mathrm{Tip}\mathrm{twist}\mathrm{angle}\\ \mathrm{Twist}\mathrm{slope}\mathrm{at}\mathrm{the}\mathrm{root}\end{array}\right]\hfill \end{array}$$

(14)

Since the system of equations in Equation (13) is non-linear due to the Gaussian kernel, $k(\overrightarrow{x},\overrightarrow{x_i})$ defined in Equation (10), it should be solved numerically. In this work, it was solved iteratively using Newton’s method.

4. Design of Experiment

Any model based on a machine learning algorithm is trained using a dataset that consists of input-output pairs. In this study, the input was a vector of four parameters, which defined the spanwise twist distribution using a cubic spline, while the output was the corresponding torque value calculated based on this twist. The number of the input-output pairs in a training dataset depends on the topology used to construct the design of experiment (DOE). A typical topology provides the list of the element-wise combinations of the input based on the minimum, maximum, and intermediate values of the input vector elements. In this work, the dataset for the SVR model was obtained using the Box–Behnken design of experiment [65]. The Box–Behnken method suggests 25 input-output pairs for a four-element input vector.

Table 1. Minimum, maximum, and intermediate input values to construct the Box–Behnken DOE.

Parameter	Minimum Value	Intermediate Value	Maximum Value
Root twist angle in degrees, ${\theta }_{root}$	5.0	15.0	25.0
Mid twist angle in degrees, ${\theta }_{mid}$	−2.5	5.0	12.5
Tip twist angle in degrees, ${\theta }_{tip}$	−5.0	0.0	5.0
Twist slope at the root in degrees per meter, $\frac{d\theta }{dr}{|}_{root}$	−20.0	0.0	10.0

gives the minimum, maximum, and intermediate values of the input vector elements for both of the baseline blade cases. According to these values, a total of 50 CFD simulations was carried out: 25 simulations for the NREL II case and 25 simulations for the NREL VI case. It should be noted that two separate SVR models were generated for two baseline blade cases.

When building the SVR model for a baseline blade case, the dataset was split into five sections (so-called folds), each of which included five CFD results. Then, each fold was used as a testing set. Finally, the average of the approximation errors of five testing sets was considered as the performance of the SVR model based on the parameters $\epsilon$, C, and $\gamma$. Such a validation method for a machine learning algorithm is called as K-fold cross-validation, where K = 5 in this study.

As stated earlier and also seen from Table 1, there were four parameters that defined the cubic spline-based smooth twist variation along the span:

• Twist angle value at the root, ${\theta }_{root}$
• Twist angle value at the mid-span, ${\theta }_{mid}$
• Twist angle value at the tip, ${\theta }_{tip}$
• Spanwise rate of change of the twist angle at the root, $\frac{d\theta }{dr}{|}_{root}$

where the root and tip locations of the blade span are shown in Figure 3. The mid-span location is the middle point between the root and the tip. r denotes the radial direction along the span.

Cubic Spline-Based Twist Distribution

A cubic spline is a spline constructed of piecewise third-degree polynomials that pass through a set of control points, which are generally called knots. First and second derivatives of a cubic spline are continuous at each knot, providing the smoothness of the data. Moreover, when constructing a cubic spline, the first (or second or third) derivatives of the spline at the endpoints must also be known. The convention is to select the second derivatives at the first and second knots as zero, which is called the natural cubic spline.

In this study, three knots ($r_{root}$, $r_{mid}$, and $r_{tip}$) were used to construct the cubic spline together with three corresponding twist angle values (${\theta }_{root}$, ${\theta }_{mid}$, and ${\theta }_{tip}$). Unlike the conventional (natural) cubic spline, the first derivatives at the root and the tip were provided. The first derivative at the first knot, that is $\frac{d\theta }{dr}{|}_{root}$, was included as the fourth parameter in the input vector, while the first derivative at the last known on, that is $\frac{d\theta }{dr}{|}_{tip}$, was taken as zero.

A mathematically-expressed summary of the cubic spline used to define the spanwise twist distribution $\theta =\theta (r)$ is as follows:

$$\theta (r)=\left\{\begin{array}{cc}{\theta }_1(r)=a_1{(r-r_{root})}^3+b_1{(r-r_{root})}^2+c_1(r-r_{root})+d_1\hfill & \mathrm{if}{r}_{root}\le r\le r_{mid}\hfill \\ {\theta }_2(r)=a_2{(r-r_{mid})}^3+b_2{(r-r_{mid})}^2+c_2(r-r_{mid})+d_2\hfill & \mathrm{if}{r}_{mid}\le r\le r_{tip}\hfill \end{array}\right.$$

(15)

where the unknown coefficients $a_1$, $b_1$, $c_1$, $d_1$, $a_2$, $b_2$, $c_2$, and $d_2$ are determined according to the following conditions:

$$\begin{array}{ccc}{\theta }_1(r_{root})& =& {\theta }_{root}\\ {\theta }_1(r_{mid})& =& {\theta }_{mid}\\ {\theta }_2(r_{mid})& =& {\theta }_{mid}\\ {\theta }_2(r_{tip})& =& {\theta }_{tip}\\ d{\theta }_1{/dr|}_{mid}& =& d{\theta }_2{/dr|}_{mid}\\ d^2{\theta }_1/d{r}^2{|}_{mid}& =& d^2{\theta }_2/d{r}^2{|}_{mid}\\ d{\theta }_1{/dr|}_{root}& =& {d\theta /dr|}_{root}\\ d{\theta }_2{/dr|}_{tip}& =& 0\end{array}$$

(16)

Equation (16) leads to a linear system for eight unknowns, which is easily solved.

5. Validation Study

The CFD software FINE/Turbo was previously validated against the results of the NREL II and NREL VI experiments [66,67]. A validation study was also conducted in this work. The results of the grid sensitivity analysis performed before the validation are given in

Table 2. Mesh study results.

Mesh	Number of Cells	Torque (Nm)	Wall Clock Time of Computation
1	$0.8\times {10}^6$	754.4	6 min
2	$1.6\times {10}^6$	776.4	12 min
3	$3.2\times {10}^6$	781.6	26 min
4	$4.8\times {10}^6$	784.0	35 min
5	$6.9\times {10}^6$	784.1	55 min
6	$10.7\times {10}^6$	784.4	78 min

. For this purpose, six different meshes at varying resolutions were studied. The number of cells on the finest mesh was about 11 million, while the coarsest mesh had around 800 thousands cells. It is seen from Table 2 that the mesh with about seven million cells gave almost the same results as the finest mesh. Therefore, in this study, the mesh settings used in the fight mesh were employed for all analyzed cases including validation. In fact, these were the same spatial and temporal resolutions stated previously by Kaya and Elfarra [67]:

The number of nodes in each mesh was around seven million. In order to have a $y^{+}$ value close to one, which was suitable for the implemented k-$\epsilon$ turbulence model and the freestream Reynolds number of Re $={10}^6$, the first nodes had a fixed distance of $3\times {10}^{-6}$ m in the orthogonal direction to the solid boundaries. One hundred forty five nodes were employed to define each two-dimensional blade section, that is airfoil profiles along the span with the normalized grid spacing of $1.5\times {10}^{-3}$ on the leading edge and $2.0\times {10}^{-4}$ on the trailing edge. These normalized values were based on the sectional chord length.

All the meshes used in this study were generated for a single blade, imposing the periodic condition to account for the other blade(s). A typical mesh generated according to these settings is shown in Figure 2, while Figure 4 gives the surface mesh around the rotor blade.

Figure 5 and Figure 6 show the chordwise variations of the pressure coefficients calculated at different span locations of the baseline blades NREL II and NREL VI. The validation study was carried out for the wind speed of $w_{\infty }=7.2$ m/s as in the NREL II experiment and for $w_{\infty }=7.0$ m/s as in the NREL VI experiment. It is seen from the figures that the computed values for both blades agreed well with the experimental data.

6. Results

The analyzed cases for the NREL II blade were investigated at a wind speed of 7.2 m/s as given in the NREL II experiment. For the NREL VI cases, the solutions were obtained at 7 m/s, which was one of the reported wind speeds in the experiment. The turbulent flowfields were solved using the k-$\epsilon$ model as stated by Elfarra et al. [14]. Low Mach number preconditioning for compressible solvers was employed to overcome the incompressible nature of the flows around the wind turbine blades. The steady-state flow solutions were obtained using cell size-based local time stepping.

The flowfields for all the cases were computed in parallel using 16 processors. A typical CFD solution with about a sixth-order convergence requires almost an hour of wall clock time. Figure 7 shows the convergence history of residual and generated torque for an analyzed case. A total of 50 CFD solutions (25 for NREL II, 25 for NREL VI) was computed as stated in Section 4. Consequently, it took about 50 h to complete the samples in the dataset. The time spent for building the SVR models and obtaining the optima was in the order of minutes, that is negligible. Therefore, the total time cost for the machine learning-based solutions of the optimization problems in this study was approximately two days. However, instead of being based on a machine learning model, a classical or modern optimum seeking technique based on iterative function evaluations would be much more expensive in terms of the time consumption. For example, the conjugate gradient method needed gradient calculations and line searches in each optimization step, which meant a large number of function evaluations. It should be also noted that as the number of optimization parameters increased, the number of the function evaluations increased as well. Similarly, since a genetic algorithm theoretically converges after infinitely many generations, the number of the function evaluations may be quite high, although there is no need for the gradient calculation. In that algorithm, the size of the population (function evaluation) in each generation is the main cause for the time consumption because it is generally in the orders of hundreds. Machine learning methods like SVR appear to be a good remedy to reduce the computational costs for the solutions of the optimization problems in which a function evaluation requires a highly long computation duration [35,36,37,38,39].

SVR models were separately constructed for the NREL II and NREL VI cases. Then, the optimum parameters were determined for the maximum torque. The maximum torque predictions at the optimized twist conditions of the NREL II and NREL VI blades are given in

Table 3. Maximum torque predictions for the optimum twist distributions.

Blade	Torque (Nm)	Increase in Torque (%)
Baseline NREL II	362	-
NREL II with optimum twist	835	131
Baseline NREL VI	784	-
NREL VI with optimum twist	856	9.2

. The baseline values are also reported to compare with the maximum ones. It is seen from the table that the maximum torque was more than twice the baseline value for the NREL II case. In fact, such a high increase was expected because the baseline NREL II was an untwisted blade. The effect of spanwise twist was significantly observed in this optimization process.

The increase in torque was almost 10% for the NREL VI case. It should be noted that the twist variation of the baseline NREL VI blade was stated to be optimum in terms of a trade-off design [58]. The result showed that it was possible to increase the torque of the NREL VI blade just by redistributing the twist. The present maximum torque value was quite significant compared to the torque reported by Economon et al. [13] in their shape optimization work for the NREL VI blade. They obtained an increase of 4% in the torque generation after they optimized the shape geometry of the NREL VI blade using a continuous adjoint method with 50 optimization variables. Another optimization study for maximizing the torque generation of the NREL VI blade was by Shrestha [15] who obtained approximately a 6% increase in the torque. Shrestha used 14 optimization variables for both taper and twist variations along the span. Applying a discrete adjoint method with 12 variables (one variable was the pitch angle) to define the twist, Dhert et al. increased the torque of the NREL VI blade by about 6% [11]. It is clearly seen that the present study, which only used four variables, was able to obtain more torque generation. It is worthwhile mentioning that Dhert et al. observed a maximum torque value about 22% higher than the baseline when they increased the number of optimization variables from 12 to 252, including variables to define the sectional airfoil profiles.

Table 4. Validation of SVR predictions by CFD computations.

Blade	Torque by SVR (Nm)	Torque by CFD (Nm)	Difference (%)
NREL II with optimum twist	835	832	0.36
NREL VI with optimum twist	856	853	0.35

provides information about the validation of the SVR predictions. In order to validate the SVR predictions, extra CFD solutions were computed for the optimized blades. As seen from the table, SVR successfully estimated the torque values for the optimized twist parameters.

The optimum parameters determined by the SVR models of both NREL II and NREL VI cases are shown in

Table 5. Optimum parameters for the cubic spline-based twist distribution.

Parameter	Optimized NREL II	Optimized NREL VI
Root twist angle in degrees, ${\theta }_{root}$	12.1	4.2
Mid twist angle in degrees, ${\theta }_{mid}$	−9.6	−2.9
Tip twist angle in degrees, ${\theta }_{tip}$	−11.4	−5.0
Twist slope at the root in degrees per meter, $\frac{d\theta }{dr}{|}_{root}$	−22.6	−2.9

. The corresponding twist distributions based on the cubic spline interpolation are depicted in Figure 8. The interesting observation about the results is the fact that the twist slope at the root was very high in the NREL II case where the blade was untapered, in contrast to the low twist slope at the root for the NREL VI case, which had a designed taper distribution. Moreover, for the NREL VI case, the twist angle at the root was significantly less compared to the baseline value.

The chordwise distributions of the pressure coefficient for the optimized NREL II blade are shown in Figure 9. The distributions are plotted at various span locations and compared to the baseline distributions. The pressure coefficient distribution at the 30% span location was almost the same for both baseline and optimized twist conditions. This was due to the fact that the optimized twist distribution gave almost zero twist at this span location as the baseline distribution (see Figure 8 for NREL II). In the other span locations of the optimized blade, the pressure coefficient difference between the upper and lower surfaces was relatively higher, as expected.

Figure 10 shows the pressure coefficient distributions for the NREL VI blade with the optimum twist. The baseline distributions are also plotted for comparison. Taking into consideration that the difference between the twist angle values of the optimized and baseline blades was approximately five degrees as a spanwise average (See Figure 8 for NREL VI), the difference between the optimized and baseline pressure coefficient distributions was not as pronounced as the NREL II case. Yet, this difference was sufficient to increase the torque by 10%.

The optimization was carried out at a single wind speed for both NREL II and NREL VI. As stated earlier, the wind speed at which the twist optimization was conducted for NREL II was 7.2 m/s, while the wind speed for the NREL VI optimization was 7 m/s. The torque performance of both optimized blades were also investigated for other wind speeds. The investigation was done for the speed values ranging from 5 m/s–9 m/s. The result is given in Figure 11. As seen from the figure, although the optimization was carried out at a single wind speed, the torque of both optimized blades increased for the whole speed range. The only exception was the wind speed of the 9 m/s case of the NREL VI blade. Elfarra [66] showed that 9 m/s was the threshold speed at which the flow separation started for the NREL VI blade. Therefore, such a discrepancy in the flow regime may be the reason for this particular result, as stated by Volikas et al., Aksenov et al., and Rogowski et al. [16,17,18].

Previous works by other researchers on the CFD-based shape optimization of wind turbine blades either maximized the torque [11,13,14,15] or minimized the thrust [12]. However, torque generation may also be maximized while maintaining the baseline thrust value. For this purpose, new SVR models were built to output the thrust in addition to the SVR models outputting the torque. Then, the optimum twist parameters providing the maximum torque at the baseline thrust condition were determined. The results are given in

Table 6. Maximum torque for the baseline thrust condition.

Blade	Torque (Nm)	Thrust (N)	Increase in Torque (%)
Baseline NREL II	362	588	-
NREL II with optimum twist	398	588	9.9
Baseline NREL VI	784	1307	-
NREL VI with optimum twist	795	1307	1.4

. It is clearly seen from the table that if the baseline thrust value was maintained when optimizing the twist distribution for the maximum torque, then the increase in the torque was not as significant as in the unconstrained maximization of the torque generation. Corresponding twist variations are plotted in Figure 12. As seen from the figure, for the NREL II case where the blade was untapered, the twist angle first increased, then decreased along the span. For the NREL VI case, the twist distribution was almost the same except for the root. It was again observed that the optimum twist angle at the root was much less than the baseline value for the NREL VI case.

Since the SVR models were ready for both torque and thrust outputs, one may also obtain the maximum torque variation along the specified values of thrust. This approach is depicted in Figure 13. The figure shows also the baseline torque and thrust values as a single point in the plot. As observed from the figure, the torque was maximized at the expense of higher thrust. Therefore, it may be concluded that a wind turbine blade designed for maximum torque must be reinforced to withstand the increased thrust.

7. Conclusions

This study presents a framework to determine the optimum twist distribution for the wind turbine blades. The suggested method can also be applied for taper distribution or any other geometric definitions of the blade shape. The optimization was carried out for maximum torque. For this purpose, a machine learning approach based on the support vector regression was implemented. Since support vector regression models provide analytical functions for the input-output relations, the optimization problem was easily solved. Support vector regression models used in this study were built based on the three-dimensional flowfield solutions computed by a CFD solver.

The twist distribution was defined using a cubic spline with three knots at which the twist angle values were given. The knots were selected as the root, mid, and tip span locations. The fourth parameter to define the twist distribution was the rate of change of the twist angle at the root. Since the twist distribution was based on the cubic spline, it was smoothly continuous along the span.

NREL II and NREL VI rotor blades were chosen as the baseline blades to be optimized. All CFD solutions were computed for the rigid blade assumption. It was observed that the optimized twist distribution for the NREL II blade provided 130% more torque compared to the baseline distribution. The torque for the optimized NREL VI blade was determined to be 10% higher relative to the baseline blade. The optimum torque value obtained for the NREL VI case in this study was quite significant compared to the results given by other researchers who even used more optimization parameters to define the blade shape.

SVR predictions were also validated using extra CFD solutions at the optimum conditions in order to assess the accuracy of the models for the input, which was not in the dataset (design of experiment). It was observed that the torque values estimated by the SVR models were almost the same as the values given by the CFD solution. This is a good demonstration for the reliability of the suggested framework.

The optimization of the NREL II blade, which is an untapered blade, led to a high twist slope at the root, while a low twist slope was observed in the optimization of the NREL VI blade, which had a designed spanwise distribution of the chord length. Moreover, the optimized NREL VI blade had a significantly less twist angle value at the root compared to the baseline value.

The torque performance of both optimized blades was also investigated for other wind speeds. It was observed that the torque of both optimized blades increased for the whole speed range except for the 9 m/s speed of wind flowing through the NREL VI blade. The reason was thought to be the separated flow conditions at this speed.

CFD-based blade shape optimization studies by other researchers either maximized the torque or minimized the thrust. In this study, torque generation was also maximized, while not allowing an increase in the thrust value. However, it was observed that if the baseline thrust value was maintained when optimizing the twist distribution for the maximum torque, then the increase in the torque was not as significant as in the unconstrained maximization of the torque generation.

The maximum torque variation along the specified values of thrust was also obtained in this study. As expected from the observation mentioned above, the torque was maximized at the expense of higher thrust. Therefore, one may say that if a wind turbine blade were designed for maximum torque, then it must be reinforced to withstand the increased thrust.

The summarizing remark for the conclusion of this study is the fact that the suggested method based on SVR and CFD was reliable, fast, and easy for determining the blade shape geometry.