A Modified Reynolds-Averaged Navier–Stokes-Based Wind Turbine Wake Model Considering Correction Modules

Increasing wind power generation has been introduced into power systems to meet the renewable energy targets in power generation. The output efficiency and output power stability are of great importance for wind turbines to be integrated into power systems. The wake effect influences the power generation efficiency and stability of wind turbines. However, few studies consider comprehensive corrections in an aerodynamic model and a turbulence model, which challenges the calculation accuracy of the velocity field and turbulence field in the wind turbine wake model, thus affecting wind power integration into power systems. To tackle this challenge, this paper proposes a modified Reynolds-averaged Navier–Stokes (MRANS)-based wind turbine wake model to simulate the wake effects. Our main aim is to add correction modules in a 3D aerodynamic model and a shear-stress transport (SST) k-ω turbulence model, which are converted into a volume source term and a Reynolds stress term for the MRANS-based wake model, respectively. A correction module including blade tip loss, hub loss, and attack angle deviation is considered in the 3D aerodynamic model, which is established by blade element momentum aerodynamic theory and an improved Cauchy fuzzy distribution. Meanwhile, another correction module, including a hold source term, regulating parameters and reducing the dissipation term, is added into the SST k-ω turbulence model. Furthermore, a structured hexahedron mesh with variable size is developed to significantly improve computational efficiency and make results smoother. Simulation results of the velocity field and turbulent field with the proposed approach are consistent with the data of real wind turbines, which verifies the effectiveness of the proposed approach. The variation law of the expansion effect and the double-hump effect are also given.


Introduction
In the wake region of a wind turbine, the energy absorption rate of the wind turbine is decreased and the fatigue load is increased due to the decreasing wind speed and increasing intensity of The existing wake models are divided into three categories: semi-empirical engineering models [9,10], vortex method models [11,12] and computational fluid dynamics (CFD) models [13]. The stateof-the-art CFD models can be divided into two classes: Reynolds-averaged Navier-Stokes (RANS) [14] and large-eddy simulation (LES) [15,16]. As the semi-empirical engineering model is onedimensional, the computation is small, and the velocity distribution of the wake region can be located quickly. This model is often used in the software of engineering design, such as Windsim and Windfarmer. However, the calculation results are far from the measured data of wind farms [16]. The vortex method uses a lift line or lift face blade to simulate for the steady wind flow and simple unsteady conditions, especially for the aerodynamic characteristics and wake tip vortex of the wind turbine in the random wind velocity field with the wind shear effect [11]. The CFD method is a general method, which describes the flow field by using hydrodynamic control equations, and it obtains the numerical results of the flow field by using computer program design. The results are affected by the design of RANS and LES models in the CFD method, which needs to improve the accuracy of turbulence models. Moreover, computational cost is also an index to assess the CFD method [17][18][19]. The advantages and disadvantages of research on wind turbine wake models are shown in Table 1 [20][21][22][23]. Table 1. The advantages and disadvantages of research on wind turbine wake models.

Semi-empirical model
Jensen [20] Widely used in industry Suitable for the far wake region Fuga [21] Widely used in industry with medium accuracy Only uses for offshore wind farms Simulation model DNS [22] Direct numerical simulation with high accuracy Requires a lot of computing resources and simulation time LES [23] Widely used in research with high accuracy Appropriate assumptions are necessary DNS (direct numerical simulation), LES (large-eddy simulation).  The existing wake models are divided into three categories: semi-empirical engineering models [9,10], vortex method models [11,12] and computational fluid dynamics (CFD) models [13]. The state-of-the-art CFD models can be divided into two classes: Reynolds-averaged Navier-Stokes (RANS) [14] and large-eddy simulation (LES) [15,16]. As the semi-empirical engineering model is one-dimensional, the computation is small, and the velocity distribution of the wake region can be located quickly. This model is often used in the software of engineering design, such as Windsim and Windfarmer. However, the calculation results are far from the measured data of wind farms [16]. The vortex method uses a lift line or lift face blade to simulate for the steady wind flow and simple unsteady conditions, especially for the aerodynamic characteristics and wake tip vortex of the wind turbine in the random wind velocity field with the wind shear effect [11]. The CFD method is a general method, which describes the flow field by using hydrodynamic control equations, and it obtains the numerical results of the flow field by using computer program design. The results are affected by the design of RANS and LES models in the CFD method, which needs to improve the accuracy of turbulence models. Moreover, computational cost is also an index to assess the CFD method [17][18][19]. The advantages and disadvantages of research on wind turbine wake models are shown in Table 1 [20][21][22][23]. Table 1. The advantages and disadvantages of research on wind turbine wake models.

