The iteration method for tower height matching in wind farm design

J. Wind Eng. Ind. Aerodyn. 132 (2014) 37–48

Contents lists available at ScienceDirect

Journal of Wind Engineering and Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia

The iteration method for tower height matching in wind farm design K. Chen, M.X. Song, X. Zhang n Key Laboratory for Thermal Science and Power Engineering of Ministry of Education, Department of Engineering Mechanics, Tsinghua University, Beijing 100084, PR China

art ic l e i nf o

a b s t r a c t

Article history: Received 18 December 2013 Received in revised form 16 June 2014 Accepted 17 June 2014

This paper studies the tower height matching problem in wind turbine positioning optimization. Various models are introduced, including the power law wind speed model with height in the wind farm, the linear wake flow model for flat terrain, the particle wake flow model for complex terrain and the power curve model with power control mechanisms. The greedy algorithm is employed to solve the wind turbine positioning optimization at a specified tower height. The optimization objective is to maximize the Turbine-Site Matching Index (TSMI), which includes both the production and the cost of wind farm. Assuming that the optimized layout for each tower height is the same, an iteration method is developed to obtain the approximated optimal height. The convergence of the proposed iteration method is discussed through the mathematical analysis. The proposed iteration method is validated through the numerical cases over both flat terrain and complex terrain. The results indicate that the proposed method can obtain better optimized height in shorter computational time than previous studies. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Tower height matching Wind turbine positioning optimization Greedy algorithm Iteration method Wind farm

1. Introduction Nowadays, energy shortage problem has become one of the serious problems all over the world. In order to mitigate this problem, people start to pay attention to all kinds of renewable energy resources, including solar energy, geothermal, biomass, tide and wind energy. Among these renewable energies, wind energy is an important alternative energy due to the advantages of clean and rich resources (Chen and Zhang, 2007). In China, there are a large amount of wind energy resources. The potential wind power at 10 m height on the land and on the sea are 253 GW and 750 GW, respectively (Tong and Dong, 2012). In the past decade years, wind energy is developing rapidly in China. The total installed capacity of wind turbines reached 75.3 GW up to 2012, which was the largest in the world (Song, 2012). Wind energy is extracted by wind turbine in wind farm. The power output of the wind turbine increases as the wind speed increases. Meanwhile, the wind turbine will generate a wake region downstream due to the extraction of the wind power and the disturbance of the wind rotor. In the wake region, the wind speed is reduced and the turbulence is increased. Therefore, the wind turbine positions should be designed to reduce the wake effect and increase the total power output of the wind farm.

n

Corresponding author. E-mail address: [email protected] (X. Zhang).

http://dx.doi.org/10.1016/j.jweia.2014.06.017 0167-6105/& 2014 Elsevier Ltd. All rights reserved.

Much investigation has been done on wind turbine positioning optimization (WTPO). Researchers introduced many optimization algorithms to solve the problem. Genetic algorithm was the first algorithm introduced to solve WTPO by Mosetti et al. (1994). This algorithm simulates the biological evolution process. Bases on the population, the algorithm searches the optimized solution through the selection, crossover and mutation operators. In Mosetti's study, binary coding method was used combined with the linear wake model and 3-order power curve model. The target was to maximize the production per unit cost. The effectiveness of genetic algorithm on WTPO was validated by three numerical cases. Based on Mosetti's study, others have used larger population and more generations (Grady et al., 2005), and more realistic models (Mora et al., 2007; Kusiak and Song, 2010) to improve the optimized results. Wan et al. (2009) used real coding genetic algorithm to optimize the wind turbine positions with the target of maximizing the total power output with the number of wind turbines fixed, obtaining better results than the ones by binary coding genetic algorithm. Another type of optimization algorithms used in WTPO is greedy algorithm. Greedy algorithm is based on a single turbine layout and the turbines are placed in the positions one by one that make the objective value maximum in each step. Compared to genetic algorithm, greedy algorithm requires less computation and the optimized result does not have randomness. Ozturk and Norman (2004) combined greedy algorithm with the adding, removing and moving operators to optimize the wind turbine positioning problem. Based on the submodular property in the optimization problem with linear wake model, Zhang et al. (2011)

38

K. Chen et al. / J. Wind Eng. Ind. Aerodyn. 132 (2014) 37–48

used the lazy algorithm to reduce the computational time of the algorithm. Song et al. (2012) developed a particle wake model that can be used over complex terrain. Based on this particle wake model, the greedy algorithm with repeated adjustment was developed to optimize wind turbine positioning problem over complex terrain (Song et al., 2013). Besides, some other optimization methods were also introduced in WTPO, including simulated annealing method (Rivas et al., 2009), Monte Carlo method (Marmidis et al., 2008) and particle swarm optimization method (Wan et al., 2010). The wind speed of the wind farm and the cost of the wind turbines will increase with the tower height. Therefore, the tower height of the turbines should match the potential site to achieve maximum power output per unit cost. In the literature, the tower height matching problem has been considered to further improve the turbine layout based on WTPO (Chen et al., 2013a). The Turbine-Site Matching Index (TSMI) was introduced as the objective function, including the production and the cost of the turbine layout. The greedy algorithm with repeat adjustment was introduced to solve WTPO. The optimal height of wind turbine can be obtained through the enumeration method. That is, apply WTPO at each optional tower height. Then the height with the maximum objective value is the optimal tower height. However, it requires a large amount of computational time, especially when the turbine has a large range of optional tower height. In previous study, the fitting method was developed to obtain the optimized height in less computation. When using the fitting method, the normalized power output (L) is defined. The optional height range of the wind turbine is divided into several parts with the same interval and the splitting points are obtained. Then apply WTPO at each splitting point and reproduce the whole L curve using the L values at these points through polynomial fitting. Finally, calculate the extreme points of TSMI using the L curve and obtain the optimized point. The extreme points with the maximum objective value is the optimized height. Three numerical cases were used to test the performance of the fitting method. The results indicated that the fitting method can obtain the approximated optimal height in fewer less times of applying WTPO than the enumeration method (Chen et al., 2013a). However, it needs at least 5-order fitting to obtain the optimized results with the error less than 5% for the multi-direction wind situations. In this paper, the tower height matching for WTPO is studied. Assuming that the optimized layout for each optional tower height is the same, an iteration method is developed to obtain the optimized tower height. The convergence of the iteration method is discussed through the mathematical analysis. The effectiveness of the proposed method is validated by the numerical cases over both flat terrain and complex terrain.

