Development of a CFD-Based Wind Turbine Rotor Optimization Tool in Considering Wake Effects

Abstract

In the present study, a computational fluid dynamic (CFD)-based blade optimization algorithm is introduced for designing single or multiple wind turbine rotors. It is shown that the CFD methods provide more detailed aerodynamics features during the design process. Because high computational cost limits the conventional CFD applications in particular for rotor optimization purposes, in the current paper, a CFD-based 2D Actuator Disc (AD) model is used to represent turbulent flows over wind turbine rotors. With the ideal case of axisymmetric flows, the simulation time is significantly reduced with the 2D method. The design variables are the shape parameters comprising the chord, twist, and relative thickness of the wind turbine rotor blades as well as the rotational speed. Due to the wake effects, the optimized blade shapes are different for the upstream and downstream turbines. The comparative aerodynamic performance is analyzed between the original and optimized reference wind turbine rotor. The results show that the present numerical optimization algorithm for multiple turbines is efficient and more advanced than conventional methods. The current method achieves the same accuracy as 3D CFD simulations, and the computational efficiency is not significantly higher than the Blade Element Momentum (BEM) theory. The paper shows that CFD for rotor design is possible using a high-performance single personal computer with multiple cores.

1. Introduction

To reduce the cost of material and improve the aerodynamic efficiency of wind energy conversion systems, aerodynamic optimization of wind turbine rotors is important. Modern horizontal axis wind turbine (HAWT) rotors are aerodynamically optimized to efficiently extract wind energy. The physical limit to extract energy from wind is known as the Betz limit. The wind energy extraction process is through the reduction of wind speed passing the wind turbine rotor. The zero wake flow scenario will yield full energy absorption from wind, which certainly cannot be achieved. Since a wake always exists behind any wind turbine, and wind turbines are mostly clustered in a wind farm, the design and optimization of wind turbines in the wake situation are indeed necessary to represent some real-world wind farm cases.

From the theoretical to practical rotor aerodynamic design, the final products suffer from an accumulation of many aerodynamic and mechanical losses, such as tip loss, root loss, blade surface roughness change loss, blade shape manufacturing accuracy loss, drive train efficiency loss and wake effect loss. Despite of mechanical and manufacturing imperfections, most of aerodynamic optimizations have considered all the above-mentioned facts. The effect of wind turbine wake has been a very popular topic which is essentially important for wind farm layout optimization. However, in terms of rotor design itself, the wake effect is seldom included due to the inherent property of aerodynamic design tools. For example, the Blade Element Momentum (BEM) method is fast for optimization of single rotors, but it cannot deal with wake flow problems; the CFD-based methods are more accurate but practically not used for rotor optimization purposes.

In the field of wind turbine design and optimization, various simulation methods are involved for single rotors and for wind farms. For single rotor aerodynamic force calculations, the wind turbine aerodynamic modeling relies mainly on the following approaches: BEM [1,2,3,4,5,6]-based methods, Vortex Wake (VW)-based methods [7,8,9,10,11,12], and CFD-based methods [13,14,15,16]. The CFD methods either apply body-fitted grid or AD/AL/AS (Actuator Disc, Actuator Line, Actuator Surface) techniques. BEM has been the most efficient method which has been widely used to solve engineering problems. As the most efficient tool, BEM coupling with an optimization algorithm is mainly adopted as the method of wind turbine rotor design [13,17,18]. The VW methods are more flexible in the sense that it deals with the turbulent wake state that cannot be done with BEM. A study of BEM and VW methods are seen in [19] where general agreements are achieved for most of the considered cases but not for the high-thrust cases. However, VW methods are still not efficient for optimization purposes, and the accuracy in wake flow predictions affects the wind turbine performance to some extent. It is worthwhile to point out that the method to be used is relied on the matter of interest. From the multi-wind turbine optimization point of view, it is important to focus on both blade aerodynamics and turbulent wake. Therefore, BEM, VM and CFD based on curvilinear grid are commonly used for load calculations. However, rotor optimization requires a high computational demand, and so far CFDs are popularly used for wind turbine aerodynamic analysis. On the other hand, rotor aerodynamics has never been isolated with wake dynamics. However, existing models are individually used either for rotor design or wake study. Due to the inherent limitations of the BEM method, it cannot be directly applied to study wind turbine wakes. The more accurate methods, such as VM and CFD methods can be good choices to investigate both rotor aerodynamics and wake dynamics, but VM is not designed for long distance wake simulation and CFD is very expansive to study long range wake propagations, such that it is very hard to combine rotor aerodynamics with wake effects for rotor optimization purpose.

For a wind turbine cluster or wind farm, the velocity deficit causes a momentum loss in its wake before reaching the downstream wind turbines. Optimal design of the most upstream turbine will certainly increase the power product of the first turbine or perhaps the entire first row of a wind farm. However, for the second and third turbines downwind in the wake of the first turbine, the original design may not be the best configuration. Although this problem is realized through wind farm wake studies, but a practical design of wind turbines with wake effect is difficult either due to the model accuracy or high computational demand. The current work relates the rotor aerodynamics to the wake dynamics of multi-wind turbines. Therefore, modeling the velocity deficit is significantly important. Numerous studies were carried out in connection of large wind farm developments [20,21]. More recent work of Göçmen et al. [22] presented a good overview of various wake models that were developed at DTU. Some selected wake modelling approaches are listed in