Semi-empirical model
Jensen [20] Widely used in industry Suitable for the far wake region Fuga [21] Widely used in industry with medium accuracy Only uses for offshore wind farms Simulation model DNS [22] Direct numerical simulation with high accuracy Requires a lot of computing resources and simulation time LES [23] Widely used in research with high accuracy Appropriate assumptions are necessary DNS (direct numerical simulation), LES (large-eddy simulation). Concerning LES, a full-scale 3D wind turbine model by using sliding tetrahedron mesh is established to calculate the wakes in the rotor diameter range of 20 times [24]. The authors of [23] designs a non-full-scale model combined with the model of the wind turbine, and they discuss the function of each sub-grid model under the LES method. Concerning the RANS, a shear-stress transport (SST) k-ω turbulence model, a 2D high-precision computational grid, is used to calculate the wakes [15]. A tree crown model is applied in the turbulence model [13]. Its effect on the wind turbine in regard to incoming airflow is similar to that of tree crown on airflow. Based on the tree crown model, a modified calculation of the turbulence model is proposed by adding a source term [25]. Its calculation results are compared with the measured data of the Danish Nibe B wind turbine. The results show that the modified model is feasible in the near wake region. However, the wake velocity deficit in the far wake region is underestimated. Furthermore, Shen et al. obtain the surface pressure and friction of airfoil by X foil software, and distribute them into the 2D actuator surface model as volume force [26]. The blade surface lift force is calculated by utilizing the vortex method in [27]. In [28] and [29], the scholars develop an aeroelastic fatigue, aerodynamics, structures, and turbulence (FAST) model of a wind turbine and offshore wind farm application (OWFA) model combined with an actuation model considering the influence of blade deformation on wind turbine wakes. An improved k-ε turbulence model is presented for the numerical simulation of wind turbine wakes. However, in the vertical direction, the simulation results do not have good agreement with the experimental results [30].
In summary, the above approaches lack the consideration of comprehensive corrections in the aerodynamic model and turbulence model. Motivated by the above discussion, the velocity field and turbulence field of wind turbines are simulated by a modified RANS (MRANS) method in this paper. The main contributions of this paper are summarized as follows: • This paper proposes a modified Reynolds-averaged Navier-Stokes (MRANS)-based wind turbine wake model to simulate the wake effects. Based on the correction module, the proposed blade element momentum (BEM) -fuzzy aerodynamic model can amend the inconsistent condition between the wake simulation and the experiment test. For the turbulence model, the turbulence attenuation is effectively avoided by adding the hold source term. The accuracy of the turbulence intensity distribution is improved by correcting the closure constant and the dissipation term.

•
Simulation results of the velocity field and turbulent field with the proposed approach are consistent with the data of real wind turbines, which verifies the effectiveness of the proposed approach. Furthermore, the computation efficiency is significantly improved by the developed mesh partition method.
The rest of this paper is organized as follows: Section 2 describes the overall proposed wind turbine wake model. The modified wind turbine 3D model with BEM and fuzzy theory is studied in Section 3. The SST k-ω turbulence model is presented in Section 4. The variable size hexahedron mesh partition method is developed in Section 5. Numerical simulations are given in Section 6. Finally, conclusions are drawn in Section 7.

Overall Wind Turbine Wake Model
In this paper, the RANS model is adopted in the wind turbine wake model. Wind turbine wakes can be regarded as an incompressible flow field, and energy equations are not considered in this paper in accordance with the Reynolds-averaged equations [5]. The continuity equation is given by Equation (1).
where u i is the average velocity component in direction x i , i = 1, 2, 3 . . . The momentum equations are given by Equation (2).
where ρ is the air density, u i is the i-th velocity component, p is positive stress, µ is the dynamic viscosity is the average strain tensor, u i is the pulsation velocity component in direction x i , and f i is the volume source term.
In the process of Reynolds time averaging, an unknown tensor additional term ρu i u j is derived, which is called the Reynolds stress tensor and represents energy transfer caused by turbulence pulsation. Due to this unknown additional term, the Reynolds time-averaged Equations (1) and (2) cannot be closed. The turbulence model can calculate the Reynolds stress tensor by the vorticity viscosity hypothesis. The Boussinesq equation is shown in Equation (3) [5].
where the turbulent kinetic energy k = 0.5u i u j , δ ij = i = j 1 i j 0 , and µ t is the eddy viscosity coefficient determined by the turbulence model. To enhance the accuracy of the abovementioned Reynolds stress tensor and eddy viscosity coefficient, we develop a correction module for an SST k-ω turbulence model. Meanwhile, to further comprehensively complete the volume source term of the RANS model, another correction module is also introduced to the BEM-fuzzy 3D aerodynamic model. By combining the turbulence model and the aerodynamic model, we propose an M-RANS wake model; the overall modeling framework is shown in Figure 2. (2 ) where ρ is the air density, ui is the i-th velocity component, p is positive stress, μ is the dynamic viscosity coefficient, where the turbulent kinetic energy  is the eddy viscosity coefficient determined by the turbulence model. To enhance the accuracy of the abovementioned Reynolds stress tensor and eddy viscosity coefficient, we develop a correction module for an SST k-ω turbulence model. Meanwhile, to further comprehensively complete the volume source term of the RANS model, another correction module is also introduced to the BEM-fuzzy 3D aerodynamic model. By combining the turbulence model and the aerodynamic model, we propose an M-RANS wake model; the overall modeling framework is shown in Figure 2.

