Actuator Disc Approach of Wind Turbine Wake Simulation Considering Balance of Turbulence Kinetic Energy

The operation of the wind turbines downstream is affected by the wake of the wind turbines upstream. Wind turbine wake flow is investigated by applying the actuator disc (AD) method. The modified k-ε turbulence model is proposed by using both the turbulent kinetic energy source term and the dissipation rate source term to improve the standard k-ε turbulence model for coordinating the generation and the dissipation of the turbulent kinetic energy. The dissipation rate parameter C4ε that obeys a parabolic distribution is used, based on theoretical analysis. The force distributed on the AD is also used instead of a constant, as used in the classical AD method. The simulation results were consistent with the measurements that correspond to different kinds of wind turbines and conditions. The nacelle and the inflow turbulence intensity have great influences on accurately simulating the wake, so it is necessary to imitate the rotor along with the nacelle and accurately measure the inflow turbulence intensity.


Introduction
The wake region forms behind the rotor when air flowing through the running wind turbine has a major effect on the safe operation of the wind turbines downstream.In particular, the variation of the turbulence intensity and the wind speed in this region are important factors to consider.The power loss of a large offshore wind farm caused by the wake effect can be as high as 10-20% [1], so, the wake is an important factor that affects the safe operation of the wind turbine and benefits of the wind farm.
In the early stage of wake effect research, wind tunnel experiments and field measurements were used as the main experimental methods.An early wind tunnel experimental study was carried out at the National Aeronautical Research Institute (Swedish: Flygtekniska Försöksanstalten Abbreviated; FFA) and there were several more wind tunnel tests in recent years [2].Some important onshore field measurements can be found in Tjaereborg [3], Nibe [4] and Sexbierum [5].The most famous offshore wind farm measurements are from Horns Rev [6] and Nysted.With the advancement of technology, advanced measurement methods are constantly appearing and the measurements of wind turbine wakes can be obtained by various methods, such as the supervisory control and data acquisition (SCADA) system, sodar and lidar [7][8][9][10].But not all the data are public.It is difficult to solve the problem with the experimental method only, because it requires huge amounts of money and time but is unable to provide flow details, although the experimental method is credible.
With the rapid development of computer technology, this problem has been solved by the computational fluid dynamics (CFD) method.Simulation methods of the wind turbine wakes have been continually developing [11].The engineering "wake model" was established earlier and has a simple structure and short calculation time.The most famous Jensen model, which is also the basic model of the commonly used commercial wind resources software WAsP, was proposed to evaluate the wind turbine wakes on flat terrain [12].The eddy viscosity model appeared as another analytical wind turbine wakes model and was used in the commercial wind resources software WindFarmer belonging to the Garrad Hassan company [13].These models are based on assumptions of the flow behavior, such as axisymmetry and self-similarity.They work fast and are suitable for engineering applications, but they cannot predict wake accurately.The blade model can be constructed directly in the flow field and used to precisely calculate the load of the wind turbine and other aspects related to wind turbine performance [14].Qian Y.R. used this method to predict the air-load with different wind speeds under yaw conditions [15].However, this method needs more grids around the blade and a lot of computing resources; therefore, it does not apply to multiple wind turbines or wind farm scales.The actuator line model simplifies the blade into a line which is rotated as the blade, then constructs a functional relationship between the aerodynamic data of the blade's local section airfoil and body force, so that it can simulate the interaction between the blade and the inflow [16].The actuator line model, combined with the large eddy simulation (LES) method, was used to investigate the near wake region and explain the reason why the velocity deficit and turbulence intensity increase nearby the rotor [17].Storey R.C. coupled the model with multi-physics method in order to research the wind turbine performance under transient conditions [18].The actuator line model also can demonstrate the wake of the wind turbine over complex terrain and reduce the cost of wind farm simulation to a certain extent [19], but it still takes a lot of computing resources.The AD method can solve this problem by equating the rotor to a penetrable disc to greatly reduce the calculation time [20].
Wind turbines work in the atmospheric boundary layer (ABL).So, not only the wind turbine but also the flow field in the ABL should be simulated well to improve the numerical accuracy of the wake.To some extent, the flow field can be precisely predicted in the ABL by using the LES method.Nilsson K. simulated the power production of the Lillgrund wind farm using the LES method [21].The wake of the horizontal and vertical axis wind turbine was imitated by Martínez-Tossas L.A. and Shamsoddin S. also using the LES method [22,23].However, this takes a long time and requires more computing resources.The Reynolds Averaged Navier Stokes (RANS) method can balance the numerical accuracy of the wake simulation and the need for time and computing resources.
The standard k-ε turbulence model is widely used in Computing Wind Engineering.In 2009, Europe organized a blind comparison based on the Bolund Island Wind Measurement Project which used high frequency data from 35 anemometers and provided a unique database designed to validate micro-scale flow models [24].With its great influence, more than 40 groups from around the world participated in the blind comparison with over 50 model predictions [25,26].The results showed that the CFD method for solving the RANS equation had as smaller calculation error and was superior to other model results, thus, it was obviously superior to the traditional linear model.In this blind comparison, the ten results of the model with the smallest errors are shown, the k-ε turbulence model accounts for more than half of them and showed higher accuracy [27], so that the k-ε turbulence model was found better for simulating in the ABL.
In summary, the wind turbine wakes can be studied by the AD method combined with the k-ε turbulence model.The combined method has been used by the commercial codes WindSim and Fuga as their basic CFD wake model [28][29][30].The flat topographic wind resources have almost been exhausted, wind power projects in complex terrain and offshore have become the mainstream of development.The flow field of wind turbines was simulated by using the combined CFD and AD method in complex terrain [31].Along with the CFD method, experiments were taken to investigate the performance of the wind farm operation in complex terrain [32].The interaction of multiple wind turbine wakes has also been studied [33].Three k-ε turbulence models were used to evaluate the wind energy for a highly complex terrain by Dhunny A.Z. [34].The interaction of atmospheric flow, wind turbine wakes and ground surface has been taken into consideration under different situations [35].The wakes and power of the offshore wind farm can also be investigated by using the combined CFD and AD method [36].The standard k-ε turbulence model overestimates the wake recovery, so a k-ε model which include a parameter sensitive to high velocity gradients was established to improve the simulation accuracy [37,38].The purpose of adding the dissipation rate source term into the k-ε turbulence model by Kasmi A.E. was to solve this problem [39].The improved model is also used for optimizing the layout of the wind farm [40].
In order to investigate the wake accurately, a modified k-ε turbulence model is proposed in which both the turbulent kinetic energy source term and the dissipation rate source term are added.The major improvement is that it coordinates the generation and dissipation of turbulent kinetic energy better than only adding the dissipation rate source term.The parameter C 4ε , which obeys a parabolic distribution, is used based on a theoretical analysis.The body force distribution on the AD is also used, instead of a constant value which is used in the classical AD method.The results are consistent, with the measurements, and better reflect the relative velocity distribution of the wake.