Table 1. Selected wake models.

Class	Methodology	References
Class 1	Jensen	Jensen (1983)
	Dynamic Wake Meandering	Larsen (2007, 2009)
Class 2	Vortex Wake Method	Bareiss (1993), Voutsinas (1993), Wang (1998)
Class 3	Fuga	Ott (2011)
	Eddy Viscosity Model	Ainslie (1988)
Class 4	RANS/LES + curvilinear	Sørensen (2014), Thé (2017), Allah (2017), Abdulqadir (2017)
	RANS/LES + AD/AL/AS	Sørensen, Shen et al. (1998, 2002, 2009, 2010, 2012), Troldborg (2014)
	LBM/IBM	Deiterding (2015), Angelidis (2015), Zhu (2013)

. The methods are divided into four groups according to their computational time demand. The model of Jensen [23] and Larsen [24,25] belong to class 1 that are the fastest engineering models. Some modifications of class 1 models become necessary to better represent the flow physics. For example, Tian et al. [26,27] proposed a new Jensen model that uses a cosine wake shape instead of the top-hat shape for the velocity deficit. The VW models and its modifications are adopted from aircraft aerodynamics and well-suited for wind turbine aerodynamics in terms of accuracy and efficiency [8,9,28]. Efforts were also made to reduce the computational time using CFD-based approaches [29,30], where the source equations are highly simplified. From the CFD point of view, many turbulence models are investigated and improved either for the purpose of Atmospheric Boundary Layer (ABL) simulations or of wake dynamics. Direct Numerical Method (DNS) is not yet applicable for rotor simulations, instead, the Reynolds Averaged Navier-Stokes/Large Eddy Simulation/Detached Eddy Simulation (RANS/LES/DES) and the hybrid turbulence models are widely used [31,32,33,34]. For wake investigations, in particular wake meandering and wake stability analysis, the traditional CFD with body-fitted mesh is very expensive even for one- or two-rotor analysis. So far, the best configuration to achieve a good quality and simulation speed is the so-called AL/AD techniques [35,36,37] which still apply the above-mentioned turbulence models but relatively faster. Using the similar mesh topology as the AL/AD techniques, the Lattice Boltzmann Method (LBM) [38] and Immersed Boundary Method (IBM) [39,40] have an advantage of reducing the computational time, but still demand a large grid number for high Reynolds number flows.

The main purpose of the current paper is to develop a wind turbine rotor optimization tool that has a numerical accuracy of Class 4 and has a computational efficiency better than Class 2. More specifically, the current study will not only focus on single wind turbine rotor optimization but also for two or more aligned wind turbine rotors. In this scenario, the 3D CFD with RANS/AL/AD approaches are suitable methods to perform rotor aerodynamic analysis, but too computationally heavy to perform aerodynamic optimization. Benefiting from the axisymmetric characteristic of a HAWT, the 3D RANS/AD equations can be reduced into a 2D formulation. The axisymmetric boundary condition applied for a 2D simulation is a good assumption in particular for rotor optimization of ideal flow cases, such as no wind shear, no yawing effects. The novelty of such method is to achieve the same accuracy as using 3D CFD simulation with a significantly higher efficiency. Then, based on the multi-objective and multi-variable optimization algorithm, the design variables are the shape parameters comprising the chord, twist, and relative thickness of wind turbine rotor blade. Finally, the blade optimization is achieved based on an in-house developed flow solver combining with an optimization algorithm. A multi-objective optimization model is proposed with maximizing the overall power efficiency for multi-wind turbines. The non-dominated sorting genetic algorithm-II (NSGA-II) is used to handle the complex process. The novel part of the present study is to combine the CFD solver with an efficient optimization algorithm, such that the single or multiple wind turbine rotors are optimized at the same time, which finally obtain a global high aerodynamic efficiency.

The paper is organized as follows: Section 2 describes the governing equations of the aerodynamic models; Section 3 presents numerical validations of the current wake model; in Section 4, the aerodynamic methods are applied for rotor optimization. Conclusions are given in the final section.

2. Numerical Methods

The recent research development on wind turbine flow modeling has largely focused on CFD approaches. With the assumption of incompressible flows over rotor blades, the NS equations are solved together with RANS, LES, DES etc. The rotor blades can either be modeled as a solid body surface or with AD/AL/AS techniques. In addition to that, more numerical issues are involved such as inflow turbulence, and atmospheric boundary layer which lead to possible modifications of turbulence models.

2.1. Governing Equations

The incompressible RANS equations in 2D form reads

$$\frac{\partial u_i}{\partial x_i}=0 (i=1,2)$$

(1)

$$\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial x_j}=-\frac{1}{\rho }\frac{\partial P}{\partial x_i}+\frac{\partial }{\partial x_j}(\nu \frac{\partial u_i}{\partial x_j}-\overline{{u^{\prime }}_i{u^{\prime }}_j})+f_{ext}$$