BEM-Fuzzy 3D Aerodynamic Model
As shown in Figure 2, the BEM-fuzzy model is established to obtain the volume source term. When compared to the traditional 2D BEM model of the actuator hypothesis, this model is modified by three corrections to the 2D model and by Cauchy fuzzy distribution of the 3D model. Figure 3 illustrates the principle of the BEM theory calculation, which can effectively obtain the lift and resistance of wind turbines [31]. According to BEM theory, the element of local velocity is calculated by Equation (4).
where Ω is the angular velocity of wind wheel rotation, r is the span position of the blade element airfoil section, V ∞ is the axial atmospheric free flow velocity, and a and b are the axial induction factor and tangential induction factor, which are determined by iterative methods [32].

BEM-Fuzzy 3D Aerodynamic Model
As shown in Figure 2, the BEM-fuzzy model is established to obtain the volume source term. When compared to the traditional 2D BEM model of the actuator hypothesis, this model is modified by three corrections to the 2D model and by Cauchy fuzzy distribution of the 3D model. Figure 3 illustrates the principle of the BEM theory calculation, which can effectively obtain the lift and resistance of wind turbines [31]. According to BEM theory, the element of local velocity is calculated by Equation (4).
where  is the angular velocity of wind wheel rotation, r is the span position of the blade element airfoil section,  V is the axial atmospheric free flow velocity, and a and b are the axial induction factor and tangential induction factor, which are determined by iterative methods [32]. The inflow angle  is calculated by Equation (5).
where     is the attack angle, and  is the local pitch angle.
The force vector of the blade element at the span position in the 2D plane is obtained by Equation (6).

Tip Loss Correction
Due to the pressure difference between the lift surface and the pressure surface of the blade, the airflow at the tip and root of the blade will flow twice along the blade radial direction. The moment that acting on the blade is reduced, the blade element force at the tip has a great influence on the aerodynamic performance of the whole blade. Therefore, the tip loss correction is considered. The tip loss coefficient is defined by Equation (7) [33]. The inflow angle φ is calculated by Equation (5).
where α = φ − γ is the attack angle, and γ is the local pitch angle. The force vector of the blade element at the span position in the 2D plane is obtained by Equation (6).
where n b is the number of blades, C L = C L (a, R e ) and C D = C D (a, R e ) (determined by the airfoil aerodynamic characteristic curve) are the lift coefficient and drag coefficient, R e is the Reynolds number of the local chord length c, and e L and e D are the unit vectors for the directions of lift and drag, respectively. Due to the poor precision of the traditional 2D model, it is necessary to add corrections to improve accuracy. The traditional 2D model is modified from the following three aspects.

Tip Loss Correction
Due to the pressure difference between the lift surface and the pressure surface of the blade, the airflow at the tip and root of the blade will flow twice along the blade radial direction. The moment that acting on the blade is reduced, the blade element force at the tip has a great influence on the aerodynamic performance of the whole blade. Therefore, the tip loss correction is considered. The tip loss coefficient is defined by Equation (7) [33].
The recommended values of the empirical parameters c 1 , c 2 , c 3 are 0.125, 21, and 0.1, respectively.

Hub Loss Correction
The hub aerodynamic performance is affected by the separation of airflow to the root of the blade. Similar to the blade tip loss, the hub loss is considered. The hub loss coefficient is defined by Equation (9) [32].
where R hub is the hub radius, and the correction factor η is η = η 1 .η 2 .

