Study on Actuator Line Modeling of Two NREL 5-MW Wind Turbine Wakes

: The wind turbine wakes impact the efﬁciency and lifespan of the wind farm. Therefore, to improve the wind plant performance, research on wind plant control is essential. The actuator line model (ALM) is proposed to simulate the wind turbine efﬁciently. This research investigates the National Renewable Energy Laboratory 5 Million Watts (NREL 5-MW) wind turbine wakes with Open Field Operation and Manipulation (OpenFOAM) using ALM. Firstly, a single NREL 5-MW turbine is simulated. The comparison of the power and thrust with Fatigue, Aerodynamics, Structures, and Turbulence (FAST) shows a good agreement below the rated wind speed. The information relating to wind turbine wakes is given in detail. The top working status is proved at the wind speed of 8 m/s and the downstream distance of more than 5 rotor diameters (5 D ). Secondly, another case with two NREL 5-MW wind turbines aligned is also carried out, in which 7 D is validated as the optimum distance between the two turbines. The result also shows that the upstream wind turbine has an obvious inﬂuence on the downstream one.


Introduction
The flow passing through a wind turbine always makes the performance of its downstream wind turbine less efficient than that which was designed. This influence can decrease 10-20% of the production of the total wind farm [1] and increase 5-15% of the loading on the wind turbine blades [2,3]. So it is significant to design the wind turbine by analyzing the wake-induced power and fatigue. Generally speaking, these wakes are very complicated and difficult to reproduce.
There are usually two methods employed to calculate the performance and loading of the wind turbines. One is the Blade-Element/Momentum (BEM) theory [4,5]. It is simple and computationally efficient, but needs to introduce some three-dimensional corrections, tip loss, dynamic stall, etc., to approximate the true situation. The other is the computational fluid dynamics (CFD) method. It can solve the physical problem accurately but always runs on a computer with consuming time. Despite this, many related researches on the simulation of the wind turbine have been documented by using Computational Fluid Dynamics (CFD) [6][7][8]. Besides, there also exist some other models that can always be used, including the actuator disc model [9], the actuator line model (ALM) [10], and the actuator surface model [11]. In these actuator models, the body forces are distributed in the flow field. Initially, the actuator disc model was used in the wind turbine far wake simulation. However, it could not reproduce the tip vortex system accurately [12]. Then, ALM was introduced by Shen and Sørensen [10] to overcome this difficulty, where the body forces were distributed radially along the blades and added to the three dimensional Navier-Stokes equations as the source item. Moreover, ALM is also an effective technique for predicting the loading on the blades with only simple structured grids, which has been proved by Troldborg [12]. Recently, Meng et al. [13] firstly proposed the concept of the elastic actuator line (EAL), which combines ALM and a finite difference structural model. In the EAL model, the elastic blades are incorporated into ALM and its simulation results are closer to the reality. Nevertheless, there still exists some limitations in EAL and this needs further study for the improvement of the theory of the actuator line.
For quite a while, most models using ALM were realized by the in-house code, and there is very little commercial ALM software before 2013. Then, the National Renewable Energy Laboratory (NREL) in America developed a CFD tool called Simulator for Off/Onshore Wind Farm Applications (SOWFA) [14], which coupled the OpenFOAM [15] and FAST [16]. It contains the ALM and the atmospheric boundary layer (ABL) model. Fleming et al. [17,18] used it to simulate the wake of the two-turbine case and the result was accurate enough for a practical purpose. Recently, Bachant et al. [19] have developed a library named turbinesFoam, which uses ALM to simulate the hydrokinetics of wind and marine turbines in OpenFOAM. Its interpolation, Gaussian projection, and vector rotation functions are all adapted from NREL's SOWFA.
In this paper, the contents of each section are organized as follows: Firstly, the concept of ALM is reviewed and the tip loss corrections [20] for the computations are given in brief. Then, an uncertainty analysis is carried out through a typical airfoil. After this, the single NREL 5-MW wind turbine case is used as the computational validation. Finally, the influence of the downstream wind turbine rotor speed and the separation distance between the two wind turbines on the whole system performance is examined. Moreover, the velocity fields in the wake flow are also analyzed. Considering the size of the domain, turbinesFoam is used to implement the computational task. The main reference data and the airfoil data can be found from the report of NREL [21].

Basic Theory of Actuator Line
The Navier-Stokes (N-S) equations and the actuator line concept are combined as follows: where V is the wind speed; t is the time; ρ is the air density, and ρ = 1.225 kg/m 3 in this paper; p is the pressure; υ is the kinematic viscosity, and it is set as 1.5 × 10 −5 m 2 /s; and f ε is the loading on the blade as the body force that is added into the momentum equations. The body forces are located along the actuator lines, which represent the blades, as shown in Figure 1, and they are calculated by the blade-element theory with two-dimensional airfoil data. The lift and drag acting on the blade element are, respectively: where L and D represent the lift and drag, respectively; C l (α) and C d (α) are the lift and drag coefficients, correspondingly; α is the local angle of attack; U rel is the local velocity relative to the blade; c is the chord length; and dr is the thickness of the blade element. Figure 2 shows that local velocity relative to the blade can be calculated by: where U z and U θ are the axial and tangential velocity, respectively; Ω is the rotating angular velocity; and r is the rotating radius. Additionally, the angle of attack is: where φ is the angle of the inflow and γ is the local pitch angle. Therefore, the force per spanwise unit length is: where → e L and → e D are the unit direction vectors of L and D, respectively. As the source term in Equation (1) is given in the form of the curl of the load, it acts as a singular vorticity source along the rotor blades. To avoid singularity, a constant ε is added to adjust the strength of the regularization kernel function, and one form can be shown as: where d i is the distance between the measured point (x i , y i , z i ) and initial force points (x, y, z) on the blade. η ε is the kernel function between the measured and initial force points, which is quite similar to the kernel function used in Smoothed Particle Hydrodynamics (SPH) [22]. Finally, f ε on the nearby mesh can be calculated by: where N is the number of neighbouring blade sections. It means that the discreted force on each blade section can be interpolated to the neigbouring mesh nodes smoothly. More details on the solution of Equation (9) can be referred to in [10]. To add the tower and hub effect into ALM, these will be achieved by following a similar technology. where z U and U θ are the axial and tangential velocity, respectively; Ω is the rotating angular velocity; and r is the rotating radius. Additionally, the angle of attack is: where φ is the angle of the inflow and γ is the local pitch angle. Therefore, the force per spanwise unit length is: where L e  and D e  are the unit direction vectors of L and D, respectively.
As the source term in Equation (1) is given in the form of the curl of the load, it acts as a singular vorticity source along the rotor blades. To avoid singularity, a constant ε is added to adjust the strength of the regularization kernel function, and one form can be shown as: where di is the distance between the measured point (xi, yi, zi) and initial force points (x, y, z) on the blade. ε η is the kernel function between the measured and initial force points, which is quite similar to the kernel function used in Smoothed Particle Hydrodynamics (SPH) [22]. Finally, f ε on the nearby mesh can be calculated by: where N is the number of neighbouring blade sections. It means that the discreted force on each blade section can be interpolated to the neigbouring mesh nodes smoothly. More details on the solution of Equation (9) can be referred to in [10]. To add the tower and hub effect into ALM, these will be achieved by following a similar technology.   where z U and U θ are the axial and tangential velocity, respectively; Ω is the rotating angular velocity; and r is the rotating radius. Additionally, the angle of attack is: where φ is the angle of the inflow and γ is the local pitch angle. Therefore, the force per spanwise unit length is: where L e  and D e  are the unit direction vectors of L and D, respectively.
As the source term in Equation (1) is given in the form of the curl of the load, it acts as a singular vorticity source along the rotor blades. To avoid singularity, a constant ε is added to adjust the strength of the regularization kernel function, and one form can be shown as: where di is the distance between the measured point (xi, yi, zi) and initial force points (x, y, z) on the blade. ε η is the kernel function between the measured and initial force points, which is quite similar to the kernel function used in Smoothed Particle Hydrodynamics (SPH) [22]. Finally, f ε on the nearby mesh can be calculated by: where N is the number of neighbouring blade sections. It means that the discreted force on each blade section can be interpolated to the neigbouring mesh nodes smoothly. More details on the solution of Equation (9) can be referred to in [10]. To add the tower and hub effect into ALM, these will be achieved by following a similar technology.