(2)

$$-\rho \overline{{u^{\prime }}_i{u^{\prime }}_j}={\mu }_t\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)-\frac{2}{3}\rho k{\delta }_{ij}$$

(3)

where $x_i$ and $u_i$ are the displacement and velocity components in the inertial coordinate system. $t$ is the time, $\rho$ is the density of fluid, $P$ is the static pressure, $v$ is the kinematic viscosity coefficient, μ_t is the turbulent eddy viscosity δ_ij is the Kronecker delta function, k is the turbulent kinetic energy, $f_{ext}$ is the momentum source terms representing possible external forces.

2.2. Turbulence Model

Herein the turbulence model is discussed with the aim to model the nonlinear Reynolds stresses terms in Equations (2) and (3). Accurate predictions of turbulent kinetic energy and dissipation are the key factors in turbulent wake modelling. In the present work, the modified k-ω turbulence model of Menter [41] is used. Menter’s baseline model combines the k-ω model by Wilcox [42] in the inner region of the boundary layer with the standard k-ε model in the outer region and the free steam outside the boundary layer. The modified transport equations are written as follows:

$$\frac{\partial k}{\partial t}+\frac{\partial (u_{j}k)}{\partial x_j}=\frac{1}{\rho }{p}_k+\frac{1}{\rho }\frac{\partial }{\partial x_j}[(\mu +{\sigma }_k{\mu }_t)\frac{\partial k}{\partial x_j}]-{\beta }^{\ast }k\omega$$

(4)

$$\frac{\partial \omega }{\partial t}+\frac{\partial (u_j\omega )}{\partial x_j}=\frac{\gamma }{\rho {\upsilon }_t}{p}_k+\frac{1}{\rho }\frac{\partial }{\partial x_j}[(\mu +{\sigma }_{\omega }{\mu }_t)\frac{\partial \omega }{\partial x_j}]+2(1-F_1)\frac{\rho {\sigma }_{\omega 2}}{\rho \omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}-\beta {\omega }^2$$

(5)

where the original model constants for the inner region are:

$${\sigma }_{k1}=0.5, {\sigma }_{\omega 1}=0.5, {\beta }_1=0.075, {\beta }^{\ast }=0.09, \kappa =0.41, {\gamma }_1=0.55$$

and for the outer region

$${\sigma }_{k2}=1.0, {\sigma }_{\omega 2}=0.856, {\beta }_1=0.0828, {\beta }^{\ast }=0.09, \kappa =0.41, {\gamma }_2=0.44$$

(6)

The re-formulation of the original turbulence model contains two parts: (1) change of model constants; (2) added the source terms to maintain the turbulence level.

(1) The modified model constants are based on the work of Prospathopoulos et al. [43] in which most of the constants remain the same except for the following ones

$${\beta }_1=0.033, {\beta }^{\ast }=0.025, {\gamma }_1=0.37$$

The new constants were suggested according to measurements. With the decrease of dissipation effect, the kinetic energy towards the rotor is effectively increased, which improves the numerical accuracy for wind turbine wake flow modeling [43].

(2) For RANS solvers, the turbulence level is usually specified from the inlet of a computational domain. For the k-ω models, it is the initial k and ω values that are responsible for the turbulence level

$$k=\frac{3}{2}{\left(u_0{I}_0\right)}^2$$

(7)

$$\omega =\frac{k}{\mu }{\left(\frac{{\mu }_t}{\mu }\right)}^{-1}$$

(8)

where I₀ represents the ambient turbulence level and ${\mu }_t/\mu$ is the eddy viscosity ratio. A similar approach is implemented in some commercial CFD programs such as the ANSYS Fluent software. In another approach proposed by Richards and Hoxey [44], the boundary values of $u_0$, k and ω are basically a function of an aerodynamic roughness height corresponding to a specific site. This approach does not apply to the current work where an axisymmetric flow is assumed, and no wall boundary is specified. No matter what boundary values are given at inlet, the free decay of turbulence intensity cannot be avoided with the current set of RANS equations. In the free-stream flow, the mean velocity gradient is zero such that the turbulence quantities decay. Therefore, the production terms and the diffusion terms in Equations (4) and (5) can be neglected. The k and ω equations are reduced to

$$\frac{\partial (u_{j}k)}{\partial x_j}=-{\beta }^{\ast }k\omega$$

(9)

$$\frac{\partial (u_j\omega )}{\partial x_j}=-\beta {\omega }^2$$

(10)

Solving the partial differential equations for k and ω, the turbulence decay in the free-stream condition is obtained as a function of downstream distance

$$k=k_{in}{\left(1+\frac{{\omega }_{in}\beta x}{u_0}\right)}^{-{\beta }^{\ast }/\beta }$$

(11)

$$\omega ={\omega }_{in}{\left(1+\frac{{\omega }_{in}\beta x}{u_0}\right)}^{-1}$$

(12)