Attack Angle Correction
The blade has a certain thickness and width, especially at the root of the blade, which makes the direction of the airflow change greatly. In the front and rear edge of the airfoil, the circumferential velocity of the airflow increases. At the same time, the cross-section area of airflow decreases, and the axial velocity of airflow increases due to the thickness of airfoil. The thickness and width of the blade affect the attack angle. Attack angle changes are given by Equation (10) [34].
where t max is the maximum thickness of blade element airfoil, ∆α 1 is the attack angle change caused by the influence of blade width on the direction of airflow, and ∆α 2 is the attack angle change caused by the influence of blade thickness on airflow direction. The attack angle is revised by Equation (11).
The modified blade element aerodynamic force in the 2D plane is expressed by Equation (12).
Here, we have obtained the 2D model, and it is necessary to establish the 3D model via the improved Cauchy fuzzy distribution.
The thickness of wind turbine volume force distribution in the axis is uncertain. To ensure the stability of numerical simulation and increase the convergence rate of calculation, it is necessary to smooth the volume force to both sides of the neutral layer. Take an observation point of blade airfoil as an example; the volume force of wind turbine firstly remains unchanged in a small range of the axial direction. Subsequently, both ends of the axis are rapidly attenuated to zero. This is essentially a fuzzy problem, as shown in Figure 4. Therefore, the mapping can be obtained by Equation (13).
where the domain U is the volume force at a distance l from the origin, the fuzzy rule A is the stable condition of the domain element, and A(l) is the membership function of domain U. where the domain U is the volume force at a distance l from the origin, the fuzzy rule A is the stable condition of the domain element, and   Al is the membership function of domain U.
(a) Actuator hypothesis model (b) Improved Cauchy fuzzy distribution is the interval element as shown in Figure 4b; the parameters α and β are 1.5 and 3. The cut set The axial distribution is treated as square intervals ΔL in Figure 4 to reduce the computation time.
The calculated thickness L under accuracy  is expressed by Equation (15).
where [ ] is the downward rounding operation.
The attenuation coefficient at any location l within the calculated thickness is obtained by Equation (16).
According to the theory of actuator hypothesis [35], the average lift and drag distribution of the wind turbine in the 2D plane is calculated by Equation (17).
The volume force of the discrete fuzzy distribution along the normal direction of the wind turbine in 3D space is obtained by Equation (18).
The volume forces at three different dimensions are calculated by Equation (19). The membership function A(l) can be determined by the improved Cauchy fuzzy distribution as follows.
where ∆L = sin φ is the interval element as shown in Figure 4b; the parameters α and β are 1.5 and 3.
The cut set A λ = l l ∈ U, A(l) ≥ λ under the accuracy λ level determines the calculation area. The axial distribution is treated as square intervals ∆L in Figure 4 to reduce the computation time.
The calculated thickness L under accuracy λ is expressed by Equation (15).
where [ ] is the downward rounding operation.
The attenuation coefficient at any location l within the calculated thickness is obtained by Equation (16).
According to the theory of actuator hypothesis [35], the average lift and drag distribution of the wind turbine in the 2D plane is calculated by Equation (17).
The volume force of the discrete fuzzy distribution along the normal direction of the wind turbine in 3D space is obtained by Equation (18).
The volume forces at three different dimensions are calculated by Equation (19).
where the original point of the right angle coordinate system is the model center, x-direction is the axial direction, z-direction is the height direction, X 2 = cos φ sin φ sin φ − cos φ is the transformation matrix for a direct coordinate system, is the transformation matrix for the normal force as well as tangential force, and Ψ is the polar angle in a polar coordinate system.

SST k-ω Turbulence Model
As shown in Figure 2, the modified SST k-ω model is established to obtain the Reynold stress term of the global RANS model.
The SST k-ω model includes a k-ω model in the outer layer (near the wall) and a k-ε model in the inner layer (the outer edge of the boundary layer, the free shear layer and the fully developed region of turbulence) [36]. This method combines the k-ω model and the k-ε model by the mixed function of weighted average. The k-ε model has less dependence on the far field condition, and the k-ω model has a high accuracy in the near wall simulation. To improve the simulation accuracy in the strong adverse pressure gradient and separation flow, a modified vortex viscosity coefficient is proposed. The SST k-ω two-equation model is expressed by Equation (20). where where y is the distance from the nearest wall to the calculation point, and the closed constants are set as: β * = 0.09, a 1 = 0.31 α 1 = 5/9, β 1 = 3/40, σ k1 = 0.085, σ ω1 = 0.5 in layer α 2 = 0.44, β 2 = 0.0828, σ k2 = 1, σ ω2 = 0.856 out layer If the two-equation model is used to simulate the free stream, the initial value of the inlet boundary will gradually decrease with the downstream flow. When the fluid approaches to the wind turbine, the local turbulent variable value is no longer the initial value of the inlet boundary, which may cause turbulence attenuation.
The turbulent kinetic energy generation term of the equations is zero due to no velocity gradient in the inlet free flow. The diffusion term and cross diffusion term can be ignored due to no gradient of turbulence variables. Solving Equation (20), we can obtain the result in the X direction as follows: where the index I is the initial inlet boundary, x is the downstream distance, and u is local wind velocity.