Tip Loss Correction Model
In order to simulate a real wind turbine, the tip loss effects must be introduced. The tip loss correction model is used here to take into account the difference between the ALM and the actual blades, where the velocity should always be zero at the tip. This always makes the actual loading on the blades different to the loading we obtain. This concept was introduced by Prandtl, and refined by many scholars. In the present research, the Prandtl-Glauert model and Shen model [20] are used. The Prandtl-Glauert tip loss function is introduced firstly by, where B is the number of the blades; R is the radius of the blades; r is the distance between the blade-element location and the root of the blade; and φ is the angle between the local relative velocity and the rotor plane. Then, Shen improved the function into: where g is a coefficient that depends on the tip speed ratio, numbers of blades, pitch setting, and chord distribution, etc. Furthermore, he also assumed g as the following form: where c 1 and c 2 are the two coefficients determined from the experimental data. They are set as c 1 = 0.125 and c 2 = 21 in this paper.

Refrence Data of NREL 5-MW
The NREL 5-MW wind turbine is a conventional three-blade upwind variable-speed variable-bladepitch-to-feather-controlled turbine. Its rotor diameter is about 126 m and the hub height is nearly 90 m. The rated wind speed is 11.4 m/s. Under the rated wind speed, the rotor speed is 12.1 rpm and the power is 5 MW. The structure of the blades used in ALM is shown in Table 1. Additionally, the radius of the tower top is 3 m and it is linearly distributed along the tower.

Numerical Model
Flow occurs around a NACA 0012 airfoil at Reynolds number based on the airfoil chord c, incoming velocity V ∞ , and fluid kinematic viscosity ν of Re = 6 × 10 6 at two angles of attack β, β = 4 • and 10 • . For each attack angle, there are nine geometrically similar grids with available refinement based on OpenFOAM. The domain for the calculation of the flow over around the NACA 0012 airfoils is shown in Figure 3, where c is the chord of the airfoil. The inlet is a semicircle with r = 12c and the outlet is placed 13c downstream of the trailing edge of the airfoil. The external boundary is located approximately 12c away from the airfoil. The Reynolds number based on the velocity of the incoming flow V ∞ , chord of the foil c.

Numerical Model
Flow occurs around a NACA 0012 airfoil at Reynolds number based on the airfoil chord c, incoming velocity V ∞ , and fluid kinematic viscosity ν of 6 6 10 Re = × at two angles of attack β, β = 4° and 10°. For each attack angle, there are nine geometrically similar grids with available refinement based on OpenFOAM. The domain for the calculation of the flow over around the NACA 0012 airfoils is shown in Figure 3, where c is the chord of the airfoil. The inlet is a semicircle with r = 12c and the outlet is placed 13c downstream of the trailing edge of the airfoil. The external boundary is located approximately 12c away from the airfoil. The Reynolds number based on the velocity of the incoming flow V ∞ , chord of the foil c. Distribution of grid nodes is defined by the one-sided stretching functions proposed by Vinokur [23]. For the sake of completeness, the largest, average, and smallest dimensionless distance of the near-wall cell centre to the airfoil surface are included in the file that contains the lift ( L C ) and drag ( d C ) coefficients of the airfoil. u τ is the friction velocity defined by u ω τ τ ρ = , where ω τ is the shear-stress at the wall and ρ is the fluid density. The grid refinement ratio ri may be obtained from the total number of cells Ncells: The number of cells is independent of the angle of attack β and the values of Ncells and ri are given in Table 2. x/c y/c Distribution of grid nodes is defined by the one-sided stretching functions proposed by Vinokur [23]. For the sake of completeness, the largest, average, and smallest dimensionless distance of the near-wall cell centre to the airfoil surface are included in the file that contains the lift (C L ) and drag (C d ) coefficients of the airfoil. u τ is the friction velocity defined by u τ = τ ω ρ , where τ ω is the shear-stress at the wall and ρ is the fluid density. The grid refinement ratio r i may be obtained from the total number of cells N cells : The number of cells is independent of the angle of attack β and the values of N cells and r i are given in Table 2.

Uncertainty Analysis
There are various existing verification methods for solution verification, and it is a common practice to retain only the first term of the Richardson extrapolation and assume that the solutions are in the asymptotic range, which leads to a grid-triplet study. The Grid Convergence Index (GCI) derived by Roache [24] can be used to estimate the uncertainties due to grid-spacing and time-step errors and is widely used and recommended by American Society of Mechanical Engineers (ASME) [25] and American Institute of Aeronautics and Astronautics (AIAA) [26]. The other two common verification methods in ship hydrodynamics are Factor of Safety methods (FS) developed by Xing and Stern [27] and Least Squared Root methods (LSR) introduced by Eça [28]. The Least Squared Root methods (LSR) are applied in this paper.
Nowadays, the most common verification methods in CFD are based on Richardson extrapolation. Then, discrete solutions ϕ are assumed to: In which ϕ 0 is the exact value, p is the order of convergence, and α is a constant. The LSR is based on power series expansions that neglect high-order terms and assume that ϕ has at least second-order finite derivatives. It is also assumed that the lowest-order schemes used in the discretization are second or first-order accurate. The present subsection will show the procedure for numerical uncertainty estimation with the latest LSR method in the following: (1) Determination of ε φ Solve the non-linear equations, Equation (15), in the least-squares sense with and without weights to obtain δ RE (δ RE = φ i − φ 0 ), p, and the standard deviations of the two fits σ.
where w is the weighted number, w i = 1, and nw i = 1 in the non-weighted approach. For the weighted approach, w i = 1/h i ng ∑ i=1 1/h i and n g is the grid number, which is supposed to be at least 4. The standard deviation is given by: If any of the fits exhibit 0.5 ≤ p ≤ 2, ε φ = δ RE . If both fits exhibit 0.5 ≤ p ≤ 2, the value of δ RE selected corresponds to the fit with the smallest standard deviation.
If the observed order of grid convergence is p > 2, solve the non-linear equations, Equations (17)- (19), in the least-squares sense with and without weights and determine the standard deviations of the six fits σ. ε φ is obtained from the fit that exhibits the smallest standard deviation.
(2) Determine the following data range parameter to assess the quality of the fit used to obtain the error estimate ε φ (3) Determine the safety factor If 0.5 ≤ p < 2.1 and σ < ∆ φ , F s = 1.25; otherwise, F s = 3.
(4) Obtain the uncertainty For σ < ∆ φ : For σ ≥ ∆ φ : The lift (C L ) and drag (C d ) coefficients with different attack angles are shown in Table 3. By using the LSR uncertainty estimator which is presented in the previous section, the convergence of lift and drag coefficients with the grid refinement are shown in Figures 4 and 5, respectively. The lift coefficient is nearly zero when the angle of attack is zero, so these lift coefficients are not analyzed. The error bars for the mesh ratio equal to 1.0, 2.0, and 4.0 are also included in those figures. Then, Table 4 displays the uncertainty of the lift and drag coefficient with the grid refinement for different angle of attacks. It is clear that the error and uncertainty are reduced with the mesh refinement. Finally, we can bound the exact value (φ exact ) and the interval is: Table 5 shows the interval of the exact lift and drag coefficient with the grid refinement for different angle of attacks. Figure 6 shows the uncertainty of the drag and lift coefficient with the grid refinement for different attack angles. Considering the experiment results of NACA 0012 by Ladson [29] and McCroskey [30], the results are well bounded by the OpenFOAM uncertainty interval.
For φ σ ≥ Δ : The lift ( L C ) and drag ( d C ) coefficients with different attack angles are shown in Table 3. By using the LSR uncertainty estimator which is presented in the previous section, the convergence of lift and drag coefficients with the grid refinement are shown in Figures 4 and 5, respectively. The lift coefficient is nearly zero when the angle of attack is zero, so these lift coefficients are not analyzed. The error bars for the mesh ratio equal to 1.0, 2.0, and 4.0 are also included in those figures. Then, Table 4 displays the uncertainty of the lift and drag coefficient with the grid refinement for different angle of attacks. It is clear that the error and uncertainty are reduced with the mesh refinement. Finally, we can bound the exact value ( exact φ ) and the interval is: Table 5 shows the interval of the exact lift and drag coefficient with the grid refinement for different angle of attacks. Figure 6 shows the uncertainty of the drag and lift coefficient with the grid refinement for different attack angles. Considering the experiment results of NACA 0012 by Ladson [29] and McCroskey [30], the results are well bounded by the OpenFOAM uncertainty interval.