where k_in and ω_in are the initial values and x is the downstream distance. It is desired that these turbulence quantities specified at inflow boundary can be maintained without unphysical decay. In the work of Spalart and Rumsey [45], measurements of turbulent airfoil flows in wind tunnel were numerically reproduced with RANS computations. Even though the turbulence intensity is very low in the wind tunnel, the numerical simulations were not in perfect agreement with measurements. Inspired from the experimental and numerical work done in [45], the idea of maintaining turbulence quantities is extended to the current external wind turbine flows. The modified transport equations are

$$\frac{\partial k}{\partial t}+\frac{\partial (u_{j}k)}{\partial x_j}=\frac{1}{\rho }{p}_k+\frac{1}{\rho }\frac{\partial }{\partial x_j}[(\mu +{\sigma }_k{\mu }_t)\frac{\partial k}{\partial x_j}]-{\beta }^{\ast }k\omega +{\beta }^{\ast }{k}_{amb}{\omega }_{amb}$$

(13)

$$\frac{\partial \omega }{\partial t}+\frac{\partial (u_j\omega )}{\partial x_j}=\frac{\gamma }{\rho {\upsilon }_t}{p}_k+\frac{1}{\rho }\frac{\partial }{\partial x_j}[(\mu +{\sigma }_{\omega }{\mu }_t)\frac{\partial \omega }{\partial x_j}]+2(1-F_1)\frac{\rho {\sigma }_{\omega 2}}{\rho \omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}-\beta {\omega }^2+\beta {\omega }_{amb}^2$$

(14)

Different from the original model, the added source terms are ${\beta }^{\ast }{k}_{amb}{\omega }_{amb}$ and $\beta {\omega }_{amb}^2$ with the subscript “amb” indicates the ambient values [46]. It is quite straight forward to see that the destruction terms can be exactly cancelled if the inflow values of k and ω are set equal to the ambient values, which means in the free-stream flow case, turbulence quantities will remain unchanged throughout the computational domain. An example is given in Figure 1 that clearly shows the effect of using the added source terms. In the figure, the downstream distance x is non-dimensionalized to the rotor size D. In the free stream condition, the decay of turbulent kinetic energy using different approaches is compared. The original k-ω turbulence model agrees well with the theory (Equation (10)), where the modified model maintains the turbulent kinetic energy everywhere in the computational domain.

2.3. Actuator Disc

As the disc is assumed permeable, the mass conservation equation remains unchanged. In Equation (2), an external force f_AD is added that represents aerodynamic loads of wind turbine blades. To avoid the numerical singularity, each volume force needs to be smoothly re-distributed. A practical way of smearing the force is to use the 1D Gaussian function such that the actual forces in the disc is redistributed along one direction. In the current model, each force element is smeared locally in the blade normal direction with a distance d away from the disc using the convolution and the smeared force f′ is computed as [14,47]

$$f{}^{\prime }=f_{AD}\otimes {\eta }^{1D}$$

(15)

$${\eta }^{1D}(d)=1/\left(\epsilon \sqrt{\pi }\right)\exp \left[-{(d/\epsilon )}^2\right]$$

(16)

The parameter use in the smearing function is $\epsilon =3\Delta z$ where $\Delta z$ is the reference grid size in the axis direction. Based on the hypothesis of blade elements, the body forces of rotor blades can be achieved through airfoil data and the flow field information.

The actuator disc solves the entire axisymmetric flow field and the induced velocity is naturally included into the formulation i.e., the relative velocity $V_{rel}$ and flow angle $\varphi$ are determined from

$$\varphi ={\tan }^{-1}\left(\frac{V_0-W_z}{\Omega r+W_{\theta }}\right), V_{rel}^2={\left(V_0-W_z\right)}^2+{\left(\Omega r+W_{\theta }\right)}^2$$

(17)

As shown in Figure 2, the z-axis and θ-axis represent the rotor normal and rotational directions, respectively. The relative velocity V_rel seen by the blade element is a combination of axial velocity and rotational velocity, as shown in Equation (17). The lift force L and the drag force D are perpendicular and parallel to V_rel, respectively, as sketched in Figure 2. The local angle of attack is given by $\alpha =\varphi -\gamma$, where $\gamma$ is the local pitch and twist angle, the axial velocity $\left(V_0-W_z\right)$, where $V_0$ is the inflow wind speed and $W_z$ is the axial induced velocity, the tangential velocity $\left(\Omega r+W_{\theta }\right)$, where $\Omega$ is the rotational speed and $r$ is the radius of wind turbine and $W_{\theta }$ is the tangential induced velocity. Lift and drag forces per spanwise length are found from tabulated airfoil data and the velocity triangles shown in Figure 2.

$$L=\frac{1}{2}\rho V_{rel}^{2}cB{C}_l, D=\frac{1}{2}\rho V_{rel}^{2}cB{C}_d$$

(18)

Therefore, projecting the force in the axial and tangential directions to the rotor gives the aerodynamic force components

$$F_z=L\cos \varphi +D\sin \varphi , F_{\theta }=L\sin \varphi -D\cos \varphi$$

(19)

The forces are calculated in every time-step such that the momentum equation produces a new velocity field which updates the flow angle of attack as shown in Figure 2.

2.4. General Flow Solver