Hold Source Term
It can be seen from Equation (22) that the turbulent kinetic energy and the specific dissipation rate will be attenuated from the initial stage of the inlet boundary to the outlet of the computational domain. Therefore, the turbulence attenuation effect must be corrected in the wake calculation. To reduce the turbulence attenuation, the dissipative term caused by turbulence attenuation can be offset by adding the hold source term in the turbulence model. The turbulent equations of the inlet free flow can be obtained by Equation (23).
where k R and ω R are the atmospheric real environmental turbulence values calculated by Equation (24).
where U 0 is the average velocity of incoming flow, I 0 is the atmospheric turbulence intensity, and µ t µ is the vortex-viscosity ratio.
The new source term is smaller than the dissipative term of the original equation for the wind turbine wakes region due to the introduced turbulence attenuation. The addition of the hold source term has little effect on the original SST k-ω model. Meanwhile, to avoid the numerical oscillation caused by the addition and switch of the hold source term in the free flow region and the non-free flow region, the hold source term is adopted in the whole computational domain.
The overestimation of the turbulent dissipation rate will result in the slow recovery rate of the wakes predicted by the RANS model [35]. The higher turbulent kinetic energy can promote the convection-diffusion between the wake region and the surrounding free fluid, and accelerate the recovery of wake velocity, which can be obtained by reducing the dissipation term. We reduce the turbulent dissipation by closure constant correction and correction factor addition, as discussed in the following subsection.

Closure Constant Correction
Based on the equilibrium flow theory, the surface friction velocity u* is introduced. The turbulent kinetic energy k is calculated by Equation (25).

Correction Factor Addition
The purpose of the dissipative correction is to modify the dissipation in the near wake region, but the previously modified equation is not universal in the whole. In this paper, a dissipative correction adaptive factor with the full computational domain is proposed. This method can highlight the effect of the correction term in the near wakes and reduce the correction in other regions. Dissipative item revision for the k-ω equation is expressed by Equation (27). where