u0

The optimized results by three methods are compared for each case, including


Enumeration method: Apply WTPO at each optional height and





take the height with maximum objective value as the optimal one. The result of this method is treated as the optimal result of the tower height matching problem. Fitting method: Obtain the optimized height through fitting the L curve, developed in previous study (Chen et al., 2013a). Iteration method: Assuming that the optimized layout for each optional tower height is the same, the optimized tower height is obtained by an iteration process, developed in present study.

The remainder of the paper is organized as follows. Section 2 presents the models introduced in WTPO. Section 3 introduces the optimization methodology. Section 4 presents the iteration method for tower height matching problem. Section 5 discusses the numerical results of the test cases. Section 6 presents the conclusions.

2. Models 2.1. Linear wake model In this paper, the linear wake model used in the study of Mosetti et al. (1994) is employed to calculate the wind turbine wake effect for the wind farm on flat terrain. The wake model is considered to be a conical area, as shown in Fig. 1. The velocity inside the wake region is calculated by the following algebraic expression: 2 3 6 u ¼ u 0 41 


1þα

x r1

7 2 5

ð1Þ

where u0 is the local wind speed without placing the turbine, x is the distance downstream the turbine rotor, r1 is the downstream rotor radius, a is the axial induction factor and α is the entrainment constant, which are expressed as follows: pffiffiffiffiffiffiffiffiffiffiffiffiffi 1  1  CT a¼ ð2Þ 2 rffiffiffiffiffiffiffiffiffiffiffiffiffi 1a ð3Þ r1 ¼ r 1  2a

α¼

0:5 lnðh=z0 Þ

ð4Þ

where C T is the trust coefficient, r is the radius of the wind rotor, h is the tower height of the wind turbine, and z0 is the surface roughness. The size of the wake region is described by the wake influenced radius R, which is the radius of the wake region at a specified section in the crosswind direction, expressed as

u0

R ¼ αx þ r u

r

2a

R=αx+r

x Fig. 1. Schematic of linear wake model (Chen et al., 2013a).

ð5Þ

Considering multiple wake interference effect, the velocity of the ith turbine is calculated by Gonzalez et al. (2010) vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 1 " u  2 # u N A u ij ij t A ð6Þ ui ¼ u0i @1  ∑ 1 2 u0j j ¼ 1 π ri where u0i and u0j are the local velocities at the ith and the jth turbines' positions without placing the turbines. They are equal to the inlet speed of wind farm over flat terrain. uij is the wind speed at the wind rotor of ithe turbine in the wake region of the jth turbine, N is the number of wind turbines, ri is the rotor radius of

K. Chen et al. / J. Wind Eng. Ind. Aerodyn. 132 (2014) 37–48

the ith turbine, and Aij is the rotor area inside the jth turbine wake, calculated by 8 2 πr ; dij rr wj r i > > > i < 1 2 1 2 Aij ¼ r ðθ1  sin θ1 Þ þ r wj ðθ2  sin θ2 Þ; r wj  r i o dij o r wj þ r i > 2 i 2 > > : 0; dij Zr wj þr i ð7Þ where 2

θ1 9 ∠AOj B ¼ 2 arccos

r 2wj þ dij  r 2i

ð8Þ

2r wj dij 2

θ2 9 ∠AOi B ¼ 2 arccos

r 2i þ dij  r 2wj

ð9Þ

2r i dij

where r wj is the wake influenced radius of the jth turbine by Eq. (5), dij is the crosswind distance between the ith and the jth wind turbines, and θ1 and θ2 are shown in Fig. 2. As the near wake flow is very complicated due to the disturbance of wind rotor, the linear wake model cannot reproduce the near wake flow. This model is only applicable to the far wake flow, which is usually four times of the wind rotor diameter away from the turbine hub in downstream. On the other hand, this wake model does not consider the effect of the increase in turbulence intensity on the power output. However, the linear wake model is a good approximation to wake flow over flat terrain. Furthermore, the wake flow is calculated using the algebraic expression, which can obtain the wake velocity in short time. Therefore, the linear wake model is commonly used combining with the optimization algorithms to solve the turbine positioning problem over flat terrain. 2.2. Particle wake model For complex terrain, the shape of the wake flow will change, so the linear wake model cannot be applied. Song et al. (2012) proposed the particle wake model that can be applied over complex terrain. Previous study has shown that particle wake flow can be combined with the optimization algorithm to optimize the wind turbine positioning problem (Song et al., 2013). In the particle wake model, the turbine wake flow is calculated by particle simulation and the velocity deficit in the wake region is represented by the particle concentration. The steps of the particle wake model are shown as follows (Song et al., 2012): 1. The velocity of the field should be obtained through the Computational Fluid Dynamics (CFD) method. The obtained flow field is named the pre-calculated flow field. 2.

39

The particles are generated within the area of the wind rotor and moved based on the pre-calculated flow field. 3. In each time step of particle simulation, the convective effect, the diffusive effect and the attenuation effect are considered. The convective displacement and the diffusive displacement satisfying Gaussian distribution are added to each particle. The total displacement is shown as

Δx ¼ ð1 þ σ rÞuΔt

ð10Þ

where u is the local velocity interpolated through the precalculated flow field and Δt is the span of a time step, σ is the diffusive coefficient of the particle wake model, and r is a Gaussian distributed random number. During the particle simulation, the particles are disappeared according to a attenuation coefficient γ. 4. A cube with the side length as same as the wind rotor diameter is used to calculate the particle concentration. The particle concentration of each time step is summed and averaged. After hundreds of time steps of simulation, the particle distribution becomes stable. The averaged number of particles inside the cube during the particle simulation is normalized to the relative particle concentration, denoted by c. 5. The wake influenced velocity is calculated by the following transformation expression: u0 ¼ uð1  βcÞ

ð11Þ

0

where u is the wake influenced velocity, u is the velocity of the pre-calculated flow field, β is the transformation constant and c is the relative particle concentration. In the particle wake model, σ, β and γ are three parameters to be determined, which depend on the characteristics of the wind turbine. The particle model still does not consider the effect of the increase in turbulence intensity on the power output. However, this model contains the influence of the terrain through the particle simulation based on pre-calculated flow field. Therefore, it can calculate the wake flow over complex terrain. 2.3. Total power output of wind farm The power output of a wind turbine is expressed by the power curve model, which gives the power output for each wind speed. The turbine power curve is characterized by the characteristic speeds of the wind turbine, expressed as (Albadi and El-Saadany, 2010) 8 u o uc or u Zuf > < 0; P ðuÞ; u P e ðuÞ ¼ P r  ð12Þ asc c ru o ur > : 1; ur r u o uf