The RANS equations are solved by the EllipSys code [48,49], which is a general-purpose incompressible Navier-Stokes code based on a second-order multi-block finite volume method. The code solves the velocity-pressure coupled flow equations employing the SIMPLE/SIMPLEC/PISO (Semi-Implicit Method for Pressure-Linked Equations/SIMPLE Consistent/Pressure Implicit with Split Operator) methods and a multi-grid strategy. The momentum equations are first solved with a guessed pressure as a predictor. Next, the continuity equation is used as a constraint to obtain an equation for the pressure correction. In the predictor step, the momentum equations are solved by a second-order accurate backward differential scheme in time. For the spatial discretization a central difference scheme is employed for the diffusive terms and the QUICK (Quadratic Upstream Interpolation for Convective Kinematics) upwind scheme is utilized for the convective terms. In the corrector step, an improved Rhie-Chow interpolation technique is applied to suppress numerical oscillations from the velocity-pressure decoupling [50]. Furthermore, an improved SIMPLEC scheme developed for collocated grids is used to ensure that the solution does not depend on the value of relaxation parameters and the time-step [51]. The Navier-Stokes equations are discretized using a finite-volume scheme and all the information from the grid geometry is directly transferred into the discretized terms.

3. Results and Discussion

One of the major difficulties for engineering wind turbine wake models is to simulate the wake superposition effect. In this section, the 2D axisymmetric flows over two rotors are numerical validated against experimental data or LES results. The Danish Nibe wind turbines and Horns Rev wind turbines are selected for the numerical study.

3.1. Computational Setup

As sketched in Figure 3a, the two rotors are placed in tandem. The uniform inflow is applied at inlet, and the distance to the first rotor is about 10D (rotor diameter). The distance between downstream wind turbine and outlet is about 30D. The inlet and outlet width are chosen as 10D. The distance between the two rotors is not necessarily fixed which provides flexibility to be further optimized.

With a non-dimensional form, the uniform flow at the inlet and far field is Vo = 1. The change of wind turbine operational condition is through the Tip Speed Ratio (TSR). For the turbulence inflow conditions, the initial k and ω are computed through Equations (6) and (7) which are determined by the selected turbulence intensity. The mesh points are clustered around the rotor as seen in Figure 3b where only half the symmetrical plan is displayed. The mesh is block-structured with a total number of 14 blocks and 64 × 64 points in each block. The thick white lines indicate the two rotors whereas the mesh lines are clustered towards the rotors. With an OpenMPI (Message Passing Interface), the parallel simulation (on a Dell work station with Intel Xeon CPU E5 series) takes about 30 s for each simulation depending on the convergence rate. The computational time varies according to different TSRs as well as initial k and ω values. The simulation time using an AD method on a full 3D mesh will be more than hundred times slower than the 2D axisymmetric solver. On the other hand, the 2D CFD method is still slower than the classical steady BEM, but the 2D CFD method provides detailed flow field around the wind turbines.

3.2. Nibe Wind Turbine Case

To validate the aerodynamic model accuracy, the full-scale measurements [52] of the Nibe wind turbines are shown in Figure 4 where the two turbines are seperated with an axial distance of 5D. The Nibe wind turbines operate at a rated power of 630 KW, and the wind turbine has a rotor diameter of D = 40 m, and the hub height is H = 45 m. The simulation is performed at an inflow wind speed of 8.55 m/s and a turbulence intensity of 10%.

The comparisons between the measured and simulated wake velocity profiles are shown in Figure 5. The measurement data are avialable for a single wind turbine (Nibe B) at down-stream positions of 2.5D, 4.0D and 7.5D as depicted in Figure 4. It is worth noting that for the present single wake study, Nibe A is not operating during the measurement period. The detailed wake deficit is clearly observed from Figure 5a–c which ranges from the near wake to the far wake region. In terms of maximum wake deficit and turbulence intensity, the new model shown relatively better agreement with the measured data. At the near wake region, x = 2.5D, the new model does not produce a better velocity deficit than the original model. At x = 4D and 7.5D, the new model shows better results, which is also the normal spacing between wind turbines.

For the double wake case, the streamwise velocity contour plot is presented in Figure 6. The axial velocity is largely reduced as the flow first passes through the Nibe B wind turbine and than the Nibe A wind turbine. The comparisons between the measured and simulated wake profiles for the two wind turbines at down-stream positions of 2.5D, 4.0D and 7.5D are shown in Figure 7. In this case, both Nibe A and Nibe B operate at a wind speed of 8.55 m/s. At the downstream locations 2.5D and 4D, the wake profiles are expected to be similar as that shown in Figure 6. As expected, only tiny differences are seen from the simulated results, which are due to the flow feedback from the Nibe A wind turbine, but the effect is very small. On the other hand, the experimental data showed a wider spread at 2.5D and 4D which is caused by different measurement conditions, Figure 7c shows the more complicated double wake case where the wake flow of two wind turbines are superimposed. At x = 7.5D, wake data are both measured behind the Nibe B and Nibe A wind turbines. The velocity deficit is seen to be relatively small as previously observed in Figure 5c. However, when the Nibe A wind turbine is operating, the x = 7.5D position is only 2.5D apart from the Nibe A wind turbine and thus it is considered to be the near wake position of turbine Nibe A. As seen from Figure 7a,c, wake profiles are similar at positions x = 2.5D and x = 7.5D. However, at x = 7.5D, a smaller velocity deficit is observed, which is due to the wake mixing effect introduced by the Nibe B wind turbine. At x = 7.5D, it is also seen that the fluctuation of turbulence intensity is quite large which cannot be captured by any RANS model.