A Single NREL 5-MW Turbine Case
In this section, a single NREL 5-MW wind turbine case is built to identify the accuracy of the model. Additionally, there are three cases with different representative wind speeds to find the reason why the wind speed is 8 m/s in some studies. These wind speeds are set as 11.4 m/s, 8 m/s, and 5 m/s, respectively. The wind turbine obtains the rated status at the wind speed of 11.4 m/s, and the tip speed ratio is 7. When the wind speed is 8 m/s, the power coefficient is up to the maximum and the corresponding tip speed ratio is 7.55. Compared with the two cases above, the last case is set, in which the wind speed is used as 5 m/s.

Numerical Model
The wind turbine is placed in a domain of 1890 m × 1260 m × 468 m. Details are shown in Figure 7, which shows that the mesh is refined by three steps in a rectangular region. The size of the cell outside the refinement zone is approximately 24 m; therefore, the inner zone mesh size is nearly 3 m. In addition, the regions behind the blade tip and root are important for ascertaining an accurate wake structure, so these parts are refined by the smaller mesh, which is indicated in Figure 8. In order to reduce the artificial reflections and oscillations, the refined ratio is set to 2 and three levels are adopted to avoid a large refinement ratio. Meshes are generated by blockMesh and snappyHexMesh in OpenFOAM. The wind turbine and the tower are set by topoSet.
The wind flow is steady. The inflow and outflow boundary conditions of U are set by the fixedValue and inletOutlet, respectively, and zeroGradient and fixedValue are used correspondingly for the pressure p. All others are defined as the slip boundary. The Reynolds-averaged Navier-Stokes (RANS) equations using the shear stress transport (SST) k − ω turbulence model are considered in this model. The pressure-velocity coupling solver application uses the combined the Pressure-Implicit split-operator (PISO) and the semi-implicit Method for Pressure-Linked Equations (SIMPLE) algorithms, which be can be abbrivated as PIMPLE. The number of times that the algorithm solves the pressure equation and momentum corrector in each step is 2. The solver for p is GAMG (Generalized geometric-Algebraic Multi-Grid) and for others is smoothSolver, which is a solver that uses a smoother. The numerical schemes are set as Euler for the time marching schemes, which is first-order implicit, with the Gauss liner for the gradient term, and Gauss linerUpwind for the divergence term. Besides, the time-step size is set corresponding to a blade rotation of 1 degree and the data is saved every five time steps. To obtain a relatively stable state, the simulation is calculated for 20 revolutions of the wind turbine. The Glauert model is chosen as the tip loss model.

A Single NREL 5-MW Turbine Case
In this section, a single NREL 5-MW wind turbine case is built to identify the accuracy of the model. Additionally, there are three cases with different representative wind speeds to find the reason why the wind speed is 8 m/s in some studies. These wind speeds are set as 11.4 m/s, 8 m/s, and 5 m/s, respectively. The wind turbine obtains the rated status at the wind speed of 11.4 m/s, and the tip speed ratio is 7. When the wind speed is 8 m/s, the power coefficient is up to the maximum and the corresponding tip speed ratio is 7.55. Compared with the two cases above, the last case is set, in which the wind speed is used as 5 m/s.

Numerical Model
The wind turbine is placed in a domain of 1890 m × 1260 m × 468 m. Details are shown in Figure  7, which shows that the mesh is refined by three steps in a rectangular region. The size of the cell outside the refinement zone is approximately 24 m; therefore, the inner zone mesh size is nearly 3 m. In addition, the regions behind the blade tip and root are important for ascertaining an accurate wake structure, so these parts are refined by the smaller mesh, which is indicated in Figure 8. In order to reduce the artificial reflections and oscillations, the refined ratio is set to 2 and three levels are adopted to avoid a large refinement ratio. Meshes are generated by blockMesh and snappyHexMesh in OpenFOAM. The wind turbine and the tower are set by topoSet.
The wind flow is steady. The inflow and outflow boundary conditions of U are set by the fixedValue and inletOutlet, respectively, and zeroGradient and fixedValue are used correspondingly for the pressure p. All others are defined as the slip boundary. The Reynolds-averaged Navier-Stokes (RANS) equations using the shear stress transport (SST) k − ω turbulence model are considered in this model. The pressure-velocity coupling solver application uses the combined the Pressure-Implicit split-operator (PISO) and the semi-implicit Method for Pressure-Linked Equations (SIMPLE) algorithms, which be can be abbrivated as PIMPLE. The number of times that the algorithm solves the pressure equation and momentum corrector in each step is 2. The solver for p is GAMG (Generalized geometric-Algebraic Multi-Grid) and for others is smoothSolver, which is a solver that uses a smoother. The numerical schemes are set as Euler for the time marching schemes, which is first-order implicit, with the Gauss liner for the gradient term, and Gauss linerUpwind for the divergence term. Besides, the time-step size is set corresponding to a blade rotation of 1 degree and the data is saved every five time steps. To obtain a relatively stable state, the simulation is calculated for 20 revolutions of the wind turbine. The Glauert model is chosen as the tip loss model.