Actuator Disc Model
The wind turbine is represented by a simplified AD method; it equates the wind turbine to a penetrable disc, the air is free to pass through, but encounters some drag [41].To some extent, the effect of the rotor on the inflow can be reflected by this hypothetical model.In this paper, the AD method is derived from one-dimensional momentum theory (Figure 1).model which include a parameter sensitive to high velocity gradients was established to improve the simulation accuracy [37,38].The purpose of adding the dissipation rate source term into the k-ε turbulence model by Kasmi A.E. was to solve this problem [39].The improved model is also used for optimizing the layout of the wind farm [40].
In order to investigate the wake accurately, a modified k-ε turbulence model is proposed in which both the turbulent kinetic energy source term and the dissipation rate source term are added.The major improvement is that it coordinates the generation and dissipation of turbulent kinetic energy better than only adding the dissipation rate source term.The parameter C4ε, which obeys a parabolic distribution, is used based on a theoretical analysis.The body force distribution on the AD is also used, instead of a constant value which is used in the classical AD method.The results are consistent, with the measurements, and better reflect the relative velocity distribution of the wake.

Actuator Disc Model
The wind turbine is represented by a simplified AD method; it equates the wind turbine to a penetrable disc, the air is free to pass through, but encounters some drag [41].To some extent, the effect of the rotor on the inflow can be reflected by this hypothetical model.In this paper, the AD method is derived from one-dimensional momentum theory (Figure 1).It is assumed that the inflow velocity is u1, and the axial thrust of the rotor is: where ρ is the air density and CT is the thrust coefficient.Therefore, the thrust of the AD in equations is: where Su is the source term of the rotor, ∆x is the thickness of the AD.
The thrust is mainly determined by the thrust coefficient and inflow velocity.It is easy to obtain the inflow velocity with a single wind turbine.However, it becomes difficult to determine in the case of multiple wind turbines, so changes need to be made as follows: ( ) where a is the induction factor, uD is the wind speed at the AD [42].The Equation ( 3) is not consistent with the experiment when the value of a is too large, so Glauert's correction and Shen's correction are popular countermeasures [43,44].The Equation (3) will have non-physical rationality when a > 0.5, in the present model, the parameter CT represents the actual performance of the rotor, so the maximum value of a calculated by CT will not exceed 0.5.It is assumed that the inflow velocity is u 1 , and the axial thrust of the rotor is: where ρ is the air density and C T is the thrust coefficient.Therefore, the thrust of the AD in equations is: where S u is the source term of the rotor, ∆x is the thickness of the AD.The thrust is mainly determined by the thrust coefficient and inflow velocity.It is easy to obtain the inflow velocity with a single wind turbine.However, it becomes difficult to determine in the case of multiple wind turbines, so changes need to be made as follows: where a is the induction factor, u D is the wind speed at the AD [42].The Equation ( 3) is not consistent with the experiment when the value of a is too large, so Glauert's correction and Shen's correction are popular countermeasures [43,44].The Equation (3) will have non-physical rationality when a > 0.5, in the present model, the parameter C T represents the actual performance of the rotor, so the maximum value of a calculated by C T will not exceed 0.5.The momentum source term of the rotor is changed as follows: where C x = 4a/(1−a) is the drag coefficient of the AD.So, the momentum source term is linked to the wind speed on the AD and it is more suitable for predicting the flow.
The nacelle can also be simplified as a part of the AD and its drag to the incoming flow can be calculated as the flow drag of the solid.The thrust on the flow is: where C D is the drag coefficient of the nacelle, which varies from 0.8 to 1.2 [45].In this paper, C D takes a value of 1.0.

Governing Equations
The steady, incompressible RANS equations are used: where ρu i u j is the Reynolds stress, S u and S d are the source terms which represent the rotor and the nacelle respectively.They are valid only in the specified area of the AD zone.
The k-ε turbulence model was selected, and it is briefly introduced below: where is the production of kinetic energy; µ τ is the turbulent viscosity coefficient which is used to calculate the fluid viscosity and then solve the governing equations; S k is the turbulent kinetic energy source term and S ε is the dissipation rate source term.They are valid only in the vicinity of the AD.The following model constants are proposed for ABL flows: C µ = 0.033, C 1ε = 1.176 and C 2ε = 1.92 are constants, σ k = 1.0 is the Prandtl number for k and σ ε = 1.3 is the Prandtl number for ε [39].