A rwj

rwj

ri

2

dij

oj

1

oj

dij

ri

oi

oi B

Fig. 2. Schematic of wake area (Chen et al., 2013a): (a) dij r r wj  r i and (b) r wj  r i o dij o r wj þ r i .

40

K. Chen et al. / J. Wind Eng. Ind. Aerodyn. 132 (2014) 37–48

where uc , ur and uf are the cut-in speed, the rated speed and the cut-out speed of the turbine, respectively. P e is the actual power output of the wind turbine, P r is the rated power output and P asc is P e as the percentage of P r . In the present study, P asc is calculated by the cubic model, expressed as (Albadi and El-Saadany, 2010) u3 P asc ðuÞ ¼ 3 ur

ð13Þ

The 2-parameter Weibull distribution is commonly used to model the wind speed characteristics, shown as (Burton et al., 2001) k
uðk  1Þ  ðu=cÞk e ð14Þ f ðuÞ ¼ c c where k is the shape parameter and c is the scale parameter. After discretized, the discrete wind speeds ui and the corresponding pi are obtained. The total power output is the sum of the outputs of all the turbines, expressed as " !# M

M

P tot ¼ ∑ ½pi P layout ðui Þ ¼ ∑ i¼1

i¼1

N

pi

∑ P e ðuij Þ

ð15Þ

j¼1

where P tot is the total power output of the wind farm, M is the number of the wind cases, pi is the probability of the ith wind case and P layout ðui Þ is the power output of the turbine layout with the incoming velocity ui. N is the number of the wind turbines and P e is the power output calculated by Eq. (12). Normally, wind speed will increase with the height from the ground which is commonly modeled by power law or logarithmic law (Burton et al., 2001). In the present study, the power law model is used, shown as  αu h u ¼ uref ð16Þ href where uref is the wind speed at the reference height href , named the reference wind speed. αu is the wind shear coefficient. Eq. (16) is substituted into Eq. (15), shown as " !# M

P tot ¼ ∑

N

i¼1

pi

∑ P e ðHðhÞðuij Þref Þ

ð17Þ

j¼1

The wind speed increases with the tower height, increasing the total power output of the wind farm. 2.4. Turbine-site matching index When constructing wind farm, the profit is always desirable for optimality, through increasing the power output and decreasing the cost. In this paper, the Turbine-Site Matching Index (TSMI) is introduced as the objective of the tower height matching problem, the same as the previous study (Chen et al., 2013a). TSMI represents the production per unit cost of the wind farm, expressed as TSMI ¼

P tot =P r ICC

ð18Þ

where the numerator of TSMI is the total power output of wind farm normalized by P r . The denominator is the normalized cost ICC, modeled by Albadi and El-Saadany (2010) ICCðhÞ ¼

ICCðhÞ ¼ aC  h þ bC ICC80 m

A potential factor (PF) is defined to describe the total height effect on TSMI, shown as PF ¼

TSMImax  TSMImin TSMImin

ð20Þ

where TSMImax and TSMImin are the maximum TSMI and the minimum TSMI, respectively. PF describes the maximum potential improvement of TSMI after optimization. The larger PF is, the larger the possibility of potential improvement is.

3. Wind turbine positioning optimization 3.1. Greedy algorithm with repeated adjustment For a specified tower height, tower height matching problem becomes a traditional WTPO problem, which can be solved by the optimization methods in the literature. In order to compare with the results of previous study, the greedy algorithm with repeated adjustment (Chen et al., 2013a) is introduced to solve the WTPO problem. The algorithm consists of two stages: locating stage and relocating stage. Locating stage starts from the empty wind farm. The steps of locating stage are listed below. 1. The wind farm is meshed by square grids and all the grids are numbered. 2. Add a turbine at the grid with the maximum evaluation value. 3. Try each empty grid with a wind turbine and calculate the evaluation value. 4. Add a turbine at the grid with the maximum evaluation value and record the order of the turbines adding in the wind farm. When there are some grids having the same maximum evaluation value, it will choose the grids with small number to place the wind turbine. 5. If the specified number of turbines is reached, locating stage is completed. Otherwise, return to step 4 and continue the procedure. When locating stage is finished, a turbine layout with the specified number of wind turbines is obtained. However, the strategy of turbine placement in locating stage cannot guarantee obtaining the optimal solution, so relocating stage is used for further optimization of the turbine layout. The steps of relocating stage are shown as follows: 1. Remove one turbine according to the order of the turbines in the adding process. 2. Try each empty grid with a turbine and calculate the evaluation value. 3. Add a turbine back to the layout at the grid with the maximum evaluation value. 4. A cycle is defined when all the turbines go through a relocating procedure. If the layout does not change during the whole cycle, relocating stage is completed. Otherwise, return to Step 1 and continue the procedure. When relocating stage is finished, the final optimized layout is obtained. In the algorithm, the evaluation function is chosen as the objective function (TSMI), which is desirable for optimum.

ð19Þ

where ICC80 m is the ICC value at the height of 80 m, which is the characteristic ICC value. In present study, aC and bC are chosen as 1:1875  10  3 and 0.905 (Albadi and El-Saadany, 2010). As the tower height increases, P tot and ICC both increase. There exists a optimal height for the largest TSMI value.

3.2. Penalization of the objective In the wind farm, the wind turbines should be placed far enough in case of damaging others when falling down. This is the safe distance condition. In order to guarantee this safe distance condition satisfied, the objective value is penalized and the

K. Chen et al. / J. Wind Eng. Ind. Aerodyn. 132 (2014) 37–48

5. Numerical study

evaluation function is shown as

E¼

8 > > > < 0; > > > : TSMI;

  Dij oλ min Rij i;j   Dij min Zλ Rij i;j

41

ð21Þ

where λ is the adjust coefficient, Rij is the sum of the heights of two turbines, including the hub height and the rotor radius, and Dij is the horizontal distance between two turbines. Through penalization, the wind turbines tend to locate at the positions satisfying the distance condition ðDij Z λRij Þ. The distances among the turbines can be adjusted through changing the value of λ. λ should be greater than 1 to guarantee the safe distance condition.