Validation of Numerical Model
In order to quantify the influence of the mesh quality, a series of computations with variable grid numbers are carried out for U0 = 11.4 m/s, which are 375,432, 744,792, 2,903,144, and 5,958,336, respectively. As there are no experimental data for a full scale NREL 5-WM wind turbine, it adopts the results of FAST as the benchmark data, which can be referred to in [21]. The simulation time is set 30 s. The power and thrust coefficients are introduced and defined as below: where P is the mechanical power; T is the thrust on the blades; and R is the radius of the rotor blade and R = 63 m, ρ = 1.225 kg/m 3 . According to reference [21], Cp = 0.4418, Ct = 0.8099 when U0 = 11.4 m/s. The errors of Cp and Ct of different cell numbers are shown in Table 6. The errors of the power and thrust are compared in Figure 9, which shows that when the numbers of the grid increase, the differences in the power and thrust begin to stabilize. Considering the computational efficiency, the total number of 2,903,144 cells should be high enough to obtain accurate results, considering the errors of Cp and Ct which are 4.858% and 5.981%, respectively.

Validation of Numerical Model
In order to quantify the influence of the mesh quality, a series of computations with variable grid numbers are carried out for U0 = 11.4 m/s, which are 375,432, 744,792, 2,903,144, and 5,958,336, respectively. As there are no experimental data for a full scale NREL 5-WM wind turbine, it adopts the results of FAST as the benchmark data, which can be referred to in [21]. The simulation time is set 30 s. The power and thrust coefficients are introduced and defined as below: where P is the mechanical power; T is the thrust on the blades; and R is the radius of the rotor blade and R = 63 m, ρ = 1.225 kg/m 3 . According to reference [21], Cp = 0.4418, Ct = 0.8099 when U0 = 11.4 m/s. The errors of Cp and Ct of different cell numbers are shown in Table 6. The errors of the power and thrust are compared in Figure 9, which shows that when the numbers of the grid increase, the differences in the power and thrust begin to stabilize. Considering the computational efficiency, the total number of 2,903,144 cells should be high enough to obtain accurate results, considering the errors of Cp and Ct which are 4.858% and 5.981%, respectively.

Validation of Numerical Model
In order to quantify the influence of the mesh quality, a series of computations with variable grid numbers are carried out for U 0 = 11.4 m/s, which are 375,432, 744,792, 2,903,144, and 5,958,336, respectively. As there are no experimental data for a full scale NREL 5-WM wind turbine, it adopts the results of FAST as the benchmark data, which can be referred to in [21]. The simulation time is set 30 s. The power and thrust coefficients are introduced and defined as below: where P is the mechanical power; T is the thrust on the blades; and R is the radius of the rotor blade and R = 63 m, ρ = 1.225 kg/m 3 . According to reference [21], C p = 0.4418, C t = 0.8099 when U 0 = 11.4 m/s. The errors of C p and C t of different cell numbers are shown in Table 6. The errors of the power and thrust are compared in Figure 9, which shows that when the numbers of the grid increase, the differences in the power and thrust begin to stabilize. Considering the computational efficiency, the total number of 2,903,144 cells should be high enough to obtain accurate results, considering the errors of C p and C t which are 4.858% and 5.981%, respectively.  To examine the accuracy of different wind speeds, the thrust and power are computed and compared with the FAST code [21], which are shown in Figure 10. The results agree with FAST excellently for the wind speeds up to 11.4 m/s. The quantitative values of interest and errors are shown in Table 7.    To examine the accuracy of different wind speeds, the thrust and power are computed and compared with the FAST code [21], which are shown in Figure 10. The results agree with FAST excellently for the wind speeds up to 11.4 m/s. The quantitative values of interest and errors are shown in Table 7.  To examine the accuracy of different wind speeds, the thrust and power are computed and compared with the FAST code [21], which are shown in Figure 10. The results agree with FAST excellently for the wind speeds up to 11.4 m/s. The quantitative values of interest and errors are shown in Table 7.    As for the high wind speed, the local Reynolds number of the blades is calculated at the rated wind speed, as defined below: where c is the chord length; ν is the air kinematic viscosity; U blade is the local velocity, and can be calculated by U blade = U 2 axial + (ωr) 2 , where U axial is the axial velocity; and ω is the rotating angular velocity. The comparison with the result of the BEM method is depicted in Figure 11. It shows that the local Reynolds number varies in different locations of the blade and the maximum value is near 1.2 × 10 7 , which is far bigger than that given in the reference. Thus, the airfoil data become unreliable for the wind speeds greater than 11.4 m/s. In addition, the size of the wind turbine is so large that the dynamic stall should be considered at the high wind speed. The dynamic stall model is not used in this paper. Additionally, the blades need to adjust their pitch angles and maintain the tip speed ratio to maximize the power output when the wind speed exceeds the rated speed. In summary, large wind speeds are not considered in this work.
As for the high wind speed, the local Reynolds number of the blades is calculated at the rated wind speed, as defined below: blade U c Re ν = (29) where c is the chord length; ν is the air kinematic viscosity; blade U is the local velocity, and can be calculated by  Figure 11. It shows that the local Reynolds number varies in different locations of the blade and the maximum value is near 7 1.2 10 × , which is far bigger than that given in the reference. Thus, the airfoil data become unreliable for the wind speeds greater than 11.4 m/s. In addition, the size of the wind turbine is so large that the dynamic stall should be considered at the high wind speed. The dynamic stall model is not used in this paper. Additionally, the blades need to adjust their pitch angles and maintain the tip speed ratio to maximize the power output when the wind speed exceeds the rated speed. In summary, large wind speeds are not considered in this work. To validate the effect of the tip loss model, the comparison of the force distributions with and without the tip loss model is given for U0 = 11.4 m/s (See Figure 12). There is no difference in the middle of the blade, but obvious disagreement is found at the tip if the model of tip loss is not included. Although the discrepancy in the force is small, the residual of power could rise up to 12% (See Table 8). Therefore, it seems to be significant to add the tip loss correction model when using ALM.  To validate the effect of the tip loss model, the comparison of the force distributions with and without the tip loss model is given for U 0 = 11.4 m/s (See Figure 12). There is no difference in the middle of the blade, but obvious disagreement is found at the tip if the model of tip loss is not included. Although the discrepancy in the force is small, the residual of power could rise up to 12% (See Table 8). Therefore, it seems to be significant to add the tip loss correction model when using ALM.
As for the high wind speed, the local Reynolds number of the blades is calculated at the rated wind speed, as defined below: blade U c Re ν = (29) where c is the chord length; ν is the air kinematic viscosity; blade U is the local velocity, and can be calculated by  Figure 11. It shows that the local Reynolds number varies in different locations of the blade and the maximum value is near 7 1.2 10 × , which is far bigger than that given in the reference. Thus, the airfoil data become unreliable for the wind speeds greater than 11.4 m/s. In addition, the size of the wind turbine is so large that the dynamic stall should be considered at the high wind speed. The dynamic stall model is not used in this paper. Additionally, the blades need to adjust their pitch angles and maintain the tip speed ratio to maximize the power output when the wind speed exceeds the rated speed. In summary, large wind speeds are not considered in this work. To validate the effect of the tip loss model, the comparison of the force distributions with and without the tip loss model is given for U0 = 11.4 m/s (See Figure 12). There is no difference in the middle of the blade, but obvious disagreement is found at the tip if the model of tip loss is not included. Although the discrepancy in the force is small, the residual of power could rise up to 12% (See Table 8). Therefore, it seems to be significant to add the tip loss correction model when using ALM.