3.3. Horns Rev I Wind Farm Case

The Horns Rev I offshore wind farm is in the eastern North Sea near Denmark. Horns Rev wind farm consists of 80 Vestas V80 2 MW wind turbines within an area of about 20 km², as shown in Figure 8. The wind turbines are arranged in a regular array of 8 by 10 turbines which is very suitable to the study of multiple wake effects. The wind turbines are installed with an internal spacing along the main directions of 7D [53]. The wind turbines have a hub height of 70 m, a rotor diameter of D = 80 m, and a rated rotational speed of 18 rpm. The detailed pitch angle, rotational speed and blade geometry data are provided from the reference paper [54]. The 2 MW V80 wind turbine blade is composed of NACA63-series airfoils between the blade tip and the mid-span, and the FFA W3-series airfoils are used between the mid-span towards the blade inboard part. To validate the aerodynamic accuracy, the results from the present simulations and from the reference LES [26,54] simulations are compared. The inflow wind speed is 5.33 m/s, and turbulence intensity is 9.4% which are selected for fair comparison with some LES results.

The mesh topology used for this study is the same as in the previous simulations. The input airfoil data are now replaced to represent the Horns Rev wind turbines. The wind turbines under investigation are depicted in Figure 8 along the 270° directions. The wake profiles are shown in Figure 9 at downstream positions of 3D, 5D, 7D and 10D behind the first wind turbine where LES data are also available. General good agreements with LES are clearly seen from the results except very near wake. The wake expansion is well captured as well as the maximum velocity deficit. For the 7D and 10D cases, a slight over prediction of velocity deficit is seen from the current simulation. This indicates that the wake recovery at far downstream is slower using the present approach.

The CFD simulation of the full column was carried out by extending the computational domain with a total number of 40 blocks. The resulted stream-wise velocity contours are shown in Figure 10 where the flow around wind turbine cluster can be observed.

In Figure 11, information about the power and thrust of the present wind turbines are also shown. Good agreements are seen from the power and thrust curves, which proves the model accuracy in a wide range of operational conditions.

3.4. Rotor Optimizations

3.4.1. Optimization Setup

This section describes the optimization problem used in present rotor aerodynamic design. With the optimization objectives of AEP and ATP, and the design variables of chord length, twist angle, relative thickness and rotational speed, the NSGA-II optimization method [55] is used to optimize the wind turbine rotor. In the NSGA-II optimization algorithm, the binary coding, crossover mutation and gene mutation are built to avoid premature phenomena during the optimization process. To better understand the optimization process, a sketch of the optimization diagram is given in Figure 12. The design variables are the blade geometrical data as well as the rotational speed. The “Bspline” curve is applied to reduce the number of control points. The CFD solver is embedded in the optimization loop which is a time-consuming part of the whole design process. The number of necessary iterations is determined by the NSGA-II optimization routine. Each evolution contains 30–50 populations and each population takes about 60 min.

The design objective is to maximize the annual energy production (AEP) and minimize the rotor thrust force (ATP, annual thrust production), thus the design objectives are to minimize both 1/AEP and ATP

$$\mathrm{subject} \mathrm{to} X\in R^n,$$

(20)

$$x_k^L\le x_k\le x_k^U, k=1,\dots ,n$$

(21)

The vector $X$ contains a total of $n$ variables that are real numbers, $R$. The design variables, $X=[x_1,x_2,\dots ,x_n]$, are the control points that define the chord, twist angle and relative thickness as a function of blade span. The upper U and lower L notations, $x_k^L\le x_k\le x_k^U$, are the upper and lower limits of the control points. The 3rd-order Bsplines [56] are used to parameterize the chord, twist angle and relative thickness distributions, see Figure 13. The design is performed from a position at a radius of 0.337R to the tip of the blade. Thus, three control points for the chord, and two control points each for the twist angle and relative thickness distributions at the outboard of the blade are optimized. The 3rd-order Bspline curve is defined by

$$\begin{gathered} q_i(u)={\displaystyle \sum _{i=0}^n{N}_{i,p}(u)P_i} u_a<u<u_b, i=1,2,\cdots ,n \\ \begin{array}{c}{N}_{i,0}(u)=\{\begin{array}{c}1 if u_i\le u\le u_{i+1}\hfill \\ 0 otherwise\hfill \end{array}\hfill \\ N_{i,p}(u)=\frac{u-u_i}{u_{i+p}-u_i}{N}_{i,p-1}(u)+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}{N}_{i+1,p-1}(u)\hfill \\ u_0=u_1=\cdots =u_p=0\hfill \\ u_{i+p}=u_{i+p-1}+\left|p_i-p_{i-1}\right|, i=1,2,\cdots ,n-1\hfill \\ u_{n+p}=u_{n+p+1}=\cdots =u_{n+p+3}=1\hfill \end{array} \end{gathered}$$