Turbulence Model Source Term
The flow field due to the effect of the rotor is complicated.Its features are a violent change in wind speed and variation in the production and dissipation of the turbulent kinetic energy.These variations mainly occurred in the vicinity of the AD, where the turbulent kinetic energy source term is added to the k-ε model (Equation ( 9)), which represents the rate of kinetic energy transformation into turbulent kinetic energy.Just like the air flow through the tree canopy, the tree elements make the inflow into wake turbulence, so the kinetic energy converted to turbulent kinetic energy [46].The present model does not consider the thrust distribution on the rotor, the thrust distribution discussed in Section 3.3.So, the equation of the turbulent kinetic energy source term is as follows: where β p is the coefficient that represents the transformation from kinetic energy into turbulent kinetic energy, and the value is usually around 1.0 [46,47].β d = 1.0 is the coefficient of the loss of turbulent kinetic energy [48].
Differing from a canopy with large density, there are only three blades on a wind turbine; the influence on the inflow is much smaller than for the canopy and the rate of kinetic energy converted to turbulent kinetic energy is also smaller.So β p is calculated using the following equation: where β 0 = 0.1.Due to such modifications, the parameter β p can take into account the change of the induction factor a (change according to the thrust of the wind turbine), and its value range is smaller than that set for canopy-flow.In this paper, the vicinity of the AD was set to a cylinder that is coaxial with the rotor, its diameter was the same as that of the rotor, and the thickness was 0.5D (0.25D, respectively, before and after the AD) according to reference [39].
The dissipation is also increased in the vicinity of the AD, because the mixing length scale becomes shorter due to the disturbance of the flow field.S ε in the k-ε model (Equation (10)) represent the dissipation rate source term.It represents the rate of the large-scale turbulence transformation into small-scale turbulence and then dissipates as thermal energy.The equation is as follows: where C 4ε is usually taken as 0.37 [39] and the effect of changing C 4ε on the calculation results is given for a comparison with the measurements (see Section 4.3.1.).The additional term S ε was first used by Chen Y.S. and the parameter C 4ε was taken to be 0.25 [49].Kasmi A.E. also used Equation ( 14), but set the parameter to be 0.37 for calculating the wind turbine wakes.Since then, this value (C 4ε = 0.37) has begun to use for wind turbine simulation in ABL [50] or thermally-stratified ABL [51], the wake of tidal turbine was also studied by using this constant [52].However, the influence of the blade on the inflow is not uniform, owing to the blade which is combined with different shapes of airfoil.The actual dissipation rate should have a certain distribution along the radial direction.Linear distribution of the parameter C 4ε was used instead of constant value by Xu C. [53].The mixed length scale becomes shorter in the vicinity of the tip and root of the blade, because of the disturbance due to the vortices that fall off the tip and root of the blade.The shorter mixed length scale leads to an increase in the dissipation, so there is a distribution which with a higher dissipation on both tip and root of the blade, but lower on the middle of the blade.This distribution happens on all kinds of blades, therefore, in this paper, the corrected C 4ε that obeys Energies 2019, 12, 16 6 of 19 a parabolic distribution was used, as it is suitable for all kinds of wind turbine.The equation is as follows: where r is the distance between the cell center and the axis of the AD, R is the radius of the rotor.

Momentum Source Correction
The force is evenly distributed throughout the disc in the classical AD method, so that the mutation of the velocity happens at the edge of the AD.However the airfoil and angle of attack of the blade along the radial direction are different, and the abilities to absorb wind energy are also different; therefore, the force should have a certain distribution on the AD.A corrective method is introduced into the AD method in some papers.In this paper, a radial distribution function was used to amend the distribution of the source term that represents the rotor [53]: where η is the correction factor, and its equation is as follows: where b = r/R, D 0 = − 4 3 and D 1 = − 2 9 are the coefficients which are given by the user to make the force distribution more realistic [53].

Computational Domain and Grid
The AD located at 3D from the entrance.Figure 2 shows the sections.The dimensions were set to 33D, 7D, and 7D in the streamwise, transverse, and vertical directions.introduced into the AD method in some papers.In this paper, a radial distribution function was used to amend the distribution of the source term that represents the rotor [53]: where η is the correction factor, and its equation is as follows: where b = r/R,

Computational Domain and Grid
The AD located at 3D from the entrance.Figure 2 shows the sections.The dimensions were set to 33D, 7D, and 7D in the streamwise, transverse, and vertical directions.The resolution at the AD was 0.025D, the resolution along the upstream and downstream regions of the AD increased by 1.1 times in every layer.The uniform resolution in the transverse direction was 0.14D.The resolution of the first layer grid on the ground was 0.025D, the resolution increased along the vertical by 1.05 times in every layer.The grid was concentrated at the AD and the near ground layer, in order to ensure accurate calculations.The grid was a regular hexahedron structured mesh with high quality.The total number of grids was 1.7 million.

Boundary Conditions
The neutral ABL model was used in this paper [54].The inlet boundary conditions were defined as follows: 3 The resolution at the AD was 0.025D, the resolution along the upstream and downstream regions of the AD increased by 1.1 times in every layer.The uniform resolution in the transverse direction was 0.14D.The resolution of the first layer grid on the ground was 0.025D, the resolution increased along the vertical by 1.05 times in every layer.The grid was concentrated at the AD and the near ground layer, in order to ensure accurate calculations.The grid was a regular hexahedron structured mesh with high quality.The total number of grids was 1.7 million.