Wake Structures
The vortex isosurfaces to illustrate the downstream development of the wake vortices are depicted in Figure 13. The wake vortices are identified by the Q criterion, where Q is adjusted as 0.003 for the three cases. It is shown that the tip and hub vortices are clear behind the rotor. They persist for nearly two turns, diffuse at U 0 = 11.4 m/s, and spread out early when the wind speed decreased, corresponding to the increase in the tip speed ratio. Additionally, they become a continuous vortex sheet at U 0 = 5 m/s.

Wake Structures
The vortex isosurfaces to illustrate the downstream development of the wake vortices are depicted in Figure 13. The wake vortices are identified by the Q criterion, where Q is adjusted as 0.003 for the three cases. It is shown that the tip and hub vortices are clear behind the rotor. They persist for nearly two turns, diffuse at U0 = 11.4 m/s, and spread out early when the wind speed decreased, corresponding to the increase in the tip speed ratio. Additionally, they become a continuous vortex sheet at U0 = 5 m/s. The axial flow interference factors represent the influence on the performance of the wind turbine and the distribution of the loading on the blades, which can be calculated by: where x U is the velocity behind the rotor and 0 U is the wind speed. Figure 14 depicts the radial distributions of averaged axial flow interference factors at different wind speeds. From Figure 14, it can be seen that the curves at the wind speed of 5 m/s and 8 m/s are more generally steady than that of 11.4 m/s, but the values are too large under the wind speed of 5 m/s, which means the loading on the blade is not balanced and the influence on the performance is worse than that under 8 m/s.  The axial flow interference factors represent the influence on the performance of the wind turbine and the distribution of the loading on the blades, which can be calculated by: where U x is the velocity behind the rotor and U 0 is the wind speed. Figure 14 depicts the radial distributions of averaged axial flow interference factors at different wind speeds. From Figure 14, it can be seen that the curves at the wind speed of 5 m/s and 8 m/s are more generally steady than that of 11.4 m/s, but the values are too large under the wind speed of 5 m/s, which means the loading on the blade is not balanced and the influence on the performance is worse than that under 8 m/s.

Wake Structures
The vortex isosurfaces to illustrate the downstream development of the wake vortices are depicted in Figure 13. The wake vortices are identified by the Q criterion, where Q is adjusted as 0.003 for the three cases. It is shown that the tip and hub vortices are clear behind the rotor. They persist for nearly two turns, diffuse at U0 = 11.4 m/s, and spread out early when the wind speed decreased, corresponding to the increase in the tip speed ratio. Additionally, they become a continuous vortex sheet at U0 = 5 m/s. The axial flow interference factors represent the influence on the performance of the wind turbine and the distribution of the loading on the blades, which can be calculated by: where x U is the velocity behind the rotor and 0 U is the wind speed. Figure 14 depicts the radial distributions of averaged axial flow interference factors at different wind speeds. From Figure 14, it can be seen that the curves at the wind speed of 5 m/s and 8 m/s are more generally steady than that of 11.4 m/s, but the values are too large under the wind speed of 5 m/s, which means the loading on the blade is not balanced and the influence on the performance is worse than that under 8 m/s.  It is common that the wind turbine turns the wind energy into mechanical energy and the reduction of the wind speed is described as the velocity deficit. It depicts the efficiency of the energy transformation and can affect the performance of the downstream wind turbine. The velocity deficit can be seen through the axial velocity distribution, which is simply shown in Figure 15. It is worthwhile to note that the axial velocities have been normalized with the corresponding wind speed and so the normalized velocities at the upstream station are 1. Figure 15 illustrates the curves of the average axial velocity ratio at the hub height. It is easy to see the velocity deficit in the wake and the velocity is recovering with the growing distances. The location of the maximum velocity deficit is around 1D~3D. With the wind speed becoming greater, it moves away from the rotor and becomes higher. That is as indication that the utilization of the wind energy is getting smaller and the rate of the recovery gets slower.
From the above analysis, it can be seen that the utilization of the wind energy and the ratio of the velocity recovery is larger at the low wind speed, but the loading on the blades is so unsteady and the power is too small. It can thus be simply concluded that the wind turbine can work in an oppositely fine status under the wind speed of 8 m/s. In that case, the power coefficient may become coincidently maxima. Thus, it is reasonable that the wind speed has been set as 8 m/s in many documented studies. For further investigation, some representative numerical results are given at 8 m/s in the next section.

Appl. Sci. 2018, 8, x 15 of 24
It is common that the wind turbine turns the wind energy into mechanical energy and the reduction of the wind speed is described as the velocity deficit. It depicts the efficiency of the energy transformation and can affect the performance of the downstream wind turbine. The velocity deficit can be seen through the axial velocity distribution, which is simply shown in Figure 15. It is worthwhile to note that the axial velocities have been normalized with the corresponding wind speed and so the normalized velocities at the upstream station are 1. Figure 15 illustrates the curves of the average axial velocity ratio at the hub height. It is easy to see the velocity deficit in the wake and the velocity is recovering with the growing distances. The location of the maximum velocity deficit is around 1D~3D. With the wind speed becoming greater, it moves away from the rotor and becomes higher. That is as indication that the utilization of the wind energy is getting smaller and the rate of the recovery gets slower.
From the above analysis, it can be seen that the utilization of the wind energy and the ratio of the velocity recovery is larger at the low wind speed, but the loading on the blades is so unsteady and the power is too small. It can thus be simply concluded that the wind turbine can work in an oppositely fine status under the wind speed of 8 m/s. In that case, the power coefficient may become coincidently maxima. Thus, it is reasonable that the wind speed has been set as 8 m/s in many documented studies. For further investigation, some representative numerical results are given at 8 m/s in the next section.

Results in 8 m/s
The computed power and thrust are 1.9026 MW and 433.47 kN, respectively, and the errors between FAST and these are 4.245% and 5.079%, correspondingly. For more details, the reader can refer to Table 9. Due to the tower's influence on the rotor, the changes in p C and T C are not unilateral but periodical (See Figure 16). It can also be seen that the wind turbine can obtain a relatively stable status. Furthermore, the interval between the two extreme points is about 2.2 s. That is when the blades rotate nearly 120 degrees, which is one-third of the cycle of the rotor revolving.

Results in 8 m/s
The computed power and thrust are 1.9026 MW and 433.47 kN, respectively, and the errors between FAST and these are 4.245% and 5.079%, correspondingly. For more details, the reader can refer to Table 9. Due to the tower's influence on the rotor, the changes in C p and C T are not unilateral but periodical (See Figure 16). It can also be seen that the wind turbine can obtain a relatively stable status. Furthermore, the interval between the two extreme points is about 2.2 s. That is when the blades rotate nearly 120 degrees, which is one-third of the cycle of the rotor revolving.