(22)

where $q_i(u)$ is a vector-valued function of the independent variable, $u$. $\left\{P_i\right\}$ are the control points, and the $\left\{N_{i,p}(u)\right\}$ are the 3rd degree Bspline basis functions. The ith Bspline basis function of p-degree (order p+1). $u_i$ is the knot vector.

To calculate the AEP and ATP, it is necessary to combine the power curve with the probability density of wind. The Weibull distribution function defining the probability density can be written in the following form

$$f\left(V_i<V<V_{i+1}\right)=\exp \left(-{\left(\frac{V_i}{A}\right)}^k\right)-\exp \left(-{\left(\frac{V_{i+1}}{A}\right)}^k\right)$$

(23)

where $A$ is the scaling parameter, $k$ is the shape factor and $V_i$ is the wind speed distribution. If the wind turbine operates about 8760 h per year, its AEP and ATP can be evaluated as

$$AEP={\displaystyle \sum _{i=1}^{n-1}\frac{1}{2}}\left(P\left(V_{i+1}\right)+P\left(V_i\right)\right)\times f\left(V_i< V<V_{i+1}\right)\times 8760$$

(24)

$$ATP={\displaystyle \sum _{i=1}^{n-1}\frac{1}{2}}\left(T\left(V_{i+1}\right)+T\left(V_i\right)\right)\times f\left(V_i<V<V_{i+1}\right)\times 8760$$

(25)

where $P\left(V_i\right)$ and $T\left(V_i\right)$ are the power and thrust force at the wind speed of $V_i$.

3.4.2. Optimization Studies

The wind turbine rotor type of the Horns Rev I wind farm is chosen for the aerodynamic optimization. The selected wind turbines are located at the first and second rows (T1, T2 as depicted in Figure 8). It is worth noting that with the present method, it is possible to optimize the whole column of 8 wind turbines together. First, a larger number of wind turbines involve more design variables and thus a larger computational demand; second, it is the first two rows of wind turbines that have the major aerodynamic differences, as widely observed from measurements in several wind farms. According to the different optimization objectives, four optimization cases are carried out during the blade aerodynamic design.

Case 1: Baseline optimization for a single rotor. Using the optimization objectives of AEP and ATP, optimize the upstream wind turbine (the 1st row of rotors).

Case 2: Optimization of the 2nd row of rotors considering the wake effect from the 1st row of wind turbines. Using the optimization objectives of AEP and ATP, perform optimization for the 2nd wind turbine optimization.

Case 3: Further optimization of the rotational speed of the 1st row of wind turbines. The optimization objectives are the AEP and ATP for the two rows of wind turbines. The blade shape of the two types of rotors is obtained from Case 1 and Case 2, respectively, but the rotational speed of the 1st row of rotors will be optimized further.

Case 4: Based on Case 3, further optimization is performed for the two rows of wind turbines. Using the optimization objectives of both AEP and ATP, the two rows of wind turbines are optimized together.

To summarize these optimization cases: Case 1 follows a classical blade optimization routine; Case 2 only focuses on a downwind rotor, its blade shape is optimized; Case 3 optimizes the rotational speed of the first row of rotors to seek the global maximum energy production for two rows of rotors; Case 4 optimizes both the rotational speed of the first row of rotors and the blade geometry of the second row of rotors. By performing the 4 steps analysis, it expected that after the Case 4 study a more realistic optimization can be achieved for such a typical multiple wind turbine system.

The optimization Pareto fronts are shown in Figure 14 with the resulted Pareto front that measures how the solutions are distributed. The four subfigures (a–d) represents the above mentioned four cases. In the figures, the areas marked with A, B and C indicate the regions of feasible optimal solutions. The trend of area A suggests the solutions with high AEP, whereas low thrust is more focused in area C. As a multi-objective method, balance is then needed to weight between the high AEP and low ATP. In Figure 14a, Case 1 represents the standard optimization of a single rotor (the 1st row of rotors) without any wake effect. The high AEP and high ATP region is marked as area A, whereas area B finds good balance between low ATP and high AEP. In Figure 14b, the second row of wind turbine is optimized. Due to the wake effect from the first row of wind turbine, both the ATP and AEP are reduced comparing with Case 1 results. In Figure 14c, a high AEP area can hardly be achieved. The reason is that in Case 3 the AEP from the upstream rotor is reduced to increase the energy output for the downwind rotor. Such that, it is necessary to further optimize the downwind rotor. In Figure 14d, the optimization is continued such that the two rows of wind turbines are fully optimized with satisfactory results obtained.

Some key numbers are shown in

Table 2. Optimization results of Case 1 and Case 2.

Case	Case 1 (1st Rotor)		Case 2 (2nd Rotor)
Objective	AEP	ATP	AEP	ATP
Reference	4.2665	0.8058	2.6559	0.6410
Optimized	4.3153	0.8058	2.7031	0.6432
Relative	1.14%	0	1.78%	0.34%

Case 1: Single wind turbine, consider both AEP&ATP; Case 2: Two wind turbines, consider both AEP and ATP of Downstream wind turbine.