In this section, the proposed iteration method is validated by numerical cases. Section 5.1 considers three typical cases over flat terrain in previous study, including the single-directional wind case with a single speed, the multi-directional wind case with uniform speeds and the multi-directional wind case with variable speeds. The proposed method is compared to the fitting methods and the enumeration method through these three cases. In Section 5.2, Weibull distribution is introduced to model the wind speed characteristic and a cases with Weibull distribution over flat terrain is considered. Section 5.3 considers a multidirectional wind case over complex terrain to further test the iteration method. 5.1. Three typical cases

4. Iteration method for tower height matching When the power curve of the wind turbine satisfies the cubic model for all wind speeds, previous study has shown that the WTPO is the same optimization problem at each optional height when ignoring the height dependency in wake effect and considering the distances among turbines constant (Chen et al., 2013a). Therefore, the optimal wind turbine layout for each tower height is the same. In this situation, through applying one WTPO at a specified height, the optimized layout for this height can be obtained. Translate the optimized layout to other tower heights and evaluate the TSMI values to reproduce the TSMI curve. Then the height with the maximum TSMI value is the optimal height for the tower height matching problem. For the real situation, the optimal wind turbine layout for each height is different due to the power control mechanism and the height effect of the models. However, as an assumption, the optimal layout for each height is still considered to be the same. Then an iteration method for tower height matching can be developed. The detail steps of the method are shown below. 1. Choose an initial standard height, denoted by h0. 2. Apply WTPO at h0, obtaining the relevant optimized turbine layout. 3. Translate the optimized layout to other heights and reproduce the TSMI curve through evaluate the layouts. 4. Take the height with the maximum TSMI value as the optimized height, denoted by hite . If hite equals to h0 or the maximum number of iteration step is reached, stop the procedure. Otherwise, let h0 equal to hite , return to Step 2 and continue the procedure. 5. When the iteration procedure is finished, the optimized height for each iteration step is obtained. Choose the height with the largest TSMI value among these optimized heights as the final optimized height.

As the iteration method cannot guarantee obtaining a convergence optimized height, a maximum number of iteration steps is selected. When the maximum number is reached, the iteration procedure is stopped. The approximated optimal height in the tower height matching problem can be obtained through the iteration method. Generally, the optimal layout depends on the optimization method. Thus, the optimal layouts for similar tower heights may be quite different, but the objective values of the optimal layout for similar tower heights are very close to each other. This is the characteristic that makes the iteration method effective. The convergence of the proposed iteration method is discussed in Appendix A.

In this section, three typical cases over flat terrain are used to test the proposed iteration method and compare to the previous study (Chen et al., 2013a). The reference wind speed distributions of the cases are shown as follows.


Case 1: Single wind direction with a single wind speed.
Case 2: Multi-directional (16 directions with intervals of 22.51) wind with uniform speeds.


Case 3: Multi-directional (16 directions with intervals of 22.51) wind with variable speeds. The parameters of the numerical study are listed in Table 1. The calculation domain is chosen as a square area with the size of 2000 m  2000 m. The ground roughness of the wind farm is 0.3 m and the wind shear coefficient is 1/7 (Chen et al., 2013a). The wind from west to east is defined as 01 and the one from south to north is defined as 901. The linear wake model is used to calculate the wind turbine wake flow. The number of wind turbines to be placed in the domain is 30 and the adjust coefficient λ in Eq. (21) is set as 1.05. Table 2 lists the properties of the wind turbine used. The optional tower height of the turbine is in the range ½40; 140. The power curve model is expressed as Eq. (12). Greedy algorithm with repeated adjustment is used to solve the wind turbine positioning problem with a specified tower height. The grids for Table 1 The parameters of the numerical study (Chen et al., 2013a). Parameter

Value

Domain Ground roughness Wind shear coefficient Grids for optimization Number of wind turbines Adjust coefficient

2000 m  2000 m 0.3 m 1/7 30  30 30 1.05

Table 2 Wind turbine properties (Chen et al., 2013a). Property

Value

Rotor diameter D (m) Trust coefficient C T Tower height (m) Cut-in speed (m/s) Cut-out speed (m/s) Rated speed (m/s) Rated power (kW)

40 0.88 40–140 2 18 12.8 630

42

K. Chen et al. / J. Wind Eng. Ind. Aerodyn. 132 (2014) 37–48

optimization is set as 30  30. The initial tower height of the iteration method is chosen as 40 m and the maximum iteration steps is 4. Note that the distances among the wind turbines may be

less than 160 m (four times of the wind rotor diameter) when the tower height is less than 57 m for λ being 1.05. In this situation, the linear wake model cannot be applied. Therefore, when the

Fig. 3. TSMI curves of Case 1: (a) uref ¼ 9 m=s; (b) uref ¼ 9:4 m=s; (c) uref ¼ 9:6 m=s; and (d) uref ¼ 10 m=s. Table 3 Comparison of results for Case 1. uref (m/s)

9 9.4 9.6 10

PF (%)

44.9 32.4 26.1 14.8

Averaged error N opt

Enumeration

Fitting method (m)

Iteration method (m)

hopt (m)

hfit;3

hfit;4

hite;1

hite;2

hite;3

hite;4

128 98 88 63

131 100 75 61

135 98 83 64

118 101 89 66

137 101 89 66

128 – – –

128 – – –

– 101

1.51% 5

1.19% 6

2.22% 2

1.13% 3

0.72% 4

0.72% 5

2000

2000

1500

1500

1000

1000

500

500

0

0

500

1000

1500

2000

0

0

500

1000

1500

Fig. 4. Optimized wind turbine layouts of Case 1 ðuref ¼ 10 m=sÞ: (a) h ¼ 63 m and (b) h ¼66 m.

2000

K. Chen et al. / J. Wind Eng. Ind. Aerodyn. 132 (2014) 37–48

43

comparison. The relative error is defined to evaluate the optimized result, shown as

tower height is less than 57 m, the distances among wind turbines are set to be at least 160 m to guarantee the validity of the linear wake model. The tower height matching can be solved by the enumeration method, which applies optimizations in the interval of 1 m for each height. The height with the largest objective value is the optimal height. The optimized results by the fitting method in previous study (Chen et al., 2013a) are also presented for