Appl. Sci. 2018, 8, x 15 of 24
It is common that the wind turbine turns the wind energy into mechanical energy and the reduction of the wind speed is described as the velocity deficit. It depicts the efficiency of the energy transformation and can affect the performance of the downstream wind turbine. The velocity deficit can be seen through the axial velocity distribution, which is simply shown in Figure 15. It is worthwhile to note that the axial velocities have been normalized with the corresponding wind speed and so the normalized velocities at the upstream station are 1. Figure 15 illustrates the curves of the average axial velocity ratio at the hub height. It is easy to see the velocity deficit in the wake and the velocity is recovering with the growing distances. The location of the maximum velocity deficit is around 1D~3D. With the wind speed becoming greater, it moves away from the rotor and becomes higher. That is as indication that the utilization of the wind energy is getting smaller and the rate of the recovery gets slower.
From the above analysis, it can be seen that the utilization of the wind energy and the ratio of the velocity recovery is larger at the low wind speed, but the loading on the blades is so unsteady and the power is too small. It can thus be simply concluded that the wind turbine can work in an oppositely fine status under the wind speed of 8 m/s. In that case, the power coefficient may become coincidently maxima. Thus, it is reasonable that the wind speed has been set as 8 m/s in many documented studies. For further investigation, some representative numerical results are given at 8 m/s in the next section.

Results in 8 m/s
The computed power and thrust are 1.9026 MW and 433.47 kN, respectively, and the errors between FAST and these are 4.245% and 5.079%, correspondingly. For more details, the reader can refer to Table 9. Due to the tower's influence on the rotor, the changes in p C and T C are not unilateral but periodical (See Figure 16). It can also be seen that the wind turbine can obtain a relatively stable status. Furthermore, the interval between the two extreme points is about 2.2 s. That is when the blades rotate nearly 120 degrees, which is one-third of the cycle of the rotor revolving.    Figure 17 shows the spanwise distribution of the time-averaged normalized axial velocity at the streamwise positions of −0.5D, 0.5D, 1D, 2D, 3D, 5D, 6D, 7D, and 8D, which represent the development of the velocity in the wake. There is a clear phenomenon of the turbine extracting momentum from the incoming flow when the airflow passes through the wind turbine. That is the velocity deficit. The asymmetrical velocity profile is observed and the velocity value is not 1.0 at y/D = 0, when the effects of the tower and the hub are included. The velocity profile changes from a single peak line to double peaks when the wind flows past the rotor, and gradually becomes a single peak again. The double peak line represents that the velocity deficit is large in the middle of the blade and small in the hub, which is an indication that the loading is focused on the middle of the blade. The rate of velocity recovery becomes slower with the increasing distance, which means that the influence of the wake becomes weaker. More details about the wake development can be seen in the research of Mo [31]. It seems that the wake effect is little from 5D to 8D, and further study is presented in the next section.   Figure 17 shows the spanwise distribution of the time-averaged normalized axial velocity at the streamwise positions of −0.5D, 0.5D, 1D, 2D, 3D, 5D, 6D, 7D, and 8D, which represent the development of the velocity in the wake. There is a clear phenomenon of the turbine extracting momentum from the incoming flow when the airflow passes through the wind turbine. That is the velocity deficit. The asymmetrical velocity profile is observed and the velocity value is not 1.0 at y/D = 0, when the effects of the tower and the hub are included. The velocity profile changes from a single peak line to double peaks when the wind flows past the rotor, and gradually becomes a single peak again. The double peak line represents that the velocity deficit is large in the middle of the blade and small in the hub, which is an indication that the loading is focused on the middle of the blade. The rate of velocity recovery becomes slower with the increasing distance, which means that the influence of the wake becomes weaker. More details about the wake development can be seen in the research of Mo [31]. It seems that the wake effect is little from 5D to 8D, and further study is presented in the next section.

Double NREL 5-MW Wind Turbines Case
In this section, two NREL 5-MW wind turbines that are aligned together are simulated. The upstream wind turbine is set as WT1 and the downstream one is set as WT2. The best tip speed ratio of WT2 is calculated and the best separation distance is simply proved to be 7D. Finally, the influence of WT1 wake on WT2 is discussed.

Numerical Model
For the two wind turbines, the length of the wake flow field should be longer and thus more computational time is needed. For this, the size of the computation domain is extended to 2772 m × 1260 m × 468 m. The details can be seen in Figure 18. The mesh refinement is the same as that shown

Double NREL 5-MW Wind Turbines Case
In this section, two NREL 5-MW wind turbines that are aligned together are simulated. The upstream wind turbine is set as WT1 and the downstream one is set as WT2. The best tip speed ratio of WT2 is calculated and the best separation distance is simply proved to be 7D. Finally, the influence of WT1 wake on WT2 is discussed.

Numerical Model
For the two wind turbines, the length of the wake flow field should be longer and thus more computational time is needed. For this, the size of the computation domain is extended to 2772 m × 1260 m × 468 m. The details can be seen in Figure 18. The mesh refinement is the same as that shown in Section 4.1. The size of the outer mesh is 24 m and the inner zone grid size is 3 m.
Moreover, the refined mesh in domain B is also similar to A. The total number of the cell is 4,722,712. The parameters of calculation are similar to those in Section 4.1. According to the conclusion above, the wind speed is set as 8 m/s and the tip loss model is changed into the Shen model. It is vital to observe the wake that is fully developed, so the simulation is calculated for 40 revolutions of WT1, which is around 262 s. Due to the lengthy time of simulation, the time step should be much bigger. The mesh independence tests have been presented previously in Section 4.2.
Appl. Sci. 2018, 8, x 17 of 24 in Section 4.1. The size of the outer mesh is 24 m and the inner zone grid size is 3 m. Moreover, the refined mesh in domain B is also similar to A. The total number of the cell is 4,722,712. The parameters of calculation are similar to those in Section 4.1. According to the conclusion above, the wind speed is set as 8 m/s and the tip loss model is changed into the Shen model. It is vital to observe the wake that is fully developed, so the simulation is calculated for 40 revolutions of WT1, which is around 262 s. Due to the lengthy time of simulation, the time step should be much bigger. The mesh independence tests have been presented previously in Section 4.2.

Effect of Different WT2 Tip Speed Ratio
Under the wind speed of 8 m/s, the tip speed ratio of WT1 is 7.55, but this is not true for WT2 because of the WT1 wake. Figure 19 depicts the situation of the power coefficient with different tip speed ratios of WT2. The power coefficient of WT1 shows almost no changes. However, the output power of WT2 varies a lot, and there is a distinct maximum in the range from 6.0 to 6.5. In one word, the tip speed ratio of WT2 is set as 6.1, and details of the power coefficient with different tip speed ratios are shown in Table 10.

Effect of Different WT2 Tip Speed Ratio
Under the wind speed of 8 m/s, the tip speed ratio of WT1 is 7.55, but this is not true for WT2 because of the WT1 wake. Figure 19 depicts the situation of the power coefficient with different tip speed ratios of WT2. The power coefficient of WT1 shows almost no changes. However, the output power of WT2 varies a lot, and there is a distinct maximum in the range from 6.0 to 6.5. In one word, the tip speed ratio of WT2 is set as 6.1, and details of the power coefficient with different tip speed ratios are shown in Table 10.

Effect of Different WT2 Tip Speed Ratio
Under the wind speed of 8 m/s, the tip speed ratio of WT1 is 7.55, but this is not true for WT2 because of the WT1 wake. Figure 19 depicts the situation of the power coefficient with different tip speed ratios of WT2. The power coefficient of WT1 shows almost no changes. However, the output power of WT2 varies a lot, and there is a distinct maximum in the range from 6.0 to 6.5. In one word, the tip speed ratio of WT2 is set as 6.1, and details of the power coefficient with different tip speed ratios are shown in Table 10.