Boundary Conditions
The neutral ABL model was used in this paper [54].The inlet boundary conditions were defined as follows: where u * is the turbulent friction velocity, κ is the von Karman universal constant, z is the vertical height above ground, and z 0 is the roughness length.Experimental data related to velocity and turbulence intensity are generally only provided at the hub height at the entrance.Therefore, the relationship between turbulence intensity and turbulent kinetic energy is given by: where λ is the constant and was set to 1.5 with an isotropic turbulence model [55], k is the turbulent kinetic energy and U is the velocity in the flow field.
The inlet boundary conditions were also used at the top of the computational domain, in order to avoid air flow out of the top.The symmetric boundary condition was used at both sides of the computational domain.The standard wall function was used to consider the effect of the roughness, and its relationship is given by: where C s is the roughness constant and was set to 1.0 [54].

Results and Discussion
The independence study of the grid was done to get the optimal resolution.Table 1 shows the different numbers of cells on the AD.The results of the mesh independence study in the downstream section (2.5D) are shown in Figure 3.The relative velocity distribution keeps changing as the numbers of the grids increase, and then achieves a steady state.In order to improve the accuracy of computation and save the computational resources, the parameter is used as in Test 3.
The accuracy and adaptability of the proposed model in this section is verified against different types of wind turbine and inflow conditions.Three sets of measurements are selected for comparison with the results of the proposed model: (a) the Nibe-B wind turbine has a rotor diameter of 40 m and its hub height is 45 m [39]; (b) the Danwin wind turbine has a rotor diameter of 31 m and its hub height is 23 m [39]; and (c) the detailed analysis of turbine wake has a rotor diameter of 43.2 m and its hub height is 50 m.The measurements were carried out by Garrad Hassan at the Marchwood Engineering Laboratories (MEL) ABL wind tunnel [56].The results of the mesh independence study in the downstream section (2.5D) are shown in Figure 3.The relative velocity distribution keeps changing as the numbers of the grids increase, and then achieves a steady state.In order to improve the accuracy of computation and save the computational resources, the parameter is used as in Test 3. The accuracy and adaptability of the proposed model in this section is verified against different types of wind turbine and inflow conditions.Three sets of measurements are selected for comparison with the results of the proposed model: (a) the Nibe-B wind turbine has a rotor diameter of 40 m and its hub height is 45 m [39]; (b) the Danwin wind turbine has a rotor diameter of 31 m and its hub height is 23 m [39]; and (c) the detailed analysis of turbine wake has a rotor diameter of 43.2 m and its hub height is 50 m.The measurements were carried out by Garrad Hassan at the Marchwood Engineering Laboratories (MEL) ABL wind tunnel [56].

