Integrated Wind Farm Power Curve and Power Curve Distribution Function Considering the Wake E ﬀ ect and Terrain Gradient

: This work presents a computational method for the simulation of wind speeds and for the calculation of the statistical distributions of wind farm (WF) power curves, where the wake e ﬀ ects and terrain features are taken into consideration. A three-parameter (3-P) logistic function is used to represent the wind turbine (WT) power curve. Wake e ﬀ ects are simulated by means of the Jensen’s wake model. Wind shear e ﬀ ect is used to simulate the inﬂuence of the terrain on the WTs located at di ﬀ erent altitudes. An analytical method is employed for deriving the probability density function (PDF) of the WF power output, based on the Weibull distribution for describing the cumulative wind speed behavior. The WF power curves for four types of terrain slopes are analyzed. Finally, simulations applying the Monte Carlo method on di ﬀ erent sample sizes are provided to validate the proposed model. The simulation results indicate that this approximated formulation is a possible substitute for WF output power estimation, especially for the scenario where WTs are built on a terrain with gradient.


Introduction
To date, wind energy has become a promising alternative and renewable energy resource to handle the problem of energy crisis and environmental deterioration.Along with the rapid development and extensive exploitation of wind energy throughout the world, calculations of the output power of wind farm (WF)s is widely required not only for the WF layout optimization problem [1,2] but for the WF output monitoring [3][4][5], analysis [6], prediction [7,8], and control [9][10][11][12][13] problems.However, in most of the literature, the wind turbine (WT) power curves are mathematically expressed by means of a piecewise function.For a certain inflow wind speed, wind speeds experienced by different WTs on a WF differ for several reasons.One of the major influencing factors is the wake effect and another may be the altitude differences of the WTs' positions.Therefore, representation of the WF power curves should not be simply acquired as the summation of power curves produced by the individual WTs.
First of all, the power curves of a WT not only indicate the performance of the WT, but also serve as a reference for prediction, if accurately modeled [14,15].Therefore, in order to obtain WF power curves, accurate models for WT power curves must be dealt with first.The WT power curve models are commonly classified as parametric or nonparametric.In [4] it is proven that the parametric least squares model and the nonparametric k-nearest neighbor model present high accuracy.In [5] a parametric power curve was developed, that could be applied for control, monitoring, prediction, and optimization of WF performance.A comparison of several WT power curve modeling methods and a critical analysis of then is available in [16].In [17] the estimation of the WT power curve was operated on the basis of cluster center fuzzy logic modelling.In 2016, Pelletier [18] used an artificial neural network to model WT power curves.In 2017, Ouyang [19] proposed a data partitioning and mining approach for WT power curve modeling.In 2018, Zhao et al. [20] developed a data-driven correction method for the refinement of a WF power curve under wind curtailment and the work in [20] presented the ability to be directly used in different cases without the need to tune any parameters, for WTs and also for WFs.
On the basis of the extensive research on WT power curves, many previous works have established WF power curve distributions taking wake effect into consideration.In 2014, Shi et al. [21] developed a solution based on a probability analysis for modeling WF power output.In the cited paper, wake effect was taken into consideration by using the Jensen's model and the WF power probability density function (PDF) was deduced with the help of an analytical method.Feijóo et al. were inspired by [21] and combined this with [22] to deduce WF power curve models and power probability density models with three and four parameters in 2017, [23] and [24] respectively.One of the advantages of using these models is that they are not piecewise, which differs from other proposals, such as [21].Moreover, the wake deficit of the whole WF is synthesized to be represented by one wake coefficient which avoids the computational complexity of calculating the wake deficit one by one.In addition, in [25], four major wind directions were taken into account for the approximation of the WF output power calculation which further refined this model.
As is pointed out in [23][24][25], it is advantageous to utilize these models in view of their simplicity as compared with other proposals and to preserve accuracy.In addition, wake effects are already synthesized in the expression of total power and WF power PDF.Nonetheless, the major deficiency of these models is that they can only be applied in the case where the WTs are installed on flat land.If the WF is located on a hilly area where the positions of the WTs form a gradient, then these models cannot be directly applied to calculate the total output power of the WF which results in an increase of the computation error.Thus, in this case, if these simplified models are to be used a modification needs to be added beforehand.
This work is based on the advantages of the aforementioned models and is organized as follows: First, Jensen's wake model is introduced and then the analytical wake model is given which takes the height differences between WTs into consideration.Secondly, an integrated three parameter (3-P) power distribution function of a WF is established which is classified into four sub-models according to the value of the hill's slope.Finally, this model is proven by simulations and comparison with the Monte Carlo simulation results.