and

Table 3. Optimization results of Case 3 and Case 4.

Case		Case 3 (1st Rotor)		Case 4 (Two Rotors)
Objective		AEP	ATP	AEP	ATP
Reference	1st	4.2653	0.8058	4.2653	0.8058
	2rd	2.6559	0.6410	2.6559	0.6410
	Overall	6.9212	1.4468	6.9212	1.4468
Optimized	1st	4.2307	0.7919	4.1177	0.7522
	2rd	2.6875	0.6304	2.9139	0.6432
	Overall	6.9182	1.4224	7.0317	1.3955
Relative	1st	−0.81%	−1.72%	−3.46%	−6.65%
	2nd	1.2%	−1.65%	9.71%	0.34%
	Overall	−0.04%	−1.69%	1.6%	−3.55%

Case 3: Upstream WT rotation speed; consider AEP of overall WT, and ATP of Downstream WT; Case4: Upstream WT rotation speed and blade geometry; consider AEP of overall WT, and ATP of Downstream WT.

after each optimization. For the Case 1 design, the upstream rotor is first optimized. As shown in Figure 14a, a trade-off between thrust and power can be found in area B. The final selected data is listed in Table 2 where the same ATP is maintained but the AEP is increased by 1.14%. For the Case 2 design, the downwind rotor is optimized. After optimization, the rotor in the wake has an increased power production of 1.78%, and there is also a small increase of thrust of 0.34%. Using the optimized solutions from Case 1 and Case 2, further investigation is carried out in Case 3. The rotational speed of the upstream wind turbine is considered to be a new design variable that may affect the total energy product. In Table 3, by comparing the AEP of each wind turbine with its reference value, it is seen that the 1st wind turbine has a negative change in AEP of −0.81% and the second wind turbine has a positive change in AEP of 1.2%. The overall gain of AEP is −0.04% because of the contribution of higher wind resource from the 1st wind turbine. The total thrust change of the two wind turbines is about 1.69%. It is seen that loads can be decreased with a certain amount by decreasing the front rotor speed, but it is not possible to gain more power from Case 3 study. A study in Case 4 is therefore performed to seek more energy from the 2nd rotor. The individual changes from each wind turbine and the overall changes are listed in Table 3. The reduction of the 1st rotor speed leads to a change of its AEP with −3.46%. The 2nd rotor is further optimized with a relatively large amount of AEP increase up to 9.71%. The total amount of AEP in the two-wind-turbine system is increased with 1.6%. The thrust is largely reduced for the 1st wind turbine under the new design condition at a relatively low TSR.

The final aerodynamic layouts of the optimized blades are gathered in Figure 15. The chord, twist and relative thickness distributions are given in Figure 15a–c, respectively. Based on the multi-objective optimizations, the chord length of the upstream wind turbine (Case 1) and downwind wind turbine (Case 2) deviates from the baseline chord distribution. The twist distribution of the two optimized blades is rather similar. The optimized rotor speed is shown in Figure 15d. With limited lower and upper bounds, the Case 4 optimization suggests a much lower rotational speed as compared with Case 3. As a result, the upwind rotor lost a certain amount of power which is gained by the downwind wind turbine. It is worth noting that the present study is only limited at aerodynamic optimization. Since the 1st wind turbine also has a decrease of thrust of 6.65%, there is a big potential to further reduce the structure weight.

4. Conclusions

This paper presented detailed numerical methodologies and convincing results for wind turbine rotor aerodynamic optimizations. To simulate the aerodynamic forces on the wind turbine rotor, a 2D actuator disc method is implemented in the incompressible flow solver favored by the axisymmetric boundary condition defined at the center line. The combined flow solver is more accurate and powerful than existing engineering models such as blade element momentum theory. On the other hand, it is order of magnitude faster than the 3D flow solver. However, the 2D and 3D solver implemented with same AD method should give same results. In addition to the AD flow simulations, two turbulent flow cases were investigated to proof the advantage of the modified k-ω turbulence model. With the fast CFD solver, the rotor optimization was carried out for single wind turbine with uniform flow and two wind turbines combined with wake interaction. The multi-objective optimization method NSGA-II is applied in the design loop and the two design objectives of high AEP and low ATP are required. The four optimization cases aimed at single rotor optimization without wake effect, single rotor optimization with wake effect, multi-wind turbine optimization with constant rotor speed and multi-wind turbine optimization with variable rotor speed. In general, single rotor optimizations are faster and they provide clear optimization trend. For the multi-wind turbine optimization cases, it was found that high aerodynamic performance is achieved when the upstream and downstream wind turbine both change their blade shape and rotational speed. These results indicate that traditional rotor optimization methods might not be sufficient due to the complex wake effect.

On a single PC with a developed parallel computation tool, the estimated time is 30 h for optimizing a single rotor, and for the coupled two-rotor optimizations, the simulation time is almost doubled. Future work will be carried out for more complicated design cases which may include the whole column of wind turbines such as the Horns Rev wind farm case. In such an optimization, the wind turbine spacing will be included as an important design parameter, to this end, the method can be applied for simple layout optimization.