Simulation of the Nibe-B Wind Turbine
The Nibe wind farm is located in the northern coastal region of Denmark, which consists of two wind turbines with a space of 200 m along the north-south line; the Nibe-B wind turbine is one of them.There are four masts in the wind farm, which are arranged coaxially with the two wind turbines.On the west side of the wind turbine is a shallow water area.The area is broad and flat on the east side, and the surface is covered with grass.Therefore, the airflow of the wind farm can be considered to be a basically unaffected uniform flow.The rated wind speed is 13 m/s and the speed of cut-in and cut-out is 6 m/s and 25 m/s respectively.Similar to the proposed model, an improved k-ω turbulence model was also developed for simulating wakes of the Nibe-B wind turbine in neutral ABL [57].The measurements of three operating conditions were compared (i) U = 8.54 m/s, C T = 0.82, I = 11%; (ii) U = 9.56 m/s, C T = 0.77, I = 11%; and (iii) U = 11.52 m/s, C T = 0.67, I = 10.5%.
Figure 4a shows that the variation of the parameter C 4ε has a great influence on the relative velocity of the wake's centerline [53].The wake recovery obviously reduces with an increase of parameter C 4ε , and the minimum wind speed is around 2.5D downstream of the AD.The momentum in the wake is not well replenished by the external fluid.The wake recovery is slow because the momentum exchange in the wake region is weak.Therefore, the predicted velocity at the 6D and 7.5D positions are less than the measurements.A great disturbance to the flow field is caused by the disturbance of fall-off vortices both, from the tip and the root of the blade, making the mixing length shorter at the positions of the blade root and tip, and increasing the dissipation rate.Therefore, taking the parameter C 4ε as a constant reduces the prediction accuracy, and the radial distribution of C 4ε must be used as shown in Equation (15).The big change of the relative velocity happened with the parameter C 4ε varying from 0.1 to 0.9 (nine times growth).The relative velocity is sensitive to the parameter C 4ε , when it takes values within the same order of magnitude as the value variation according to Equation (15).
Figure 4b shows that the variation of parameter β 0 has an influence on the relative velocity [53].Parameter β 0 is used in Equation (13).The wake recovery is almost unchanged, as parameter β 0 varies from 0.025 to 0.1 (four times growth) and minor changes happen as the parameter varies from 0.1 to 1.0 (ten times growth).Big changes happen as the parameter varies from 1.0 to 5.0 (five times growth).This means that only the order of magnitude change of the parameter will obviously change the relative velocity.In contrast to the sensitivity of the relative velocity to the changes in parameter C 4ε , the relative velocity is not obviously sensitive to parameter β 0 , when β 0 with values taken within the same order of magnitude as the value is set as 0.1 in this paper.The relative velocity is also sensitive to the changes in the inflow turbulence intensity and discussed in Section 5.2.
shorter at the positions of the blade root and tip, and increasing the dissipation rate.Therefore, taking the parameter C4ε as a constant reduces the prediction accuracy, and the radial distribution of C4ε must be used as shown in Equation (15).The big change of the relative velocity happened with the parameter C4ε varying from 0.1 to 0.9 (nine times growth).The relative velocity is sensitive to the parameter C4ε, when it takes values within the same order of magnitude as the value variation according to Equation (15).Figure 4b shows that the variation of parameter β0 has an influence on the relative velocity [53].Parameter β0 is used in Equation (13).The wake recovery is almost unchanged, as parameter β0 varies from 0.025 to 0.1 (four times growth) and minor changes happen as the parameter varies from 0.1 to 1.0 (ten times growth).Big changes happen as the parameter varies from 1.0 to 5.0 (five times growth).This means that only the order of magnitude change of the parameter will obviously change the relative velocity.In contrast to the sensitivity of the relative velocity to the changes in parameter C4ε, the relative velocity is not obviously sensitive to parameter β0, when β0 with values taken within the same order of magnitude as the value is set as 0.1 in this paper.The relative velocity is also sensitive to the changes in the inflow turbulence intensity and discussed in Section 5.2.
Comparisons between the proposed model, other models and the measurements in three inflow conditions are shown in Figures 5-7 (four downstream sections which are perpendicular to the AD axis).The results calculated by Kasmi A.E. are also shown in Figures 5-7 [39].The parameter C4ε is set to a constant by Kasmi A.E.The results of the LES method are also shown in Figure 5a,b,d [37].Comparisons between the proposed model, other models and the measurements in three inflow conditions are shown in Figures 5-7 (four downstream sections which are perpendicular to the AD axis).The results calculated by Kasmi A.E. are also shown in Figures 5-7 [39].The parameter C 4ε is set to a constant by Kasmi A.E.The results of the LES method are also shown in Figure 5a,b,d [37].Currently, the power generation is calculated by using the velocity at the hub position, and the error analysis of the velocity at the hub position based on measurements is done and shown in Tables 2 and 3.
Currently, the power generation is calculated by using the velocity at the hub position, and the error analysis of the velocity at the hub position based on measurements is done and shown in Tables 2 and 3.  [53], Laan [37] and Kasmi [39]).
Table 2. Error analysis of the relative velocity of the Nibe-B wind turbine based on the results of the proposed model, the LES method and the standard k-ε turbulence model (EXP date and LES date obtained from Xu [53] and Laan [37]).[53], Laan [37] and Kasmi [39]).The results of the LES method are much closer to the measurements (average error is 4%), although they have slight differences (errors of 3% and −9%) in Figure 5a,b.The results of the proposed model are second only to the results of the LES method (average error is 9%), better than the results calculated by Kasmi A.E. and the standard k-ε model, due to adding the source terms.The added source terms make it possible to coordinate the generation and dissipation of turbulent kinetic energy better than the standard k-ε turbulence model.Therefore, it is necessary to correct the standard k-ε turbulence model.
The relative velocity in the section center coincides well with the measurements at position 2.5D, although it has a slight difference in Figure 7a (error is 7%).The flow field in this near wake region is so complicated that it is hard to simulate precisely.The largest differences appeared in Figure 5b; both the LES method and the proposed model all underestimated the velocity recovery, but the results of the LES method (error is −9%) are better than the results of the proposed model (error is −17%).Differences between the results and the measurements in the section center appear in Figures 5d and 6c.This is because the lower turbulent kinetic energy reduces the momentum exchange, so the wake recovery is slower.A certain difference of the relative velocity distribution along the radial direction is also illustrated in Figures 5b, 6b and 7b.This is because the turbulence intensity calculated by the turbulence model along the radial direction is higher than the actual situation, resulting in a larger momentum exchange in this area; therefore, the velocity recovery is faster than the measurements, and the error is acceptable in a certain range (average error is less than 10%).
The wake recovery gets faster with decrease of the thrust coefficient caused by variation in the inflow velocity.For example, at position 2.5D, when U = 8.54 m/s, the relative velocity in the wake center is 0.46; when U = 9.56 m/s, the relative velocity in the wake center is 0.49; while when U = 11.52 m/s, it is 0.53.The tendency is also captured by the results of the proposed model.
The effect on inflow velocity is also different along the radial direction, due to the change of the blade geometry along that direction.Therefore, in Figure 8, the turbulence intensity increases first from the center, and then changes along the radial direction, and finally reaches the environmental level.The measurements show a maximum turbulence intensity of 0.194 at the 0.8R position.The highest turbulence intensity forms at this position, because the vortices formed at the tip of the blade are "shed" downstream, then broken, therefore leading to an increase in the turbulence intensity.
The results of the proposed model have a maximum turbulence intensity of 0.205 at the 0.6R position, that is, only a position deviation of 0.2R and a value deviation of 0.011 from the measurements.This is also the reason for a certain difference between the relative velocity distribution and the measurements along the radial direction Figures 5b, 6b and 7b.The value of the maximum turbulence intensity is close to the measurements.
both the LES method and the proposed model all underestimated the velocity recovery, but the results of the LES method (error is −9%) are better than the results of the proposed model (error is −17%).Differences between the results and the measurements in the section center appear in Figures 5d and 6c.This is because the lower turbulent kinetic energy reduces the momentum exchange, so the wake recovery is slower.A certain difference of the relative velocity distribution along the radial direction is also illustrated in Figures 5b, 6b and 7b.This is because the turbulence intensity calculated by the turbulence model along the radial direction is higher than the actual situation, resulting in a larger momentum exchange in this area; therefore, the velocity recovery is faster than the measurements, and the error is acceptable in a certain range (average error is less than 10%).[53] and Kasmi [39]).[53] and Kasmi [39]).
The wake recovery gets faster with decrease of the thrust coefficient caused by variation in the inflow velocity.For example, at position 2.5D, when U = 8.54 m/s, the relative velocity in the wake center is 0.46; when U = 9.56 m/s, the relative velocity in the wake center is 0.49; while when U = 11.52 m/s, it is 0.53.The tendency is also captured by the results of the proposed model.[53] and Kasmi [39]).
both the LES method and the proposed model all underestimated the velocity recovery, but the results of the LES method (error is −9%) are better than the results of the proposed model (error is −17%).Differences between the results and the measurements in the section center appear in Figures 5d and 6c.This is because the lower turbulent kinetic energy reduces the momentum exchange, so the wake recovery is slower.A certain difference of the relative velocity distribution along the radial direction is also illustrated in Figures 5b, 6b and 7b.This is because the turbulence intensity calculated by the turbulence model along the radial direction is higher than the actual situation, resulting in a larger momentum exchange in this area; therefore, the velocity recovery is faster than the measurements, and the error is acceptable in a certain range (average error is less than 10%).[53] and Kasmi [39]).[53] and Kasmi [39]).
The wake recovery gets faster with decrease of the thrust coefficient caused by variation in the inflow velocity.For example, at position 2.5D, when U = 8.54 m/s, the relative velocity in the wake center is 0.46; when U = 9.56 m/s, the relative velocity in the wake center is 0.49; while when U = 11.52 m/s, it is 0.53.The tendency is also captured by the results of the proposed model.The effect on inflow velocity is also different along the radial direction, due to the change of the blade geometry along that direction.Therefore, in Figure 8, the turbulence intensity increases first from the center, and then changes along the radial direction, and finally reaches the environmental level.The measurements show a maximum turbulence intensity of 0.194 at the 0.8R position.The highest turbulence intensity forms at this position, because the vortices formed at the tip of the blade are "shed" downstream, then broken, therefore leading to an increase in the turbulence intensity.
The results of the proposed model have a maximum turbulence intensity of 0.205 at the 0.6R position, that is, only a position deviation of 0.2R and a value deviation of 0.011 from the measurements.This is also the reason for a certain difference between the relative velocity distribution and the measurements along the radial direction Figures 5b, 6b and 7b.The value of the maximum turbulence intensity is close to the measurements.