Wake Model
On a WF with multiple rows of WTs, the presence of a WT affects the wind speed downstream.This aerodynamic fact is called wake effect.The wake effect of a WT is depicted by using the Jensen's model [26,27] as shown in Figure 1.
where, R r is the length of a WT blade and R w is the radius of the wake area caused by WT1 at distance x, ρ is the air density, v 0 is the ambient wind speed, and K is the wake decay constant as is shown in Figure 1.
Energies 2019, 12, 2482 3 of 14 where, Rr is the length of a WT blade and Rw is the radius of the wake area caused by WT1 at distance x, ρ is the air density, v0 is the ambient wind speed, and K is the wake decay constant as is shown in Figure 1.The wind speed, vx, on WT2 is calculated according to Equation 2.
where,  is the wake coefficient and Ct are the thrust coefficients of the WT.

Shear Effect
In real situations, wind speed depends on height, which is known as wind shear, according to the influence of the boundary layer fluid flow.As depicted and shown in Figure 2, the wind vertical distribution curve exhibits a certain pattern, such as a logarithmic curve or a cubic curve.Many factors influence the shape of this curve, majorly the roughness of ground.The rougher the ground is, the slower the wind close to the ground is.The wind speed, v x, on WT2 is calculated according to Equation (2).
where, ξ is the wake coefficient and C t are the thrust coefficients of the WT.

Shear Effect
In real situations, wind speed depends on height, which is known as wind shear, according to the influence of the boundary layer fluid flow.As depicted and shown in Figure 2, the wind vertical distribution curve exhibits a certain pattern, such as a logarithmic curve or a cubic curve.Many factors influence the shape of this curve, majorly the roughness of ground.The rougher the ground is, the slower the wind close to the ground is.
where, Rr is the length of a WT blade and Rw is the radius of the wake area caused by WT1 at distance x, ρ is the air density, v0 is the ambient wind speed, and K is the wake decay constant as is shown in Figure 1.The wind speed, vx, on WT2 is calculated according to Equation 2.
where,  is the wake coefficient and Ct are the thrust coefficients of the WT.

Shear Effect
In real situations, wind speed depends on height, which is known as wind shear, according to the influence of the boundary layer fluid flow.As depicted and shown in Figure 2, the wind vertical distribution curve exhibits a certain pattern, such as a logarithmic curve or a cubic curve.Many factors influence the shape of this curve, majorly the roughness of ground.The rougher the ground is, the slower the wind close to the ground is.The distribution of the obtained wind shear values is effectively approximated by the following expression [28]: where, v(z) is the wind speed at height z, v(z ref ) is the wind speed at the height z ref where it is measured, and α h is the approximation coefficient and is obtained according to Equation ( 5) for 0.001 m < z 0 < 10 m.

Wake Analysis and Terrain Gradient
First, let us suppose that the WF with a by b WTs is located at a terrain with a gradient of s.The distance between two adjacent rows is d row .There are four scenarios in total with different gradient values: Scenario 1 WF on a flat area.
The WF is considered to be located on a flat area (Figure 3), if the slope of the terrain, s, is equal to zero.This scenario was already analyzed in [21] and the conclusion of [21] is directly applied.


where, v(z) is the wind speed at height z, v(zref) is the wind speed at the height zref where it is measured, and h  is the approximation coefficient and is obtained according to Equation 5    2 00 0.24 0.096 log 0.016 log for 0.001 m < z0 < 10 m.