Errorðf i Þ ¼

jf optimized  f optimal j maxi ðf i Þ  mini ðf i Þ

ð22Þ

where f optimized is the optimized objective value using the iteration method or the fitting method and f optimal is the optimal objective value using the enumeration method. The denominator of the relative error is the discrepancy of the best objective value and the worst objective value. The relative error ranges from 0 to 1. When the relative error equals to 0, it obtains the best optimized result. When the relative error equals to 1, it obtains the worst optimized result. The smaller the relative error is, the better the optimized result is.

5.1.1. Case 1 In Case 1, the single wind direction with one single speed is considered and various values of uref are studied. Fig. 3 shows the TSMI curve with the height for each uref , which contains the results of the enumeration method and the proposed iteration method. In Fig. 3, hn represents that the standard height of the nth iteration step. The results show that the TSMI curves of various standard heights are close to the accurate TSMI curve. Especially around hn, the errors of the curves are small. As the iteration steps increases, hn gets closer and closer to the optimal height of the accurate TSMI curve. Table 3 shows the optimized results of the three methods for Case 1. hopt is the optimal tower height by the enumeration method, which is treated as the optimal height. hfit;n represents the optimized height of n-order fitting using the fitting method. hite;n represents the optimized height of the nth iteration step by the proposed iteration method. N opt is the times of applying WTPO and is used to evaluate the total computation of the methods. It can be seen that the iteration method can obtain the convergence optimized results in 4 steps for various uref . hite;2 for various values of uref are close to hopt . The average error of objective value is only 1.13%, which is similar to the ones by 4-order fitting (1.19%) for the fitting method. Furthermore, the iteration method only needs to apply WTPO three times to obtain hite;2 , which is less than the enumeration method and the fitting method. As the number of iteration steps increases, hite;n is closer to hopt and the average error of objective value decreases. For the 4th iteration step, the average error of objective value is only 0.72%. Generally, the proposed iteration method can obtain better optimized height in less computational time than the fitting method. For the selected values of uref , the maximum PF is 44.9% and the minimum is 14.8%. Therefore, the height effect should be taken into consideration for more power output and less cost. The value of PF increases as the wind speed decreases. For a wind farm with lower wind speed, it has more potential improvement in the tower height matching problem. Fig. 4 shows the optimized layouts for the enumeration method and the iteration method for the situation of uref ¼ 10 m=s. It can be seen that the objective values are close, which is 1.019 for 63 m

Fig. 5. The wind rose of Case 2 (Chen et al., 2013a).

Fig. 6. TSMI curves of Case 2.

Table 4 Comparison of results for Case 2. uref (m/s)

10

PF

18.0%

Optimized height Relative error N opt

Enumeration

Fitting method (m)

Iteration method (m)

hopt (m)

hfit;4

hfit;5

hite;1

hite;2

hite;3

hite;4

74 – 101

78 2.60% 6

77 5.38% 7

72 3.33% 2

73 6.54% 3

74 0% 4

72 3.33% 5

44

K. Chen et al. / J. Wind Eng. Ind. Aerodyn. 132 (2014) 37–48

WTPOs for the iteration method to obtain the optimized height than other methods. Fig. 7 shows the optimized wind turbine layout of Case 2.

and is 1.017 for 66 m. The relative error calculated by Eq. (22) is only 0.2%. 5.1.2. Case 2 In Case 2, 16-directional wind with single speed of 10 m/s is studied. The wind rose of Case 2 is presented in Fig. 5. Fig. 6 shows the TSMI curves and Table 4 lists the optimized results by the three methods. For the fitting method, the best optimized height is 78 m and the relative error is 2.60%. The error of the TSMI curves around the optimal tower height between the enumeration method and the iteration method decreases as the iteration step increases. As hite;4 equals to hite;1 , the iteration method does not obtain a convergence optimized height in this situation. The best height among hite;n of each step is the final optimized height, which is the one of the third step 74 m. For Case 2, the iteration method finds the optimal height. Furthermore, it needs to apply fewer less

5.1.3. Case 3 In Case 3, multi-directional wind with variable speeds is considered. The velocities are in the range from 9 m/s to 11 m/s. The wind rose of Case 3 is shown in Fig. 8. Fig. 9 shows the TSMI curves of Case 3. Similar to the above cases, the error of the TSMI curves around the optimal height between the iteration method and the enumeration method decreases as the iteration step increases. Table 5 lists the optimized results by the three methods. For the fitting method, the best optimized height is 96 m and the error is 0.48%. While for the iteration method, the optimized height is 93 m when the iteration is convergence, which is very close to the optimal one 92 m. However, though the optimized height for the iteration method is closer to hopt than the one of the fitting method, the error of the best TSMI value for the iteration method is 1.67%, a little larger than the one for the fitting method. Due to the influences of the grids and the safe distances among the turbines, there are few extreme points around the maximum point in the TSMI curves, which can be seen in Fig. 9. Therefore, the result by the iteration method may converge to a certain local optimum. Through increasing the number of the grids in WTPO, the TSMI curves may be more smooth and the result by the iteration method may converge to the global optimum. Again, the value of Nopt for the iteration method is less than other methods. Fig. 10 shows the optimized layouts for the enumeration method and the iteration method.

2000

1500

1000

500

0

5.2. Case 4: a case with realistic models

0

500

1000

1500

2000

Fig. 7. Optimized wind turbine layout of Case 2.

In this section, a case with realistic models is considered to further test the proposed iteration method. Wind condition with 16 wind directions is considered and each direction satisfy a

Fig. 8. The wind rose of Case 3 (Chen et al., 2013a).

Fig. 9. TSMI curves of Case 3.

Table 5 Comparison of results for Case 3. uref (m/s)

9–11

PF

9.8%

Optimized height Relative error N opt

Enumeration

Fitting method (m)

Iteration method (m)

hopt (m)

hfit;4

hfit;5

hite;1

hite;2

hite;3

hite;4

92 – 101

79 5.13% 6

96 0.48% 8

87 4.47% 2

91 2.24% 3

93 1.67% 4

93 1.67% 5

K. Chen et al. / J. Wind Eng. Ind. Aerodyn. 132 (2014) 37–48

2000

2000

1500

1500

1000

1000

500

500

0

0

500

1000

1500

2000

0

0

500

1000

45

1500

2000