Simulation of the Danwin Wind Turbine
The Danwin wind turbine is operated in the Alsvik wind farm, which has four wind turbines, and three masts and is near the west coast of Gotland.The wind farm is in a flat coastal terrain covered with grazed grass [58].The rated wind speed is 12 m/s and the cut-in speed is 5 m/s.The measurements of two operating conditions were compared: (i) U = 8 m/s, CT = 0.82, I = 7%; and (ii) U = 11 m/s, CT = 0.65, I =6%.The results calculated with the full rotor model [14] and by Kasmi A.E. [39] Figure 8. Wake turbulence intensity distribution of the Nibe-B wind turbine corresponding to U = 11.52 m/s at the downstream section of x = 2.5D (EXP date obtained from Xu [53]).

Simulation of the Danwin Wind Turbine
The Danwin wind turbine is operated in the Alsvik wind farm, which has four wind turbines, and three masts and is near the west coast of Gotland.The wind farm is in a flat coastal terrain covered with grazed grass [58].The rated wind speed is 12 m/s and the cut-in speed is 5 m/s.The measurements of two operating conditions were compared: (i) U = 8 m/s, C T = 0.82, I = 7%; and (ii) U = 11 m/s, C T = 0.65, I =6%.The results calculated with the full rotor model [14] and by Kasmi A.E. [39] are shown in Figures 9 and 10.The full rotor model needs to construct more grids around the blade and takes a long time to calculate.The error analysis of velocity at the hub position was also done and is shown in Table 4.

Simulation of the Danwin Wind Turbine
The Danwin wind turbine is operated in the Alsvik wind farm, which has four wind turbines, and three masts and is near the west coast of Gotland.The wind farm is in a flat coastal terrain covered with grazed grass [58].The rated wind speed is 12 m/s and the cut-in speed is 5 m/s.The measurements of two operating conditions were compared: (i) U = 8 m/s, CT = 0.82, I = 7%; and (ii) U = 11 m/s, CT = 0.65, I =6%.The results calculated with the full rotor model [14] and by Kasmi A.E. [39] are shown in Figures 9 and 10.The full rotor model needs to construct more grids around the blade and takes a long time to calculate.The error analysis of velocity at the hub position was also done and is shown in Table 4.  [53], Kasmi [39] and Abdelsalam [14]).m/s at the downstream position of x = 6.2D (EXP date, Kasmi date and full-rotor date obtained from Xu [53], Kasmi [39] and Abdelsalam [14]).Table 4. Error analysis of the relative velocity of the Danwin wind turbine based on the results of the proposed model, the standard k-ε turbulence model, the full rotor model and the results of Kasmi A.E (full-rotor date and Kasmi date obtained from Abdelsalam [14] and Kasmi [39]).Figures 9 and 10 show that the results of the proposed model are more consistent with the measurements than that of Kasmi A.E., and are also similar to the results calculated by the full rotor model.The results of the proposed model are accurate with the two conditions.The disturbance at the 1D position in the wake is strong, as can also be seen from the scattered measured data.Kasmi [39] and Abdelsalam [14]).Table 4. Error analysis of the relative velocity of the Danwin wind turbine based on the results of the proposed model, the standard k-ε turbulence model, the full rotor model and the results of Kasmi A.E (full-rotor date and Kasmi date obtained from Abdelsalam [14] and Kasmi [39]).Figures 9 and 10 show that the results of the proposed model are more consistent with the measurements than that of Kasmi A.E., and are also similar to the results calculated by the full rotor model.The results of the proposed model are accurate with the two conditions.The disturbance at the 1D position in the wake is strong, as can also be seen from the scattered measured data.