Wake Analysis and Terrain Gradient
First, let us suppose that the WF with a by b WTs is located at a terrain with a gradient of s.The distance between two adjacent rows is drow.There are four scenarios in total with different gradient values: Scenario 1 WF on a flat area.The WF is considered to be located on a flat area (Figure 3), if the slope of the terrain, s, is equal to zero.This scenario was already analyzed in [21] and the conclusion of [21] is directly applied.
In this case, the WTs share the same height.Thus, vx is directly calculated using Equation 2and the wake coefficient of the ith WT is: The wind speed at the ith WT is: The mean value 0 w  for all the a rows of the WTs' wake coefficients Scenario 2 WF over a hill with gentle gradient.In this case, the WTs share the same height.Thus, v x is directly calculated using Equation ( 2) and the wake coefficient of the ith WT is: The wind speed at the ith WT is: The mean value ξ w0 for all the a rows of the WTs' wake coefficients 1, ξ, ξ 2 , • • • , ξ a−1 is calculated as: Scenario 2 WF over a hill with gentle gradient.
When the slope of the terrain is gentle, the WTs located at downstream positions are totally affected by the wake of upstream WTs (Figure 4).When the slope of the terrain is gentle, the WTs located at downstream positions are totally affected by the wake of upstream WTs (Figure 4).
In other words, the slope must meet the following condition: - Namely, That is to say, if the slope of the hill, s, is lower or equal to the wake decay constant, K, then the wind speed of the downstream WT is calculated as in Equation 11.In other words, the slope must meet the following condition: Namely, Energies 2019, 12, 2482 5 of 14 That is to say, if the slope of the hill, s, is lower or equal to the wake decay constant, K, then the wind speed of the downstream WT is calculated as in Equation (11).
where, x = d row and z = d row s The wake coefficient of the ith WT is: Thus, the wind speed of the ith WT is: For simplicity, let where, Scenario 3 WF over a hill with moderate gradient.When the slope of the terrain is moderate, the WTs located at downstream positions are partially affected by the wake of upstream WTs (Figure 5).When the slope of the terrain is moderate, the WTs located at downstream positions are partially affected by the wake of upstream WTs (Figure 5).
This means the slope must meet the following condition: - In this scenario, the downstream WT2 vertically deviates the upstream WT1 due to the shape of the terrain.WT2 may be incompletely covered by WT1′s wake.The area of the shadow on WT2 covered by WT1′s wake is shown in Figure 5 when considering terrain gradient.This means the slope must meet the following condition: Namely, In this scenario, the downstream WT2 vertically deviates the upstream WT1 due to the shape of the terrain.WT2 may be incompletely covered by WT1 s wake.The area of the shadow on WT2 covered by WT1 s wake is shown in Figure 5 when considering terrain gradient.
The area of the shadow in Figure 6 is calculated as in [29]: Energies 2019, 12, 2482 6 of 14 In this scenario, the downstream WT2 vertically deviates the upstream WT1 due to the shape of the terrain.WT2 may be incompletely covered by WT1′s wake.The area of the shadow on WT2 covered by WT1′s wake is shown in Figure 5 when considering terrain gradient.The area of the shadow in Figure 6 is calculated as in [29]: arccos arccos sin arccos 2 2 2 The height difference of two WTs is calculated as follows: The wind speed of the downstream WT is calculated as follows: where, Thus, the wind speed of the ith WT is: The height difference of two WTs is calculated as follows: The wind speed of the downstream WT is calculated as follows: Thus, the wind speed of the ith WT is: The wake coefficient of the ith WT is: For simplicity, let where, Scenario 4 WF over a hill with steep gradient.When the slope of the terrain is steep, wakes do not affect the WTs located at downstream positions (Figure 7).
This means the slope must meet the following condition: Namely, In this case, there is no wake influence between any two WTs.Thus, the only factor influencing the wind speed is shear effect.
The wind speed of the ith WT is calculated as follows: Energies 2019, 12, 2482 7 of 14 It should be noted that this scenario excludes the case where the gradient is extremely steep and the mountain itself blocks airflow while wake effect between WTs at different heights does not occur.
For simplicity, let where, Scenario 4 WF over a hill with steep gradient.When the slope of the terrain is steep, wakes do not affect the WTs located at downstream positions (Figure 7).
This means the slope must meet the following condition: Namely, In this case, there is no wake influence between any two WTs.Thus, the only factor influencing the wind speed is shear effect.
The wind speed of the ith WT is calculated as follows: According to the Betz limit [30], if the WT extracts 100% of the wind energy, the air behind the WT does not move, which creates a "wall" and stops further wind from passing through the WT.At a fixed free wind speed, the distance of the WT's blade to the "wall" determines the inflow wind speed on the WT blades which is reduced by a restrained path of outflow.

3-P WT PDF Model
Three models for WT power curves and WT PDFs were proposed in [22], on the basis of the logistic function, which is continuous.The one called 3-P-DP is mathematically expressed by means of Equation (26).
where, P r is the rated power of the WT, , v ip the value of the wind speed at the inflection point of the power curve, and v is the inflow wind speed of the WT.
If the inflow wind speed is regarded as a random variable which follows a Weibull distribution with a scale parameter C and a shape parameter k, i.e., , then the PDF of the WF output power is obtained using the change of variable technique [31]: The advantage of utilizing Equations ( 26) and ( 27) is that they are continuous functions, unlike the piecewise function which has been proposed in other papers.