Fig. 10. Optimized wind turbine layouts of Case 3: (a) h ¼ 92 m and (b) h ¼ 93 m.

Table 6 Parameters of the wind turbine in Case 4.

k

3 2 1

0

100

200

300

400

0

100

200

300

400

c(m/s)

15 10 5

p

0.2

Parameter

Value

Diameter (m) Rated power output (kW) Rated speed (m/s) Cut-in speed (m/s) Cut-out speed (m/s) Thrust coefficient Tower height (m) Turbine cost ð€Þ Tower cost ð€=mÞ

40 680 13.0158 2.5 28 0.8888 46–78 593,867 1500

0.1 0

Table 7 Comparison of results for Case 4.

0

100

200

300

400

Direction(degrees) Fig. 11. Wind data for Case 4.

Weibull distribution. The parameters of the Weibull distributions are generated randomly, as shown in Fig. 11. The wind turbine in the study of Chen et al. (2013b) is used, which has a tower height ranging from 46 m to 78 m. The parameters of the turbines are listed in Table 6. The power curve of the turbine are shown as 8 u o2:5 or u Z 28 > < 0; PðuÞ ¼ 0:3084  u3 ; 2:5 r u o 13:0185 ð23Þ > : 680; 13:0185 r u o 28 The cost model in the study of Chen et al. (2013b) can be expressed as below Cost ¼ C 1 þC 2 h

ð24Þ

where C1 is turbine cost and C2 is tower cost per meter. The objective is chosen as Objective ¼

P tot ∑N i Costi

Parameter

Enumeration

Iteration method

Optimized height (m) Relative error N opt

55 – 33

55 0% 2

wind directions. The initial standard height is chosen as 46 m. Table 7 shows the optimized results. Through the enumeration method, the maximum objective value corresponds to the tower height of 55 m and the minimum one corresponds to 78 m. The PF value is 1.2%. It only needs to apply twice of optimizations to obtain the convergence optimized height for iteration method. The optimized height is 55 m, equals to the optimal one. Fig. 12 shows the optimized layouts for the tower heights of 55 m and 78 m. When the tower height is large, the increase of the cost overcomes the increase of the power output. Thus, the objective value decreases. It can be seen that the iteration method has good performance for Case 4. 5.3. Case 5: a case over complex terrain

ð25Þ

where N is the number of the wind turbine placed in the wind farm. The domain of the wind farm is chosen as a square area with the size of 1000 m  1000 m, meshed by 20  20. The ground roughness of the wind farm is 0.3 m. The number of wind turbines to be placed in the domain is 10. The distance among the wind turbine is fixed as 200 m, which is five times of the wind rotor diameter. The tower height matching is solved by the enumeration method and the proposed iteration method. For the enumeration method, the height with the best objective value is treated as the optimal height. The relative error is calculated by Eq. (22). For Case 4, wind speeds satisfy various Weibull distributions for various

In this section, a case over complex terrain is used to validate the proposed iteration method. A typical complex terrain with the size of 2000 m  2000 m is introduced and the contour of altitude is shown in Fig. 13. Consider the same wind cases and the same wind turbine in Case 3. The number of wind turbines to be placed in the site is 30. The flow fields of the complex terrain of various wind directions are obtained through CFD calculation. The particle wake flow model is introduced to calculate the turbine wake flow. For the wind turbine in Case 3, the parameters β, σ and γ of the particle wake model are valued as 0.65, 0.3 and 0.005, respectively (Song et al., 2012). The greedy algorithm is used to optimize the wind turbine positioning problem for a specified tower height. Since the particle wake model contains randomness and needs

46

K. Chen et al. / J. Wind Eng. Ind. Aerodyn. 132 (2014) 37–48

1000

1000

800

800

600

600

400

400

200

200

0

0

200

400

600

800

1000

0

0

200

400

600

800

1000

Fig. 12. Optimized wind turbine layouts of Case 4: (a) h ¼ 55 m and (b) h ¼ 78 m.

Fig. 13. The contour of altitude of the complex terrain. Fig. 14. The comparison of TSMI curves of Cases 3 and 5.

more computational time to obtain the wake flow result, only the locating stage of greedy algorithm is considered. The number of optimization grids is chosen as 50  50. Other parameters are the same as Case 3. Fig. 14 shows the comparison of the TSMI curves of the flat terrain and the complex terrain for the wind cases of Case 3. For complex terrain, the wind speed for a specified height in the wind farm is not uniform. Large wind speeds are obtained at some parts of low height due to the disturbance of the terrain. When the height is low, the TSMI value of the complex terrain is larger than the one of the flat terrain. When the height is large, the power outputs and the TSMI values of both terrains are similar. Therefore, the optimal tower height of the complex terrain is lower and the relevant TSMI value is larger than the ones of the flat terrain. PF value is 6% for Case 5. As the minimum TSMI value of complex terrain increases, PF value of complex terrain decreases compared to flat terrain. Table 8 shows the optimized results of the enumeration method and the iteration method for Case 5. It can be seen that the optimized tower height for the 4th iteration steps is 82 m, close to the optimal one 76 m by the enumeration method. The relative error of the TSMI value is only 0.08%. Fig. 15 shows the optimized layouts of 76 m and 82 m. In the wind farm, the location with higher altitude usually has higher wind speed, so some wind turbines are placed on the top of the hill. Besides, most of the wind turbines are located at the edge of the wind farm to enlarge the distances from others, reducing the influence of the wake effect.

Table 8 Optimized results for Case 5. uref (m/s)

9–11

PF

6.0%

Optimized height Relative error N opt

Enumeration

Iteration method (m)

hopt (m)

hite;1

hite;2

hite;3

hite;4

76 – 101

75 0.33% 2

77 0.33% 3

82 0.08% 4

83 0.34% 5

It can be concluded from Case 5 that the iteration method can also be applied to complex terrain.

6. Conclusions In this paper, the tower height matching for wind turbine positioning optimization is studied. The wind speed in wind farm is modeled by the power law with 1/7 order to the height. The wake flow is calculated by the linear wake model over flat terrain and by the particle wake model over complex terrain. The power curve of the wind turbine is characterized by a simplified 3-order power law with power control mechanism. The wind turbine

K. Chen et al. / J. Wind Eng. Ind. Aerodyn. 132 (2014) 37–48

47

Fig. 15. Optimized wind turbine layouts of Case 5: (a) h ¼ 76 m and (b) h ¼ 82 m.