Experimental Data by Garrad Hassan
Garrad Hassan published a research report of the wake at the Marchwood Engineering Laboratories (MEL) ABL wind tunnel [56] in 1989.The horizontal axis wind turbine which operated at the hub height of 50 m with a 43.2 m diameter were used.
The report shows the distribution of relative velocity for several downstream sections for (i) U = 5.3 m/s, C T = 0.607 and z 0 = 0.075 m; and (ii) U = 5.3 m/s, C T = 0.858 and z 0 = 0.075 m.The Park results are also shown in the report [56].Figures 11 and 12 show that the results of the proposed model compare well with the measurements and better reflect the relative velocity distribution.Table 5 shows that the results of the standard k-ε turbulence model and the Park model have big differences to the measurements.

Influence of Nacelle
As an important component of the wind turbine, the nacelle is bulky, because it contains heavy power generation equipment.Some researchers ignore the influence of the nacelle in the numerical simulation of the wake.A comparison study was performed to study the impact on the wake with or without the nacelle.In the report [56], the measurements of the relative vertical profiles in three sections (downstream positions at 2.5D, 5D, and 7.5D) are given, corresponding to Tip Speed Rate (TSR) = 4.0, U = 5.3 m/s and C T = 0.793.The relative vertical profile of the ABL are also supplemented in order to observe the wake deficit.
As shown in Figure 13, the wake without the nacelle resulted in smaller vertical velocity loss in 2.5D sections.The deviation between the two cases decreased rapidly as the distance increased.In the 5D section, the velocity of the wake center without the nacelle deviated from the measurements by nearly 5%.The effect of the nacelle on relative vertical profile almost disappeared in the far downstream region.This shows that the nacelle has a greater influence on the vertical velocity profile in the near wake region and has almost no effect in the far wake region.
Energies 2018, 11, x FOR PEER REVIEW 14 of 19 the 5D section, the velocity of the wake center without the nacelle deviated from the measurements by nearly 5%.The effect of the nacelle on relative vertical profile almost disappeared in the far downstream region.This shows that the nacelle has a greater influence on the vertical velocity profile in the near wake region and has almost no effect in the far wake region.Two more operating conditions: (i) TSR = 2.9, U = 5.3 m/s, CT = 0.607 and (ii) TSR = 5.1, U = 5.3 m/s, CT = 0.858 were taken in order to confirm the impact of the nacelle again.As shown in Figures 14 and 15, the nacelle also has a greater influence on the vertical velocity profile in the near wake region but has almost no effect in the far wake region.Two more operating conditions: (i) TSR = 2.9, U = 5.3 m/s, C T = 0.607 and (ii) TSR = 5.1, U = 5.3 m/s, C T = 0.858 were taken in order to confirm the impact of the nacelle again.As shown in Figures 14 and 15, the nacelle also has a greater influence on the vertical velocity profile in the near wake region but has almost no effect in the far wake region.downstream positions as: (a) x = 2.5D (b) x = 5D (c) x = 7.5D (EXP date obtained from Report [56]).
Two more operating conditions: (i) TSR = 2.9, U = 5.3 m/s, CT = 0.607 and (ii) TSR = 5.1, U = 5.3 m/s, CT = 0.858 were taken in order to confirm the impact of the nacelle again.As shown in Figures 14 and 15, the nacelle also has a greater influence on the vertical velocity profile in the near wake region but has almost no effect in the far wake region.

Influence of Turbulence Intensity
The inflow turbulence intensity has a great impact on the flow field of a wind farm.The relative velocity measurements at the hub height in different downstream sections were given in the report corresponding to three conditions: (i) TSR = 2.9, U = 5.3 m/s, CT = 0.607; (ii) TSR = 4.0, U = 5.3 m/s, CT = 0.793; and (iii) TSR = 5.1, U = 5.3 m/s, CT = 0.858.The actual inflow turbulence intensity is 0.12.
Inflow turbulence intensity values of 0.1 and 0.08 were selected to analyze the impact on the flow field.As shown in Figure 16, the overall velocity in the wake region decreased significantly, and the influence at the 25D distance still could not be eliminated.This is because the decrease in the turbulence intensity reduced the momentum exchange, and the momentum was not well replenished by the external fluid, so the wake recovered slowly.Therefore, accurate measurement of the inflow turbulence intensity has a great influence on the accuracy of the simulation of the wake.In the ABL, the actual inflow turbulence intensity varies little.The relative velocity changed obviously with inflow turbulence intensity varying from 0.08 to 0.12 (only 1.5 times growth), so the relative velocity is sensitive to the variation of inflow turbulence intensity.

Influence of Turbulence Intensity
The inflow turbulence intensity has a great impact on the flow field of a wind farm.The relative velocity measurements at the hub height in different downstream sections were given in the report corresponding to three conditions: (i) TSR = 2.9, U = 5.3 m/s, C T = 0.607; (ii) TSR = 4.0, U = 5.3 m/s, C T = 0.793; and (iii) TSR = 5.1, U = 5.3 m/s, C T = 0.858.The actual inflow turbulence intensity is 0.12.
Inflow turbulence intensity values of 0.1 and 0.08 were selected to analyze the impact on the flow field.As shown in Figure 16, the overall velocity in the wake region decreased significantly, and the influence at the 25D distance still could not be eliminated.This is because the decrease in the turbulence intensity reduced the momentum exchange, and the momentum was not well replenished by the external fluid, so the wake recovered slowly.Therefore, accurate measurement of the inflow turbulence intensity has a great influence on the accuracy of the simulation of the wake.In the ABL, the actual inflow turbulence intensity varies little.The relative velocity changed obviously with inflow turbulence intensity varying from 0.08 to 0.12 (only 1.5 times growth), so the relative velocity is sensitive to the variation of inflow turbulence intensity.
turbulence intensity reduced the momentum exchange, and the momentum was not well replenished by the external fluid, so the wake recovered slowly.Therefore, accurate measurement of the inflow turbulence intensity has a great influence on the accuracy of the simulation of the wake.In the ABL, the actual inflow turbulence intensity varies little.The relative velocity changed obviously with inflow turbulence intensity varying from 0.08 to 0.12 (only 1.5 times growth), so the relative velocity is sensitive to the variation of inflow turbulence intensity.