3-P WF PDF Model
Consider a WF with a row, b column layout for a WF, and let d row be the space existing between two different WTs.In this case, an equivalent model for the WF power is proposed by assuming all the WTs of one column are located in a row.The proposal supposes a unidirectional wind speed.
Scenario 1 WF on a flat area (s = 0).When the WF is built on a flat area, the total output power of the whole WF is calculated as follows: where, ξ w0 is calculated using Equation ( 8).Scenario 2 WF over a hill with gentle gradient (0 < s ≤ K).
The total output power of the WF is calculated as: where, ξ w1 is the equivalent wake effect coefficient of the WF with gentle gradient, which is calculated as follows: where, C 1 is the same parameter in Equation ( 14).Scenario 3 WF over a hill with medium gradient (K < s ≤ 2R r d row + K).Similar to Scenario 2, the total output power of the WF is calculated as: where, ξ w2 is the equivalent wake effect coefficient of the WF over a hill with medium gradient, which is calculated as follows: where, C 2 is the same parameter in Equation (22).Scenario 4 WF over a hill with steep gradient (s > 2R r d row + K).The power function is different from the above scenarios since there is no wake effect in this scenario.The total output power of the WF is calculated as: Energies 2019, 12, 2482 9 of 14

WF Power PDF
The PDF of the WF output power, now considering wake effect and terrain gradient, for Scenario 1 to Scenario 3 is expressed in Equation (34), by applying the change of variable technique [31].The PDF of the total WF output power for Scenario 4 is difficult to deduce.Thus, it is not displayed here.

Implementation of the WF Output Power Calculation Function
The implementation of this integrated power distribution function is illustrated by the flow chart in Figure 8.There are five main steps in total: Step 1 Initialization of the program and input of all the parameters.
Step 2 Calculation of the terrain's slope to determine the scenario.
Step 3 Calculation of the whole WF's output power using the model for the scenario.
Step 4 Calculation of the WF output power's PDF.
Step 5 Output of the calculated power and PDF of the whole WF.

WF Power PDF
The PDF of the WF output power, now considering wake effect and terrain gradient, for Scenario 1 to Scenario 3 is expressed in Equation 34, by applying the change of variable technique [31].The PDF of the total WF output power for Scenario 4 is difficult to deduce.Thus, it is not displayed here.

 
 

Implementation of the WF Output Power Calculation Function
The implementation of this integrated power distribution function is illustrated by the flow chart in Figure 8.There are five main steps in total: Step 1 Initialization of the program and input of all the parameters.
Step 2 Calculation of the terrain's slope to determine the scenario.
Step 3 Calculation of the whole WF's output power using the model for the scenario.
Step 4 Calculation of the WF output power's PDF.
Step 5 Output of the calculated power and PDF of the whole WF.

Case Study
Simulations of a 50 MW WF with 25 two MW WTs, in a five by five layout, are presented in this section.The constants are given in Table 1: The slope of the hill, s, changes from s = 0 for Scenario 1; s = 0.02, s = 0.05, and s = 0.07 for Scenario 2; s = 0.10, s = 0.12, and s = 0.15 for Scenario 3; and s = 1.8 for Scenario 4.

Simulation of the WF Power Curve Functions
The simulation results are given below: It is clearly seen from Figure 9 that, generally, when the gradient of the WF's location changes, at the same inflow wind speed, the WF with the same layout and same WT type will produce different output power, that is to say, when s grows, P WF also increases in accordance.This phenomenon is explained from two aspects.First, with an increase of the gradient, the coverage areas of wake gradually decrease which gives rise to the reduction of wake loss.Secondly, with the rise of s, an increase of the WTs' elevation brings about higher wind speed.These two factors synthetically lead to an increase of the WF power output.Another phenomenon is that when s changes from 0 to 0.07, the curves of P WF0 and P WF1 vary dramatically.Nevertheless, when s changes from 0.10 to 1.8, the curves of P WF2 and P WF3 vary to a much lower degree.This implies that the power curve of the WF has a higher sensitivity to the terrain's gradient if s is small.When the wake is partially covered or not covered, the power curve of the WF is altered much less through a change in s.

Case Study
Simulations of a 50 MW WF with 25 two MW WTs, in a five by five layout, are presented in this section.The constants are given in Table 1: The slope of the hill, s, changes from s = 0 for Scenario 1; s = 0.02, s = 0.05, and s = 0.07 for Scenario 2; s = 0.10, s = 0.12, and s = 0.15 for Scenario 3; and s = 1.8 for Scenario 4.