positioning optimization at a specified height is solved by greedy algorithm, with the target of maximizing the Turbine-Site Matching Index (TSMI). Assuming that the optimized layout for each tower height is the same, the iteration method is developed to obtain the approximated optimal height. The convergence of the proposed iteration method is discussed through the mathematical analysis. The numerical cases over flat terrain and complex terrain are introduced to test the proposed method. The results show that the TSMI curve generated by a selected standard height is close to the accurate TSMI curve, especially around the selected height. For most wind situations, the proposed iteration method can obtain the optimized tower height in less computational time and the errors of the results are less than the one in previous studies. As a conclusion, the iteration method is an effective way to solve tower height matching problem in wind farm design.

Assumption 1. All the functions in Eq. (A.1) are at least secondorder differentiable, that is TðhÞ;

T a ðhÞ;

This research is supported by the International Scientific and Technological Cooperation Program of China (No. 2011DFG13020), the China Postdoctoral Science Foundation (2013M530043) and the National High-Tech R&D Program (863 Program) of China (No. 2007AA05Z426).

Appendix A. Convergence of the iteration method

PðhÞ ; TðhÞ ¼ CðhÞ

P a ðhÞ T a ðhÞ ¼ ; CðhÞ

P a ðhÞ ¼ PðhÞ  pðhÞ

ðA:1Þ

CðhÞ A C 2

pðhÞ;

pðhÞ Z 0; p″ðhÞ Z 0; h A S p0 ðh0 Þ ¼ 0; p″ðh0 Þ Z 0

ðA:5Þ

when h0 rhm when h0 4hm

Assumption 4. Let FðhÞ ¼ pðhÞC 0 ðhÞ  p0 ðhÞCðhÞ, which monotonically decreases in S, that is F 0 ðhÞ ¼ pðhÞC″ðhÞ p″ðhÞCðhÞ r0;

hAS

P 0 ðh0 ÞCðh0 Þ  Pðh0 ÞC 0 ðh0 Þ ½Cðh0 Þ2

40

ðA:7Þ

Consider Eq. (A.2) and Assumption 3, it obtains T 0a ðh0 Þ ¼

½P 0 ðh0 Þ  p0 ðh0 ÞCðh0 Þ  ½Pðh0 Þ  pðh0 ÞC 0 ðh0 Þ ½Cðh0 Þ2 P ðh0 ÞCðh0 Þ  Pðh0 ÞC 0 ðh0 Þ

¼

¼ T 0 ðh0 Þ 4 0

where P(h) and Pa(h) are the accurate power output and the approximated power output through translating the optimized layout of the standard height, respectively. C(h) is the cost function, expressed by Eq. (19). p(h) is the error between P a ðhÞ and P(h). Denote the maxima of T(h) and T a ðhÞ as hm and ham, respectively. The standard height of the optimized turbine layout for P a ðhÞ is h0. We have

T 0a ðhm Þ ¼

Tðh0 Þ ¼ T a ðh0 Þ;

According to Assumption 4, F 0 ðhÞ r0. As h0 o hm , then

T 0 ðhm Þ ¼

pðh0 Þ ¼ 0

P 0 ðhm ÞCðhm Þ  Pðhm ÞC 0 ðhm Þ 2

½Cðhm Þ

¼ 0;

T 0a ðham Þ ¼ 0

ðA:2Þ

In order to prove the convergence of the iteration method, the following assumptions are introduced.

ðA:6Þ

Based on these four assumptions, the convergence of the iteration method can be proved as follows: (1) When h0 o hm , then

0

Pðh0 Þ ¼ P a ðh0 Þ;

ðA:3Þ

Assumption 3. p(h) is a strictly convex function and h0 is the minima of the p(h), that is

T 0 ðh0 Þ ¼ Denote the accurate TSMI curve and the approximated TSMI curve as T(h) and T a ðhÞ, respectively, which are expressed as

P a ðhÞ;