Conclusions
The wake prediction is investigated using a variety of methods in this paper.Three kinds of wind turbine are selected for comparison with the proposed model.The results are consistent with the measurements.Adding the source terms can make the prediction of the wake distribution more accurate by coordinating the generation and dissipation of the turbulent kinetic energy.The correction of the parameter C4ε that obeys a parabolic distribution along the radial direction can make the dissipation rate more realistic in the vicinity of the AD and also improve the predicted accuracy  [56]).

Conclusions
The wake prediction is investigated using a variety of methods in this paper.Three kinds of wind turbine are selected for comparison with the proposed model.The results are consistent with the measurements.Adding the source terms can make the prediction of the wake distribution more accurate by coordinating the generation and dissipation of the turbulent kinetic energy.The correction of the parameter C 4ε that obeys a parabolic distribution along the radial direction can make the dissipation rate more realistic in the vicinity of the AD and also improve the predicted accuracy of the wake.The nacelle and the inflow turbulence intensity have a great influence on the vertical velocity distribution in the wake, especially in the near wake region.Therefore, it is very important to have accurate measurements of the inflow turbulence intensity.Wind turbines work in an unsteady-state condition, of which the steady-state condition is a special case and an important part.Therefore, the steady-state approach is used in this paper and is helpful for understanding and analyzing the unsteady-state calculation.

D
= − are the coefficients which are given by the user to make the force distribution more realistic[53].

Figure 2 .
Figure 2. Computational domain and grid: (a) view in the xz plane; (b) view in the yz plane; (c) view of the grid in the xz plane.

Figure 2 .
Figure 2. Computational domain and grid: (a) view in the xz plane; (b) view in the yz plane; (c) view of the grid in the xz plane.

Figure 3 .
Figure 3.The results of mesh independence study.

Figure 3 .
Figure 3.The results of mesh independence study.

Figure 4 .
Figure 4. Effects of parameters on the relative velocity on the wake's centerline: (a) parameter C4ε; (b) parameter β0 (EXP stands for experimental date obtained from Xu [53]).

Figure 4 .
Figure 4. Effects of parameters on the relative velocity on the wake's centerline: (a) parameter C 4ε ; (b) parameter β 0 (EXP stands for experimental date obtained from Xu [53]).

Figure 8 .
Figure8.Wake turbulence intensity distribution of the Nibe-B wind turbine corresponding to U = 11.52 m/s at the downstream section of x = 2.5D (EXP date obtained from Xu[53]).

Figure 8 .
Figure8.Wake turbulence intensity distribution of the Nibe-B wind turbine corresponding to U = 11.52 m/s at the downstream section of x = 2.5D (EXP date obtained from Xu[53]).

Figure 11 .
Figure 11.Wake relative velocity distribution of the Garrad Hassan wind turbine corresponding to CT = 0.607 at several downstream sections as: (a) x = 2.5D; (b) x = 5D; (c) x = 7.5D (EXP date and Park date obtained from Report[56]).

Figure 12 .
Figure 12.Wake relative velocity distribution of the Garrad Hassan wind turbine corresponding to C T = 0.857 at several downstream sections as: (a) x = 2.5D; (b) x = 5D; (c) x = 7.5D (EXP date and Park date obtained from Report[56]).
Drag coefficient of the actuator disc C µ , C 1ε , C 2ε Turbulence model constants D Diameter of the actuator disc D 0 , D 1 Coefficients of the correction factor η I Turbulence intensity K s Roughness height k Turbulent kinetic energy P k Production of kinetic energy r Distance between the cell center and the axis of the actuator disc R Radius of the rotor S d Source term of the nacelle S ε Dissipation rate source term S k Turbulent kinetic energy source term S u Source term of the actuator disc T Axial thrust of the rotor u 1 Inflow velocity u D Wind speed on the actuator disc u * Turbulent friction velocity z 0 Roughness length β p , β d , β 0 Coefficients in the turbulent kinetic energy source term ∆x Thickness of the actuator disc Prandtl number for k and ε

Table 2 .
[37]r analysis of the relative velocity of the Nibe-B wind turbine based on the results of the proposed model, the LES method and the standard k-ε turbulence model (EXP date and LES date obtained from Xu[53]and Laan[37]).

Table 3 .
[39]r analysis of the relative velocity of the Nibe-B wind turbine based on the results of the proposed model, the standard k-ε turbulence model and the results of Kasmi A.E (EXP date and Kasmi date obtained from Xu[53]and Kasmi[39]).

Table 5 .
[56]r analysis of the relative velocity of the Garrad Hassan (GH) wind turbine based on the results of the proposed model, the standard k-ε turbulence model, and the Park model (Park date obtained from Report[56]).

Table 5 .
[56]r analysis of the relative velocity of the Garrad Hassan (GH) wind turbine based on the results of the proposed model, the standard k-ε turbulence model, and the Park model (Park date obtained from Report[56]).

Table 5 .
[56]r analysis of the relative velocity of the Garrad Hassan (GH) wind turbine based on the results of the proposed model, the standard k-ε turbulence model, and the Park model (Park date obtained from Report[56]).