Effect of Separation
In many studies, the separation of two wind turbines is usually defined 7D as a default. From Section 4.4, it could be found that the separation distance should be 6D, 7D, or even 8D. In addition, here, a parameter is set to represent the ability of production in the whole wind farm as: where P is the power of a single downstream wind turbine and d is the average distance between the two wind turbines. The wind speed is assumed to be a constant and this is also true for the production of the downstream wind turbine. Table 11 gives the influence for the distance from 6D to 8D on the power production. It can be seen that the mean power of WT1 is almost the same, which means that the power of the upstream wind turbine is not influenced by the downstream one. Additionally, the ability of power production has little difference in either the 7D or 8D case. It is worthwhile to note that the cost of the power transmission is expensive and the loss of electricity transportation is large with the distance. Hence, this distance should be as short as possible and 7D is chosen as the realistic value.

Result of 7D Separation
Although changing the tip loss model from Glauert to Shen and other adjustments make the errors of the mean power larger, it is still with the acceptable range. The comparison of the power with improved ALM by Jha et al.
[32] is presented in Table 12, in which the power of WT2 reduces a lot and is only 39.8% of WT1 due to the influence of the WT1 wake effect. The power and thrust coefficients of WT1 and WT2 are described in Figures 20 and 21. The tendency of the WT1 power and thrust curves shows almost no changes. Though the curves of WT2 continue to become lower, the reduction is no more than 2%. Besides, the interval between the two extremum points in WT2 is bigger than WT1. The reason for this is due to the fact that the rotational speed of WT2 is smaller than that of WT1. As can be seen in Figures 20 and 21, the curves between the two extremum points fluctuate significantly. It proves that the WT2 is located in the wake flow field of WT1 and subjected to the flow unsteadiness. extremum points in WT2 is bigger than WT1. The reason for this is due to the fact that the rotational speed of WT2 is smaller than that of WT1. As can be seen in Figures 20 and 21, the curves between the two extremum points fluctuate significantly. It proves that the WT2 is located in the wake flow field of WT1 and subjected to the flow unsteadiness.   Figure 22 shows the radial distributions of average tangential and normal forces of WT1 and WT2. The distribution patterns of WT2 are very similar to those of WT1, and the force amplitude of WT2 is nearly 56% of WT1. It is noted that Figure 22 shows the averaged value of the time but not the instantaneous one, which may not be sufficient to establish a sound conclusion. Considering this, five-points data on the blade are picked up in different times, and their standard deviations are calculated. The outcomes are shown in Figure 23. It can be observed that the standard deviations of WT2 are greater than those of WT1, especially in the case of normal force. This is an indication that the force varies widely, and the forces get their maximums from r/R = 0.6 to r/R = 0.8. Consequently, the loading distribution of WT2 is lower than that of WT1 across the range. According to the discussion above, the distribution of the forces on the blades and the stabilization of the force both affect the production of the wind turbine. Because of the power of WT1 extremum points in WT2 is bigger than WT1. The reason for this is due to the fact that the rotational speed of WT2 is smaller than that of WT1. As can be seen in Figures 20 and 21, the curves between the two extremum points fluctuate significantly. It proves that the WT2 is located in the wake flow field of WT1 and subjected to the flow unsteadiness.   Figure 22 shows the radial distributions of average tangential and normal forces of WT1 and WT2. The distribution patterns of WT2 are very similar to those of WT1, and the force amplitude of WT2 is nearly 56% of WT1. It is noted that Figure 22 shows the averaged value of the time but not the instantaneous one, which may not be sufficient to establish a sound conclusion. Considering this, five-points data on the blade are picked up in different times, and their standard deviations are calculated. The outcomes are shown in Figure 23. It can be observed that the standard deviations of WT2 are greater than those of WT1, especially in the case of normal force. This is an indication that the force varies widely, and the forces get their maximums from r/R = 0.6 to r/R = 0.8. Consequently, the loading distribution of WT2 is lower than that of WT1 across the range. According to the discussion above, the distribution of the forces on the blades and the stabilization of the force both affect the production of the wind turbine. Because of the power of WT1  Figure 22 shows the radial distributions of average tangential and normal forces of WT1 and WT2. The distribution patterns of WT2 are very similar to those of WT1, and the force amplitude of WT2 is nearly 56% of WT1. It is noted that Figure 22 shows the averaged value of the time but not the instantaneous one, which may not be sufficient to establish a sound conclusion. Considering this, five-points data on the blade are picked up in different times, and their standard deviations are calculated. The outcomes are shown in Figure 23. It can be observed that the standard deviations of WT2 are greater than those of WT1, especially in the case of normal force. This is an indication that the force varies widely, and the forces get their maximums from r/R = 0.6 to r/R = 0.8. Consequently, the loading distribution of WT2 is lower than that of WT1 across the range. extremum points in WT2 is bigger than WT1. The reason for this is due to the fact that the rotational speed of WT2 is smaller than that of WT1. As can be seen in Figures 20 and 21, the curves between the two extremum points fluctuate significantly. It proves that the WT2 is located in the wake flow field of WT1 and subjected to the flow unsteadiness.   Figure 22 shows the radial distributions of average tangential and normal forces of WT1 and WT2. The distribution patterns of WT2 are very similar to those of WT1, and the force amplitude of WT2 is nearly 56% of WT1. It is noted that Figure 22 shows the averaged value of the time but not the instantaneous one, which may not be sufficient to establish a sound conclusion. Considering this, five-points data on the blade are picked up in different times, and their standard deviations are calculated. The outcomes are shown in Figure 23. It can be observed that the standard deviations of WT2 are greater than those of WT1, especially in the case of normal force. This is an indication that the force varies widely, and the forces get their maximums from r/R = 0.6 to r/R = 0.8. Consequently, the loading distribution of WT2 is lower than that of WT1 across the range. According to the discussion above, the distribution of the forces on the blades and the stabilization of the force both affect the production of the wind turbine. Because of the power of WT1 According to the discussion above, the distribution of the forces on the blades and the stabilization of the force both affect the production of the wind turbine. Because of the power of WT1 being higher, the stabilization effect is more crucial than the distribution one. Besides, the distribution is more highly concentrated and the range of change is larger in the middle of the blades, where the blades are easily damaged and need to be paid much more attention during the design. being higher, the stabilization effect is more crucial than the distribution one. Besides, the distribution is more highly concentrated and the range of change is larger in the middle of the blades, where the blades are easily damaged and need to be paid much more attention during the design. Furthermore, Figure 24 depicts the radial distributions of averaged axial flow interference factors. It can be seen that the WT1 distribution of a is similar to that in Section 4.3. The tendency of a is almost the same between WT1 and WT2, but the value is much bigger in WT2. This means that the performance of WT2 is worse than that of WT1, in spite of the fact that the distribution of force on the blades is similar. To disclose the primary trends of wake-to-wake velocity effect and study the ability of the wind turbine to extract useful momentum from the incoming flow, it is important to analyze the velocity profiles in the wake. The distribution of the mean axial velocity ratio at the hub height is given in Figure 25. From the figure, there are two pronounced velocity deficits, and the first maximum velocity deficit is bigger than the second one. That caused the quantities of momentum extraction by the rotor of WT1 from the incoming wind to be higher and thus its power is larger. From 2D to 6D, the velocity ratio increased by 0.2, and it grew nearly 0.25 from 10D to 14D. Thus, the velocities recovered much faster behind WT2 rather than WT1. This phenomenon may be supposed to be related to the structure of the complex overlapped vortices. After 14D, the rate of velocity recovery is then slower, but the velocity is below 8 m/s, and this should be caused by the arrival upstream wake combined with the wake of WT2.
Finally, to evaluate the effects of the velocity recovery in the wake region, the spanwise distributions of the average axial normalized velocity at different locations are shown in Figure 26. Before 5D, the distribution of velocity in the wake flow field is similar to that found in Section 4.4, where the distinct velocity deficit can be observed and the velocity profiles change from a double peak into a single peak. After the wind flow passing through WT2, there is another velocity deficit, Furthermore, Figure 24 depicts the radial distributions of averaged axial flow interference factors. It can be seen that the WT1 distribution of a is similar to that in Section 4.3. The tendency of a is almost the same between WT1 and WT2, but the value is much bigger in WT2. This means that the performance of WT2 is worse than that of WT1, in spite of the fact that the distribution of force on the blades is similar. being higher, the stabilization effect is more crucial than the distribution one. Besides, the distribution is more highly concentrated and the range of change is larger in the middle of the blades, where the blades are easily damaged and need to be paid much more attention during the design. Furthermore, Figure 24 depicts the radial distributions of averaged axial flow interference factors. It can be seen that the WT1 distribution of a is similar to that in Section 4.3. The tendency of a is almost the same between WT1 and WT2, but the value is much bigger in WT2. This means that the performance of WT2 is worse than that of WT1, in spite of the fact that the distribution of force on the blades is similar. To disclose the primary trends of wake-to-wake velocity effect and study the ability of the wind turbine to extract useful momentum from the incoming flow, it is important to analyze the velocity profiles in the wake. The distribution of the mean axial velocity ratio at the hub height is given in Figure 25. From the figure, there are two pronounced velocity deficits, and the first maximum velocity deficit is bigger than the second one. That caused the quantities of momentum extraction by the rotor of WT1 from the incoming wind to be higher and thus its power is larger. From 2D to 6D, the velocity ratio increased by 0.2, and it grew nearly 0.25 from 10D to 14D. Thus, the velocities recovered much faster behind WT2 rather than WT1. This phenomenon may be supposed to be related to the structure of the complex overlapped vortices. After 14D, the rate of velocity recovery is then slower, but the velocity is below 8 m/s, and this should be caused by the arrival upstream wake combined with the wake of WT2.
Finally, to evaluate the effects of the velocity recovery in the wake region, the spanwise distributions of the average axial normalized velocity at different locations are shown in Figure 26. Before 5D, the distribution of velocity in the wake flow field is similar to that found in Section 4.4, where the distinct velocity deficit can be observed and the velocity profiles change from a double peak into a single peak. After the wind flow passing through WT2, there is another velocity deficit, To disclose the primary trends of wake-to-wake velocity effect and study the ability of the wind turbine to extract useful momentum from the incoming flow, it is important to analyze the velocity profiles in the wake. The distribution of the mean axial velocity ratio at the hub height is given in Figure 25. From the figure, there are two pronounced velocity deficits, and the first maximum velocity deficit is bigger than the second one. That caused the quantities of momentum extraction by the rotor of WT1 from the incoming wind to be higher and thus its power is larger. From 2D to 6D, the velocity ratio increased by 0.2, and it grew nearly 0.25 from 10D to 14D. Thus, the velocities recovered much faster behind WT2 rather than WT1. This phenomenon may be supposed to be related to the structure of the complex overlapped vortices. After 14D, the rate of velocity recovery is then slower, but the velocity is below 8 m/s, and this should be caused by the arrival upstream wake combined with the wake of WT2.
Finally, to evaluate the effects of the velocity recovery in the wake region, the spanwise distributions of the average axial normalized velocity at different locations are shown in Figure 26. Before 5D, the distribution of velocity in the wake flow field is similar to that found in Section 4.4, where the distinct velocity deficit can be observed and the velocity profiles change from a double peak into a single peak. After the wind flow passing through WT2, there is another velocity deficit, which is more serious than before. From 8D to 10D, the magnitude of the velocity deficit is approximately 0.3, but it is only about 0.15 from 1D to 3D. That quantitatively proves that the level of the velocity recovery behind WT2 is higher again. At the region after 11D, the velocity recovers only a little and the value of the velocity ratio is not in unity. That is mainly due to the overlap of the upstream wake from WT1 and the wake generated by WT2. which is more serious than before. From 8D to 10D, the magnitude of the velocity deficit is approximately 0.3, but it is only about 0.15 from 1D to 3D. That quantitatively proves that the level of the velocity recovery behind WT2 is higher again. At the region after 11D, the velocity recovers only a little and the value of the velocity ratio is not in unity. That is mainly due to the overlap of the upstream wake from WT1 and the wake generated by WT2. Overall, the existence of WT2 has no influence on the performance of WT1. As a result, it can be concluded that only the upstream wind turbine can work at the designed status, and the downstream wind turbine cannot achieve this due to the influence of the wake, which has been generated by either a single or multiple upstream wind turbines. Due to the effect of the wake generated by the upstream wind turbine, the loading on the blades of the downstream one is more unsteady and as a result, its fatigue damage deserves further study.  Overall, the existence of WT2 has no influence on the performance of WT1. As a result, it can be concluded that only the upstream wind turbine can work at the designed status, and the downstream wind turbine cannot achieve this due to the influence of the wake, which has been generated by either a single or multiple upstream wind turbines. Due to the effect of the wake generated by the upstream wind turbine, the loading on the blades of the downstream one is more unsteady and as a result, its fatigue damage deserves further study. which is more serious than before. From 8D to 10D, the magnitude of the velocity deficit is approximately 0.3, but it is only about 0.15 from 1D to 3D. That quantitatively proves that the level of the velocity recovery behind WT2 is higher again. At the region after 11D, the velocity recovers only a little and the value of the velocity ratio is not in unity. That is mainly due to the overlap of the upstream wake from WT1 and the wake generated by WT2. Overall, the existence of WT2 has no influence on the performance of WT1. As a result, it can be concluded that only the upstream wind turbine can work at the designed status, and the downstream wind turbine cannot achieve this due to the influence of the wake, which has been generated by either a single or multiple upstream wind turbines. Due to the effect of the wake generated by the upstream wind turbine, the loading on the blades of the downstream one is more unsteady and as a result, its fatigue damage deserves further study.

Conclusions and Discussions
The purpose of the present research is to investigate the harmful effect caused by the wake flow, so a single wind turbine and double wind turbines aligned in a line are modeled using ALM. The major results are summarized as follows: 1.
It is reasonable to set the wind speed as 8 m/s. The NREL 5-MW wind turbine could obtain a better working status at this speed, the damage to the blades is minimal, and its wake recovery is fast.

2.
It is proved through the velocity deficit that the separation distance of two turbines should not be less than 5D, because of the great influence of the wake. This conclusion is the same as that found by Choi et al. [33]. A more in-depth investigation found that the best separation is nearly 7D. 3.
The upstream wind turbine has a great effect on the downstream one. Due to the wake generated by the upstream wind turbine, the output power and the loading of the downstream one are heavily influenced.
ALM is used in the model, and the result is shown to be accurate enough, but some limitations exist. For example, the elasticity of the blades is not considered adequately and added into the fully coupled system of the floating offshore wind turbine. Further study is required to refine the model and make it more widely acceptable.