Assumption 2. T(h) is a strictly concave function. Therefore, around hm, it obtains ( 0 T ðh0 Þ 40; when h0 o hm ðA:4Þ T 0 ðh0 Þ o0; when h0 4 hm

where ( ½h0 ; hm ; S: ½hm ; h0 ;

Acknowledgments

PðhÞ;

½Cðh0 Þ2

½P 0 ðhm Þ  p0 ðhm ÞCðhm Þ  ½Pðhm Þ  pðhm ÞC 0 ðhm Þ ½Cðhm Þ2 pðhm ÞC ðhm Þ  p ðhm ÞCðhm Þ 0

¼

T 0a ðhm Þ ¼

ðA:8Þ

0

½Cðhm Þ2 Fðhm Þ ½Cðhm Þ2

r

Fðh0 Þ ½Cðhm Þ2

¼

Fðhm Þ ½Cðhm Þ2

¼0

ðA:9Þ

ðA:10Þ

According to the analysis above, we have T 0a ðh0 Þ 40;

T 0a ðhm Þ r 0

ðA:11Þ

48

K. Chen et al. / J. Wind Eng. Ind. Aerodyn. 132 (2014) 37–48

Therefore, it obtains h0 o ham r hm ;

when h0 o hm

ðA:12Þ

(2) When h0 4 hm , then Tðh0 Þ ¼

P 0 ðh0 ÞCðh0 Þ  Pðh0 ÞC 0 ðh0 Þ ½Cðh0 Þ2

o0

ðA:13Þ

Consider Eq. (A.2) and Assumption 3, it obtains T 0a ðh0 Þ ¼

½P 0 ðh0 Þ p0 ðh0 ÞCðh0 Þ ½Pðh0 Þ pðh0 ÞC 0 ðh0 Þ ½Cðh0 Þ2 P ðh0 ÞCðh0 Þ  Pðh0 ÞC 0 ðh0 Þ 0

¼

½Cðh0 Þ2

T 0a ðhm Þ ¼

ðA:14Þ

½P 0 ðhm Þ  p0 ðhm ÞCðhm Þ  ½Pðhm Þ pðhm ÞC 0 ðhm Þ

½Cðhm Þ2 pðhm ÞC ðhm Þ  p ðhm ÞCðhm Þ 0

¼

¼ T 0 ðh0 Þ o 0

0

½Cðhm Þ2

ðA:15Þ

According to Assumption 4, F 0 ðhÞ r 0. As h0 4 hm , then T a 0 ðhm Þ ¼

Fðhm Þ ½Cðhm Þ2

Z

Fðh0 Þ ½Cðhm Þ2

¼0

ðA:16Þ

According to the analysis above, we have T 0a ðhm Þ Z0;

T 0a ðh0 Þ o 0

ðA:17Þ

Therefore, it obtains hm r ham o h0 ;

when h0 4 hm

ðA:18Þ

(3) When h0 ¼ hm , then T 0a ðh0 Þ ¼ ¼

½P 0 ðh0 Þ p0 ðh0 ÞCðh0 Þ ½Pðh0 Þ pðh0 ÞC 0 ðh0 Þ ½Cðh0 Þ2 P 0 ðh0 ÞCðh0 Þ  Pðh0 ÞC 0 ðh0 Þ ½Cðh0 Þ2

¼ T 0 ðh0 Þ ¼ T 0 ðhm Þ ¼ 0

ðA:19Þ

Therefore, ham ¼ h0 ¼ hm . According to the analysis above, based on the four assumptions, we have 8 > < h0 o ham r hm ; when h0 o hm ham ¼ hm ; when h0 ¼ hm ðA:20Þ > : h r h o h ; when h 4 h m am m 0 0 Therefore, from the process of the iteration, a sequence can be constructed as follows: ham;k þ 1 ¼ φðham;k Þ;

ham;1 ¼ h0

which satisfies 8 h o ham;k þ 1 r hm ; > < am;k ham;k ¼ ham;k þ 1 ¼ hm ; > : h rh oh ; m

am;k þ 1

am;k

ðA:21Þ

when h0 ohm when h0 ¼ hm

ðA:22Þ

when h0 4hm

where φ represents the iteration process. It can be seen that fham;k g is bounded and monotonous. According to calculus theory, fham;k g is convergence to a real number which is within the range ½h0 ; hm  or ½hm ; h0 . However, it can only guarantee the convergence, but cannot guarantee that it can converge to hm. Nevertheless, when

h0 is chosen to be close to hm, the convergence result can be close to hm. The convergence is based on the four assumptions. Theoretically, the tower height and the positions of turbines are continuous, so Assumption 1 is satisfied. A large amount of test cases show that T(h) is a strictly concave function for most situations. That is, Assumption 2 is satisfied for most situations. In the present study, C″ðhÞ ¼ 0 and CðhÞ 4 0. According to Assumption 3, p″ðhÞ Z 0. Therefore, F 0 ðhÞ ¼  p″ðhÞCðhÞ r 0. Assumption 4 is satisfied when Assumption 3 is satisfied. In reality, p(h) is not a strictly convex function and Assumption 3 is not satisfied for all the situations. However, after performing a large amount of wind cases, it is found that the final height of the iteration method can converge to the optimal one or at least around the optimal one for most situations. References Albadi, M.H., El-Saadany, E.F., 2010. Optimum turbine-site matching. Energy 35, 3593–3602. Burton, T., Sharpe, D., jenkins, N., Bossanyi, E., 2001. Wind Energy Handbook. John Wiley and Sons, Ltd., New York, USA. Chen, D., Zhang, W., 2007. Exploitation and research on wind energy (in Chinese). Energy Conserv. Technol. 4, 339–343. Chen, K., Song, M., Zhang, X., 2013a. The investigation of tower height matching optimization for wind turbine positioning in the wind farm. J. Wind Eng. Ind. Aerodyn. 114, 83–95. Chen, Y., Li, H., Jin, K., Song, Q., 2013b. Wind farm layout optimization using genetic algorithm with different hub height wind turbines. Energy Convers. Manag. 70, 56–65. Gonzalez, J.S., Rodriguez, A., Mora, J.C., Santos, J.R., Payan, M.B., 2010. Optimization of wind farm turbines layout using an evolutive algorithm. Renew. Energy 35, 1671–1681. Grady, S.A., Hussaini, M.Y., Abdullah, M.M., 2005. Placement of wind turbines using genetic algorithms. Renew. Energy 30, 259–270. Kusiak, A., Song, Z., 2010. Design of wind farm layout for maximum wind energy capture. Renew. Energy 35, 685–694. Marmidis, G., Lazarou, S., Pyrgioti, E., 2008. Optimal placement of wind turbines in a wind park using monte carlo simulation. Renew. Energy 33, 1455–1460. Mora, J.C., Baron, J.M.C., Santos, J.M.R., Payan, M.B., 2007. An evolutive algorithm for wind farm optimal design. Neurocomputing 70, 2651–2658. Mosetti, G., Poloni, C., Diviacco, B., 1994. Optimization of wind turbine positioning in large windfarms by means of a genetic algorithm. J. Wind Eng. Ind. Aerodyn. 51, 105–116. Ozturk, U.A., Norman, B.A., 2004. Heuristic methods for wind energy conversion system positioning. Electr. Power Syst. Res. 70, 179–185. Rivas, R.A., Clausen, J., Hansen, K.S., Jensen, L.E., 2009. Solving the turbine positioning problem for large offshore wind farms by simulated annealing. Wind Eng. 33, 287–298. Song, M., Chen, K., He, Z., Zhang, X., 2012. Wake flow model of wind turbine using particle simulation. Renew. Energy 41, 185–190. Song, M., Chen, K., He, Z., Zhang, X., 2013. Bionic optimization for micro-siting of wind farm on complex terrain. Renew. Energy 50, 551–557. Song, Y., 2012. The investigation of the development status and countermeasure of wind energy resources in china (in Chinese). J. Electr. Power 27, 494–497. Tong, X., Dong, Y., 2012. Wind energy resources and development of wind power industry in china (in Chinese). Energy Res. Util. 6, 25–27. Wan, C., Wang, J., Yang, G., Li, X., Zhang, X., 2009. Optimal micro-siting of wind turbines by genetic algorithms based on improved wind and turbine models. In: Joint 48th IEEE Conference on Decision and Control and 28th Chinese Control Conference, Shanghai, P.R. China. pp. 5092–5096. Wan, C., Wang, J., Yang, G., Zhang, X., 2010. Optimal micro-siting of wind farms by particle swarm optimization. Lecture Notes in Computer Science vol. 6145, 198–205. Zhang, C., Hou, G., Wang, J., 2011. A fast algorithm based on the submodular property for optimization of wind turbine positioning. Renew. Energy 36, 2951–2959.