Simulation of the WF Power Curve Functions
The simulation results are given below: It is clearly seen from Figure 9 that, generally, when the gradient of the WF's location changes, at the same inflow wind speed, the WF with the same layout and same WT type will produce different output power, that is to say, when s grows, PWF also increases in accordance.This phenomenon is explained from two aspects.First, with an increase of the gradient, the coverage areas of wake gradually decrease which gives rise to the reduction of wake loss.Secondly, with the rise of s, an increase of the WTs' elevation brings about higher wind speed.These two factors synthetically lead to an increase of the WF power output.Another phenomenon is that when s changes from 0 to 0.07, the curves of PWF0 and PWF1 vary dramatically.Nevertheless, when s changes from 0.10 to 1.8, the curves of PWF2 and PWF3 vary to a much lower degree.This implies that the power curve of the WF has a higher sensitivity to the terrain's gradient if s is small.When the wake is partially covered or not covered, the power curve of the WF is altered much less through a change in s.

Validation by Means of Monte Carlo Simulations
The Monte Carlo simulation method [32] is utilized for comparisons of precision and computation time.The parameters of the Weilbull distribution function are set to be C = 8m/s and k = 2.The cumulative distribution function (CDF) is chosen instead of the PDF in Equation 34, because it facilitates the comparison.

Validation by Means of Monte Carlo Simulations
The Monte Carlo simulation method [32] is utilized for comparisons of precision and computation time.The parameters of the Weilbull distribution function are set to be C = 8m/s and k = 2.The cumulative distribution function (CDF) is chosen instead of the PDF in Equation (34), because it facilitates the comparison.
The basic steps of the Monte Carlo simulation [32] are as follows: Step 1 Formulation of the probabilistic model according to the problem.The constructed model should be consistent with the actual problem or system in terms of the main characteristic parameters.This must generate a probability distribution function, and this distribution is described by its PDF or CDF.The formulation requires knowing the CDF.
Then, repeat the following process, as many times as needed to have a representative result.
Step 2 Generation of a random number in the interval (0,1) according to a uniform distribution.
Step 3 Use of the inverse of CDF to obtain the value in the desired distribution.Figure 10 illustrates the results of the cumulative distribution of output power of the whole WF on different slopes (0 for a and b, 0.05 for c and d, and 0.12 for e and f) using the proposed analytical method under different sample sizes (100 for a, c and e; and 50,000 for b, d, and f).The curves in Figure 10 correspond to "calculated power including wake effects".
parameters.This must generate a probability distribution function, and this distribution is described by its PDF or CDF.The formulation requires knowing the CDF.Then, repeat the following process, as many times as needed to have a representative result.
Step 2 Generation of a random number in the interval (0,1) according to a uniform distribution.
Step 3 Use of the inverse of CDF to obtain the value in the desired distribution.Figure 10 illustrates the results of the cumulative distribution of output power of the whole WF on different slopes (0 for a and b, 0.05 for c and d, and 0.12 for e and f) using the proposed analytical method under different sample sizes (100 for a, c and e; and 50,000 for b, d, and f).The curves in Figure 10 correspond to "calculated power including wake effects".From Figure 10, it is revealed that with the growth of sample size, particularly when the sample size is set to be 50,000, the results from the analytical model and the Monte Carlo simulation are almost identical.

Figure 3 .
Figure 3. Illustration of WTs' wake on a flat area.

Figure 3 .
Figure 3. Illustration of WTs' wake on a flat area.

Figure 4 .
Figure 4. Illustration of WTs' wake over a hill with gentle gradient.

Figure 4 .
Figure 4. Illustration of WTs' wake over a hill with gentle gradient.

Figure 5 .
Figure 5. Illustration of WTs' wake over a hill with moderate gradient.

Figure 5 .
Figure 5. Illustration of WTs' wake over a hill with moderate gradient.

Figure 6 .
Figure 6.Two scenarios of wake overlap for a hilly WF.

Figure 6 .
Figure 6.Two scenarios of wake overlap for a hilly WF.

Figure 7 .
Figure 7. Illustration of WTs' wake over a hill with steep gradient.

Figure 7 .
Figure 7. Illustration of WTs' wake over a hill with steep gradient.

Figure 8 .
Figure 8. Flowchart of the implementation of the power calculation function.

Figure 9 .
Figure 9. Simulation results of the WF power curve functions with different slopes.

Table 1 .
Parameters of the simulation cases.

Table 1 .
Parameters of the simulation cases.
Simulation results of the WF power curve functions with different slopes.