Mesh Partition Method
The Semi-implicit method for pressure-linked equations (SIMPLE) solution is used under Fluent 6.0 solver considering two initial conditions. One is the neutral atmospheric condition, the other is the wind shear effect. The source term and turbulence are set by the user-defined function. The shear effect is represented by the logarithmic function.
where u * = τ ω /ρ is the surface friction velocity, τ ω is the shear stress, K is the Karman constant and the value is 0.41, z is the altitude from the ground, and z 0 is the surface roughness.
The 3D wind turbine model of the computational domain is completed by SolidWorks. The model is imported to the Ansys ICEM program to achieve mesh partition as shown in Figure 5. With the variable size regular hexahedron method, three different mesh average sizes (1/3D, 1/10D and 1/40D) are applied. The cartesian coordinate system is used, the center of the wind turbine is the coordinate origin, and the number of meshes is about 0.5 million. Compared with the general hexahedron sweep method, the developed method saves three-quarters of the mesh number.
The purpose of the dissipative correction is to modify the dissipation in the near wake region, but the previously modified equation is not universal in the whole. In this paper, a dissipative correction adaptive factor with the full computational domain is proposed. This method can highlight the effect of the correction term in the near wakes and reduce the correction in other regions. Dissipative item revision for the k-ω equation is expressed by Equation (27

Mesh Partition Method
The Semi-implicit method for pressure-linked equations (SIMPLE) solution is used under Fluent 6.0 solver considering two initial conditions. One is the neutral atmospheric condition, the other is the wind shear effect. The source term and turbulence are set by the user-defined function. The shear effect is represented by the logarithmic function.
where * /  u   is the surface friction velocity,   is the shear stress, K is the Karman constant and the value is 0.41, z is the altitude from the ground, and z0 is the surface roughness.
The 3D wind turbine model of the computational domain is completed by SolidWorks. The model is imported to the Ansys ICEM program to achieve mesh partition as shown in Figure 5. With the variable size regular hexahedron method, three different mesh average sizes (1/3D, 1/10D and 1/40D) are applied. The cartesian coordinate system is used, the center of the wind turbine is the coordinate origin, and the number of meshes is about 0.5 million. Compared with the general hexahedron sweep method, the developed method saves three-quarters of the mesh number.  Figure 6 shows the calculation domain dimension design. The computational domain is extended 25D in the streamwise direction (x-direction), 6D in the lateral direction (y-direction) and 6D in the vertical direction (z-direction). The wind turbine (mesh size is 1/40D) is located at 0D in the xdirection. The inner layer encryption region (mesh size is 1/10D) is close to the wind turbine. In the front view, the inlet of airflow direction is on the left end, and the exit is on the right.  Figure 6 shows the calculation domain dimension design. The computational domain is extended 25D in the streamwise direction (x-direction), 6D in the lateral direction (y-direction) and 6D in the vertical direction (z-direction). The wind turbine (mesh size is 1/40D) is located at 0D in the x-direction. The inner layer encryption region (mesh size is 1/10D) is close to the wind turbine. In the front view, the inlet of airflow direction is on the left end, and the exit is on the right.

Simulation Setup
To verify the effectiveness of the proposed method, experimental data are acquired from two wind turbines at Nibe, in Northern Denmark. The experiment data are tested by the actual measurement of the wind farm [38]. As shown in Figure 7, four meteorological masts (MMs) are

Simulation Setup
To verify the effectiveness of the proposed method, experimental data are acquired from two wind turbines at Nibe, in Northern Denmark. The experiment data are tested by the actual measurement of the wind farm [38]. As shown in Figure 7, four meteorological masts (MMs) are placed in a line at downstream distances, 2.5D, 4D, 6D and 7.5D, concerning the Nibe B wind turbine. The data consist of average values every 1 min for two years. The specific operating parameters of the wind turbines are shown in Table 2. The neutral atmospheric condition with wind shear effect is used as the boundary condition in this paper.

Simulation Setup
To verify the effectiveness of the proposed method, experimental data are acquired from two wind turbines at Nibe, in Northern Denmark. The experiment data are tested by the actual measurement of the wind farm [38]. As shown in Figure 7, four meteorological masts (MMs) are placed in a line at downstream distances, 2.5D, 4D, 6D and 7.5D, concerning the Nibe B wind turbine. The data consist of average values every 1 min for two years. The specific operating parameters of the wind turbines are shown in Table 2. The neutral atmospheric condition with wind shear effect is used as the boundary condition in this paper.  Six modified items in two correction modules are proposed to improve the accuracy of the wake calculation, which are tip loss correction (Equation (7)), hub loss correction (Equation (9)), attack angle correction (Equation (10)), turbulence attenuation correction (Equation (23)), closure constant correction (Equation (26)) and dissipative item correction (Equation (27)). Thus, six models are designed and compared as shown in Table 3.  Six modified items in two correction modules are proposed to improve the accuracy of the wake calculation, which are tip loss correction (Equation (7)), hub loss correction (Equation (9)), attack angle correction (Equation (10)), turbulence attenuation correction (Equation (23)), closure constant correction (Equation (26)) and dissipative item correction (Equation (27)). Thus, six models are designed and compared as shown in Table 3.  Figure 8 shows the influence of each correction item on the change in wake velocity. The comparison between model 1 and exe1 (measured at Nibe) data shows that all the six correction items can ensure the agreement between the calculated results and the experimental values. The overall velocity changes to the lowest value of 0.48 at the position of 2.5D behind the wind turbine, increases rapidly in the near wake region, and reaches 0.64 and 0.75 at 4D and 6D, respectively. The overall velocity slowly increases in the far wake region and reaches 0.85 and 0.97 at 7.5D and 20D, respectively.  Figure 8 shows the influence of each correction item on the change in wake velocity. The comparison between model 1 and exe1 (measured at Nibe) data shows that all the six correction items can ensure the agreement between the calculated results and the experimental values. The overall velocity changes to the lowest value of 0.48 at the position of 2.5D behind the wind turbine, increases rapidly in the near wake region, and reaches 0.64 and 0.75 at 4D and 6D, respectively. The overall velocity slowly increases in the far wake region and reaches 0.85 and 0.97 at 7.5D and 20D, respectively. The comparison between models 4-6 and models 1-3 shows the effects of the modified turbulence equation on the wake velocity recovery, especially in the far wake region. At 20D, the velocity recovery of models 4-6 is slower than that of model 1 in the far wake region with 0.79, 0.81, 0.87, respectively, while models 2 and 3 can still show a faster recovery velocity in the far wake region with 0.88, 0.89, respectively, which indicates that the three turbulence correction terms could affect the far wakes. The recovery of far wakes is significantly higher in model 1 than models 5 and 6. This shows that increasing the sealing parameter and reducing the dissipative term affect the increase in velocity recovery. The comparison between models 4-6 and models 1-3 shows the effects of the modified turbulence equation on the wake velocity recovery, especially in the far wake region. At 20D, the velocity recovery of models 4-6 is slower than that of model 1 in the far wake region with 0.79, 0.81, 0.87, respectively, while models 2 and 3 can still show a faster recovery velocity in the far wake region with 0.88, 0.89, respectively, which indicates that the three turbulence correction terms could affect the far wakes. The recovery of far wakes is significantly higher in model 1 than models 5 and 6. This shows that increasing the sealing parameter and reducing the dissipative term affect the increase in velocity recovery. Figure 9a shows that the correction of the wind turbine model wake field is obvious in the near wake region. The comparison between model 2 and model 3 shows that the modification of tip loss and hub loss can effectively narrow the profile velocity change curve and reduce the expansion range of the wake flow field. At the same time, the velocity valley value is increased from 0.38 to 0.42 due to the decreasing of the volume source term by the blade tip and hub modification. Meanwhile, the valley difference between model 1 and model 4 is only 2%. The whole curvature is nearly closed to model 1, which indicates that the correction of the turbulence equation in the near wakes has less effect. Figure 9b indicates that the effect of turbulence correction increases gradually with the increase in horizontal distance, while the effect of wind turbine correction gradually decreases. The comparison between model 1 and model 4 indicates that the wake velocity recovery without turbulence correction is not satisfactory, as the velocity changes from higher than 2% at 2.5D to lower than 5% at 4D. Figure 9c shows that the correction of the turbulence equation in the far wakes of the wind turbine is obvious. Compared with models 4 and 5, the correction effect of the decreasing dissipation term in Energies 2020, 13, 4430 13 of 19 model 6 is obvious in the range of −1.5D to 1.5D in the y-direction, which causes the wind velocity curve to reach to the maximum. On the whole, comparing the 7.5D section with the 2.5D section, the velocity from the 0.5D to 1D in the y-direction gradually decreases in the axial direction. Figure 9a shows that the correction of the wind turbine model wake field is obvious in the near wake region. The comparison between model 2 and model 3 shows that the modification of tip loss and hub loss can effectively narrow the profile velocity change curve and reduce the expansion range of the wake flow field. At the same time, the velocity valley value is increased from 0.38 to 0.42 due to the decreasing of the volume source term by the blade tip and hub modification. Meanwhile, the valley difference between model 1 and model 4 is only 2%. The whole curvature is nearly closed to model 1, which indicates that the correction of the turbulence equation in the near wakes has less effect.  Figure 9b indicates that the effect of turbulence correction increases gradually with the increase in horizontal distance, while the effect of wind turbine correction gradually decreases. The comparison between model 1 and model 4 indicates that the wake velocity recovery without turbulence correction is not satisfactory, as the velocity changes from higher than 2% at 2.5D to lower than 5% at 4D. Figure 9c shows that the correction of the turbulence equation in the far wakes of the wind turbine is obvious. Compared with models 4 and 5, the correction effect of the decreasing dissipation term in model 6 is obvious in the range of −1.5D to 1.5D in the y-direction, which causes the wind  Figure 10 demonstrates that the valley position moves forward slightly from 2.5D to the horizontal distance with the increasing height. There is no significant change in the near and far wake regions. The wake drop slows down after 0.3D and increases from 3% to 10-20% per 0.1D in the direction of height. Figure 11 indicates that the expansion effect is obvious in the 10D range behind the wake, and shows a nonlinear boundary curve. The expansion range decreases gradually, and the nonlinearity of the expansion boundary curve decreases with the increasing elevation. When compared with Figure 11a,b, the velocity valley decreases from 0.48 at 0D to 0.58 at 0.3D, and the velocity gradient decreases gradually. The range of wake expansion at 0D is a reference for the wake velocity model, which is used in the microscopic site selection of a flat surface. Figure 10 demonstrates that the valley position moves forward slightly from 2.5D to the horizontal distance with the increasing height. There is no significant change in the near and far wake regions. The wake drop slows down after 0.3D and increases from 3% to 10-20% per 0.1D in the direction of height.  Figure 11 indicates that the expansion effect is obvious in the 10D range behind the wake, and shows a nonlinear boundary curve. The expansion range decreases gradually, and the nonlinearity of the expansion boundary curve decreases with the increasing elevation. When compared with Figure 11a,b, the velocity valley decreases from 0.48 at 0D to 0.58 at 0.3D, and the velocity gradient decreases gradually. The range of wake expansion at 0D is a reference for the wake velocity model, which is used in the microscopic site selection of a flat surface.   Figure 11 indicates that the expansion effect is obvious in the 10D range behind the wake, and shows a nonlinear boundary curve. The expansion range decreases gradually, and the nonlinearity of the expansion boundary curve decreases with the increasing elevation. When compared with Figure 11a,b, the velocity valley decreases from 0.48 at 0D to 0.58 at 0.3D, and the velocity gradient decreases gradually. The range of wake expansion at 0D is a reference for the wake velocity model, which is used in the microscopic site selection of a flat surface. (b) 0.3D Figure 11. Wake velocity clouds on the neutral layer of different heights.

Turbulence Field Result under the Uniform Inflow Condition
It can be seen from Figure 12 that the initial turbulence intensity of models 1, 5 and 6 is maintained in the free flow region. The result of model 4 shows that the turbulence attenuation is not corrected. In the near wake region of the wind turbine, the turbulence intensity simulated by models 1, 5 and 6 is higher than that of model 4. Notably, the result of model 5 is about 30% higher than that of model 4. With the increasing position, the wake effect starts to weaken, and the turbulence intensity

Turbulence Field Result under the Uniform Inflow Condition
It can be seen from Figure 12 that the initial turbulence intensity of models 1, 5 and 6 is maintained in the free flow region. The result of model 4 shows that the turbulence attenuation is not corrected. In the near wake region of the wind turbine, the turbulence intensity simulated by models 1, 5 and 6 is higher than that of model 4. Notably, the result of model 5 is about 30% higher than that of model 4. With the increasing position, the wake effect starts to weaken, and the turbulence intensity decreases correspondingly. The performance of model 1 is the best, as the result of the modified model is also extremely close to the experimental value with 16% at the peak. As the dissipative term of model 1 is about twice that of the original one, the specific dissipation rate is lower, and the turbulence intensity is larger.
(b) 0.3D Figure 11. Wake velocity clouds on the neutral layer of different heights.

Turbulence Field Result under the Uniform Inflow Condition
It can be seen from Figure 12 that the initial turbulence intensity of models 1, 5 and 6 is maintained in the free flow region. The result of model 4 shows that the turbulence attenuation is not corrected. In the near wake region of the wind turbine, the turbulence intensity simulated by models 1, 5 and 6 is higher than that of model 4. Notably, the result of model 5 is about 30% higher than that of model 4. With the increasing position, the wake effect starts to weaken, and the turbulence intensity decreases correspondingly. The performance of model 1 is the best, as the result of the modified model is also extremely close to the experimental value with 16% at the peak. As the dissipative term of model 1 is about twice that of the original one, the specific dissipation rate is lower, and the turbulence intensity is larger. As shown in Figure 13, the result of the 2.5D position shows that the standard SST k-ω model underestimates the magnitude of turbulence intensity qualitatively, but the model can well predict the "double peak" effect of turbulence intensity. The "double peak" effect of turbulent intensity weakens gradually with the development and fragmentation of the tip vortex convection-diffusion As shown in Figure 13, the result of the 2.5D position shows that the standard SST k-ω model underestimates the magnitude of turbulence intensity qualitatively, but the model can well predict the "double peak" effect of turbulence intensity. The "double peak" effect of turbulent intensity weakens gradually with the development and fragmentation of the tip vortex convection-diffusion of wakes, and the shear mixing layer continues to expand, which can be seen from the distribution of turbulence intensity at 4D. Through the comparison of models 1 and 3-6, it can be seen that the modified wind turbine model has no obvious effect on the turbulence field. At the same time, the effect of each turbulence correction term is obvious. The proposed model is more effective than model 5.
The turbulence intensity above the Z = 0D plane is shown in Figure 14. The maximum turbulence intensity is 21.8%. The whole variation of the "double peak" effect can be seen. In the y-direction, the turbulent peak occurs mainly at ±0.5D (the blade tip), and the diffusion range is from −1D to 1D.
In the x-direction, a rapid decay state that approaches the inlet turbulence value after 2.5 D is observed. of wakes, and the shear mixing layer continues to expand, which can be seen from the distribution of turbulence intensity at 4D. Through the comparison of models 1 and 3-6, it can be seen that the modified wind turbine model has no obvious effect on the turbulence field. At the same time, the effect of each turbulence correction term is obvious. The proposed model is more effective than model 5. The turbulence intensity above the Z = 0D plane is shown in Figure 14. The maximum turbulence intensity is 21.8%. The whole variation of the "double peak" effect can be seen. In the y-direction, the turbulent peak occurs mainly at ±0.5D (the blade tip), and the diffusion range is from −1D to 1D. In the x-direction, a rapid decay state that approaches the inlet turbulence value after 2.5 D is observed.  Table 4 shows the computational cost of different schemes.  The turbulence intensity above the Z = 0D plane is shown in Figure 14. The maximum turbulence intensity is 21.8%. The whole variation of the "double peak" effect can be seen. In the y-direction, the turbulent peak occurs mainly at ±0.5D (the blade tip), and the diffusion range is from −1D to 1D. In the x-direction, a rapid decay state that approaches the inlet turbulence value after 2.5 D is observed.  Table 4 shows the computational cost of different schemes.  Table 4 shows the computational cost of different schemes. Compared with Denmark's 2D mesh scheme (Scheme A) [25], the mesh amount of our scheme is slightly larger than scheme A. However, it should be noted that scheme A sacrifices the third dimension of mesh modeling. Compared with Xu's half-axis hypothesis 3D scheme (Scheme B) [34], the mesh amount of our scheme is substantially less than scheme B under the same average size, although scheme B uses axisymmetric assumptions to build only a half 3D model.

Conclusions
This paper proposes a modified Reynolds-averaged Navier-Stokes (MRANS)-based wind turbine wake model to simulate wake effects. Based on the correction module, the proposed BEM-fuzzy aerodynamic model can amend the inconsistent condition between wake simulation and experiment tests, which affects the simulation of the near wake. For the turbulence model, the turbulence attenuation is effectively avoided by adding the hold source term to ensure the global correctness of boundary conditions. The accuracy of the turbulence intensity distribution is improved by correcting the closure constant and the dissipation term. The turbulence model affects the whole wakes, especially the far wakes.
The simulation results of the velocity field and turbulent field with the proposed approach are consistent with the data of real wind turbines, which verifies the effectiveness of the proposed approach. Furthermore, the computation efficiency is significantly improved by the developed mesh partition method.