Aeroelastic Performance Analysis of Wind Turbine in the Wake with a New Elastic Actuator Line Model

: The scale of a wind turbine is getting larger with the development of wind energy recently. Therefore, the e ﬀ ect of the wind turbine blades deformation on its performances and lifespan has become obvious. In order to solve this research rapidly, a new elastic actuator line model (EALM) is proposed in this study, which is based on turbinesFoam in OpenFOAM (Open Source Field Operation and Manipulation, a free, open source computational ﬂuid dynamics (CFD) software package released by the OpenFOAM Foundation, which was incorporated as a company limited by guarantee in England and Wales). The model combines the actuator line model (ALM) and a beam solver, which is used in the wind turbine blade design. The aeroelastic performances of the NREL (National Renewable Energy Laboratory) 5 MW wind turbine like power, thrust, and blade tip displacement are investigated. These results are compared with some research to prove the new model. Additionally, the inﬂuence caused by blade deﬂections on the aerodynamic performance is discussed. It is demonstrated that the tower shadow e ﬀ ect becomes more obvious and causes the power and thrust to get a bit lower and unsteady. Finally, this variety is analyzed in the wake of upstream wind turbine and it is found that the inﬂuence on the performance and wake ﬂow ﬁeld of downstream wind turbine becomes more serious.


Introduction
With the improvement of wind power technology and the demand of high-power generation, the target of wind turbine design turns to large scale and offshore [1][2][3][4][5][6][7]. In 2009, the National Renewable Energy Laboratory (NREL) in America defined a 5 MW reference wind turbine for offshore system, in which the rotor diameter is 126 m [8]. Technical University of Denmark described a 10 MW reference wind turbine whose rotor diameter is 178.3 m [9] in 2013. As the blade of wind turbine gets longer, it becomes less stiff, more susceptible, and easily deformable, which will lead to increased fatigue damage and reduced production. When simulating the aeroelastic performances of a floating offshore wind turbine or a wind farm, there are many challenges to solve.
There are mainly three methods to study the aerodynamic performances of wind turbines. The first one is the Blade Element Momentum (BEM) theory, which combines the blade element theory and momentum theory. It has high efficiency and is widely used in the industrial application, but the information of flow field is not considered. The second one is the Computational Fluid Dynamics (CFD) method, which calculates the velocity and pressure fields by solving the Navier-Stokes equations. It can obtain quite accurate results and is gradually used in recent years, e.g., Tran et al. [10] analyzed In the actuator line model, the blades of wind turbine are replaced by lines with the aerodynamic force distributed on it, as shown in Figure 1. The lines are discretized into actuator sections, in which there are different airfoils, chords and structure twists. The loading in each section, which is usually named as the aerodynamic force, can be calculated by Equation (1).
where f is the aerodynamic force; L and D are lift and drag, respectively; ρ is the air density, and ρ = 1.225 kg/m 3 in this paper; U is the relative velocity of the blade element in actuator sections; c is the chord length; C l and C d are the coefficients of lift and drag, which can be looked up from airfoil data tables derived from physical experiment; → e L and → e D represent the unit direction vectors of L and D. Besides, the aerodynamic loads are concentrated on the aerodynamic center of the blade element, which are located in 0.25 chord length. To describe the influence of wind turbine blade on the flow field, the aerodynamic loads should be dispersed on the grid point. There are many distributional ways, and three-dimensional Gaussian distribution η ε is used in this paper, whose function is shown in Equation (2). This equation is a smooth function and can avoid numerical singularity.
where ε is a constant to adjust the strength of the distribution formula, and it is two times of local grid scale in this study according to Shives [26]; and d is the distance between the force point and the grid point. Therefore, the forces f ε on the nearby mesh can be calculated by Equation (3) and added into the momentum equation to solve the flow field as the body forces.
∂V ∂t where V and p are the wind speed and pressure in the flow; υ is the kinematic viscosity, and it is set as 1.5 × 10 −5 m 2 /s ; f ε is the source item calculated through Equation (3). Moreover, the Prandtl-Glauert tip loss model is used to the consideration that the velocity is zero at the tip. The function of tip loss model is shown as Equation (5).
where F tip is the aerodynamic corrected force at the tip; B is the blades number; R is the radius of the blades; r is the distance between the root of blades and the action point of aerodynamic force; and φ is the angle between U and the rotor plane [27]. Tower and hub effect are also considered in this research by similar technology, and their airfoils are cylinder.

Rotating Beam Solver
The blades of wind turbine are considered as rotating variable cross-section projecting beams in this study. In this section, an approach of solving these beams in actual design is introduced [24]. Although this method is faulty and not suitable for structural solution, it still could catch some phenomena in vibration. Figure 2 shows the diagram of the blade deflection, in which the blades have been substituted by actuator lines. In the picture, direction 0 is out of the rotating plane and called flapwise. Similarly, direction 1 is in the rotating plane and called edgewise. In this research, the subscript 0 means the component of the physical quantity in 0 direction and the subscript 1 represents component of the corresponding parameter in 1 direction. The diagrammatic sketch of the forces on the blade, which is considered as a cantilever beam, is given in Figure 3. According to the Newton's second law of a microelement, the sheer stresses 0 T , 1 T and the bending moments can be calculated from Equation (6) to Equation (9).
where p(x) is the aerodynamic force distribution; m(x) is the blade mass distribution; g is the gravitational acceleration, and it is set as 2 9.8 / ms in this paper;  is the shaft tilt angle, and it is 5 degree here; y is the displacement y to the second derivative of time, and it can be calculated by prediction and correction of modified Euler method, given from Equation (10) to Equation (13); N0, N1 are the component of centrifugal force in 0 direction and 1 direction. 11 2 n n n n y y y y dtdt

Rotating Beam Solver
The blades of wind turbine are considered as rotating variable cross-section projecting beams in this study. In this section, an approach of solving these beams in actual design is introduced [24]. Although this method is faulty and not suitable for structural solution, it still could catch some phenomena in vibration. Figure 2 shows the diagram of the blade deflection, in which the blades have been substituted by actuator lines. In the picture, direction 0 is out of the rotating plane and called flapwise. Similarly, direction 1 is in the rotating plane and called edgewise. In this research, the subscript 0 means the component of the physical quantity in 0 direction and the subscript 1 represents component of the corresponding parameter in 1 direction. The diagrammatic sketch of the forces on the blade, which is considered as a cantilever beam, is given in Figure 3. According to the Newton's second law of a microelement, the sheer stresses T 0 , T 1 and the bending moments M 0 , M 1 of each blade element can be calculated from Equation (6) to Equation (9). .. ..
where p(x) is the aerodynamic force distribution; m(x) is the blade mass distribution; g is the gravitational acceleration, and it is set as 9.8 m/s 2 in this paper; θ is the shaft tilt angle, and it is 5 degree here; ..
y is the displacement y to the second derivative of time, and it can be calculated by prediction and correction of modified Euler method, given from Equation (10) to Equation (13); N 0 , N 1 are the component of centrifugal force in 0 direction and 1 direction. ..
where the subscript n + 1 represents at time t n+1 ; the subscript n represents at time t n ; y is predicted value; function h can be calculated from transposition of structure equations, as shown in Equation (15).
where Equation (14) is the structure equations without damping; M(x) is mass of blade; K(x) is the stiffness of blade; F(x) is the external loading. According to Figure 4, the bending moments M 0 and M 1 can be transformed into principal axes direction by Equation (16) and Equation (17).
where M 11 and M 22 are the bending moments on the first principal axis and the second principal axis; β is the structure twist angel.
Water 2020, 12, x FOR PEER REVIEW 6 of 24

New Elastic Actuator Line Model
In this section, the ALM is improved into a new elastic actuator line model (EALM) based on turbinesFoam library, which is developed by Bachant et al. [25]. This library uses ALM to simulate wind and marine hydrokinetic turbines in OpenFOAM. Its interpolation, Gaussian projection, and vector rotation functions are all adapted from NREL's Simulator for Off/Onshore Wind Farm Applications (SOWFA).
In EALM, the body forces are calculated by traditional ALM and the blade deflection is computed by rotating beam solver, which is defined in Section 2.2. The computation process of EALM is given in Figure 5, in which the difference between EALM and ALM is marked by dashed rectangle. The part of structure solver is added into the actuatorLineSource class, and the actuator point is changed in the actuatorLineElement class. From Figure 5, it can be found that this combination of ALM and structure model is one-way coupling. Compared with the research by Meng et al. [23], the part of aerodynamics solver is similar, but the way of dealing with structural solver is different (see Section 2.2). The advantages of this model are that the EALM allows large time step when the position of actuator line changes every time and computes long terms of wind Based on the beam theory, the curvature of principal axes κ 11 and κ 22 can be obtained from Equation (18) and Equation (19).
where EI is the stiffness of the blade element. Equation (18) and Equation (19) are converted back to direction 0 and direction 1 through the formulas as following: Therefore, the angular deformation θ and deflection y can be calculated by: The root of the blade is clamped, so the boundary condition at root is: The tip is free end, so the boundary condition at tip is:

New Elastic Actuator Line Model
In this section, the ALM is improved into a new elastic actuator line model (EALM) based on turbinesFoam library, which is developed by Bachant et al. [25]. This library uses ALM to simulate wind and marine hydrokinetic turbines in OpenFOAM. Its interpolation, Gaussian projection, and vector rotation functions are all adapted from NREL's Simulator for Off/Onshore Wind Farm Applications (SOWFA).
In EALM, the body forces are calculated by traditional ALM and the blade deflection is computed by rotating beam solver, which is defined in Section 2.2. The computation process of EALM is given in Figure 5, in which the difference between EALM and ALM is marked by dashed rectangle. The part of structure solver is added into the actuatorLineSource class, and the actuator point is changed in the actuatorLineElement class. From Figure 5, it can be found that this combination of ALM and structure model is one-way coupling. Compared with the research by Meng et al. [23], the part of aerodynamics solver is similar, but the way of dealing with structural solver is different (see Section 2.2). The advantages of this model are that the EALM allows large time step when the position of actuator line changes every time and computes long terms of wind turbine working. This technology will be further improved and used in the simulation of the floating offshore wind turbine in the sea, which needs more simulation time to keep the floating foundation stability under waves.
Water 2020, 12, x FOR PEER REVIEW 7 of 24 turbine working. This technology will be further improved and used in the simulation of the floating offshore wind turbine in the sea, which needs more simulation time to keep the floating foundation stability under waves.

Computational Wind Turbine Model
In this study, the aerodynamic performances, which are power and thrust, and the structural responses, which are represented mainly by the blade tip displacement, will be compared with different research using varieties methods to validate EALM. According to the theory above, all these three physical quantities are obtained from the aerodynamic forces along the blade. Different cases use varieties method to achieve these forces. Only if the blade properties, including the aerodynamic and structural properties, and the computational condition such as wind speed are all the same, these three quantities could make equivalent comparisons. Thus, the computational condition will be given next. The model used in this research is NREL 5 MW wind turbine, and its gross properties are listed in Table 1. The details of the blade structural and aerodynamic properties can be referenced in Jonkman's technical report [8]. To be exact, the shaft tilt angle is considered in the structural part, which will make the external force of blade increased because of the gravity component.
The simulation physical domains used in this paper are the same as Yu et al. [28]. Their three-dimensional sketches are shown in Figures 6 and 7. In these pictures, A is the outer mesh region and B is the refined mesh region, where the grid is refined by 4 levels. The mesh refined ratio is 0.5 between each level. The mesh independence test has been completed in that research and 1.5 m is chosen as the minimum size of grid. Besides, it is confirmed that the wind turbine gets a better working status at the wind speed of 8 m/s and its corresponding tip speed ratio is 7.55. As a consequence of this model, the cases studied in present research are mainly in two situations, 8 m/s

Computational Wind Turbine Model
In this study, the aerodynamic performances, which are power and thrust, and the structural responses, which are represented mainly by the blade tip displacement, will be compared with different research using varieties methods to validate EALM. According to the theory above, all these three physical quantities are obtained from the aerodynamic forces along the blade. Different cases use varieties method to achieve these forces. Only if the blade properties, including the aerodynamic and structural properties, and the computational condition such as wind speed are all the same, these three quantities could make equivalent comparisons. Thus, the computational condition will be given next. The model used in this research is NREL 5 MW wind turbine, and its gross properties are listed in Table 1. The details of the blade structural and aerodynamic properties can be referenced in Jonkman's technical report [8]. To be exact, the shaft tilt angle is considered in the structural part, which will make the external force of blade increased because of the gravity component. The simulation physical domains used in this paper are the same as Yu et al. [28]. Their three-dimensional sketches are shown in Figures 6 and 7. In these pictures, A is the outer mesh region and B is the refined mesh region, where the grid is refined by 4 levels. The mesh refined ratio is 0.5 between each level. The mesh independence test has been completed in that research and 1.5 m is chosen as the minimum size of grid. Besides, it is confirmed that the wind turbine gets a better working status at the wind speed of 8 m/s and its corresponding tip speed ratio is 7.55. As a consequence of this model, the cases studied in present research are mainly in two situations, 8 m/s and 11.4 m/s. To analyze the aeroelastic performance in the wake, double NREL 5 MW wind turbines set in a line are studied, and its simulation domain is shown in Figure 7. The main contents of this research are focused on the downstream wind turbine.
Compared to the work of Yu et al. [28], the parameter settings are all the same expect that the LES (Large Eddy Simulation) model is used instead of RANS (Reynolds Average Navier-Stokes). That is because the LES model will get more accurate results about the vorticities than RANS and the influence of the wake flows to elastic blade is studied in this research. In addition, the wind turbine blade deflection solver is added in ALM, named EALM.
Water 2020, 12, x FOR PEER REVIEW 8 of 24 Compared to the work of Yu et al. [28], the parameter settings are all the same expect that the LES (Large Eddy Simulation) model is used instead of RANS (Reynolds Average Navier-Stokes). That is because the LES model will get more accurate results about the vorticities than RANS and the influence of the wake flows to elastic blade is studied in this research. In addition, the wind turbine blade deflection solver is added in ALM, named EALM.

Verification and Analysis
In this section, the mesh independence and uncertainty analysis of actuator point number are given first. Then, the aerodynamic performance and the structural responses are compared to validate EALM. The aerodynamic performance contains the power and the thrust. The structural responses are represented mainly by the blade tip displacement.

Mesh Independence and Uncertainty Analysis
Four levels of grids are tested to prove the mesh independence under the rated wind speed, and the non-dimension coefficients of power and thrust are compared in Table 2. They are defined as power coefficient C p and thrust coefficient C t as below: where P is the mechanical power; ρ is the air density; V is the wind speed in the flow; T is the thrust on the blades. It shows that the results drop slow when the grid levels up especially from level 3 to level 4. Therefore, the level 3 mesh is used in the following research. According to Shives' research [26], the number of actuator point has the rule that the maximum of the distance between the adjacent point should not more than the size of blade region. Different actuator point numbers are also tested in Table 3 and Figure 8. Besides, the mesh size and computational condition are kept as constant during this test. To improve the accuracy the uncertainty coefficient U ξ is calculated as follow steps.
(1) Calculate the difference of aerodynamic coefficients between neighbor level, which are defined as ε ζ21 and ε ζ32 ; (2) The convergence ratio R ξ can be computed by Equation (30).
(3) The power coefficient convergence ratio is 0.125 and the thrust coefficient convergence ratio is 0.167. They are both located in the interval which greater than 0 and less than 1. So, the Richardson extrapolation method is used to get the order of accuracy p ξ and the estimated value of error δ Reξ1 in Equation (31) and Equation (32).
where r ξ is refinement ratio of actuator point number.
where C ξ is correction factor and defined by Equation (34).
Following the calculation steps above, the uncertainty coefficient of C p and C t can be achieved. They are 0.129% D and 0.126% D, where D is the reference data from the NREL technical report [8]. Because both of them are less than 1% D, that proves the results are credible. Through the verification above, the level 2 actuator point number will be adopted, and it will obtain reliable results. where C ξ is correction factor and defined by Equation (34).
Following the calculation steps above, the uncertainty coefficient of Cp and Ct can be achieved. They are 0.129% D and 0.126% D, where D is the reference data from the NREL technical report [8]. Because both of them are less than 1% D, that proves the results are credible. Through the verification above, the level 2 actuator point number will be adopted, and it will obtain reliable results.  In addition, one of this method's advantages is that in the part of the aerodynamic solver it will cost less computational resources than CFD. The BEM theory is not considered here because the wake flow cannot be achieved in this way. The comparison of the computational cost between ALM and CFD method are shown in Table 4. It can be concluded that ALM uses less grids and computes faster than the CFD method under similar accuracy. In addition, one of this method's advantages is that in the part of the aerodynamic solver it will cost less computational resources than CFD. The BEM theory is not considered here because the wake flow cannot be achieved in this way. The comparison of the computational cost between ALM and CFD method are shown in Table 4. It can be concluded that ALM uses less grids and computes faster than the CFD method under similar accuracy.

Comparison of the Power and the Thrust
In this part, the results of aerodynamic performance are compared with 8 cases, in which 7 cases are the previous research and 1 case is the present one. Case 1 is the simulation of EALM and it is the present study. Case 2 is from the official technical report even though some data are given by NREL's FAST code [8]. In this report, the blade mode is derived by the mode's program and then passes a best-fit polynomial to get the equivalent polynomial representations of the mode shapes needed by FAST. Case 3 is simulated by HAWC2 (Horizontal Axis Wind turbine simulation Code 2nd generation) [29], which is an in-house nonlinear aeroelastic model developed by Technical University of Denmark. BEM is used as aerodynamic model and the multi-body dynamics method (MBD) is its structural model, in which each body is a linear Timoshenko beam element. The data of Case 4 are calculated by Li et al. [17] using CFD and MBD. This approach uses a dynamic overset CFD code for aerodynamics and MBD code for motion responses. The coupled way is done by exchanging the information between the fluid and structure solver in explicit form. The results of Case 5 are given by Jeong et al. [30] with BEM. In their research, the FSI (Fluid-Structure interaction) applies a strong coupling method which is using a first order implicit-explicit coupling scheme. Case 6 is set up in the research of Ponta et al. [31] by BEM and dynamic rotor deformation model, where the effects of rotor deformation are incorporates in the computation of aerodynamic loads. Yu et al. [32,33] combined CFD and CSD (Computational Structural Dynamics) in Case 7. The coupling methodology in this research is made by adopting the delta-airload loose-coupling technology. The last Case 8 is the results of the actuator line finite-element beam method (ALFBM) by Ma et al. [13]. The process of this research is three parts: the CFD solver given by OpenFOAM, the aerodynamic solver calculated by ALM, and the structural solver by finite-element beam method. The coupled steps are reading the velocity in the flow field from t-∆t and adding it to source term, which is computed by aerodynamic solver and structural solver. To clearly compare the varieties, the detailed information about the solved method on the aerodynamic performance and the structural responses in different cases are given in Table 5. In this research, BEM theory cannot describe the detailed flow field, and the CFD method may cost too much time to calculate. Therefore, ALM is chosen to resolve aerodynamic performance. Compared to ALFBM, this research concentrates on the variation of the aeroelastic characteristics in wake flow caused by upstream wind turbine.  Table 6 lists the results of 8 cases with velocities of 8 m/s and 11.4 m/s, and their mean output power and thrust differences are described in a histogram to visualize disparity of these cases, as shown in Figures 9 and 10. According to the results, the NREL's reported results are chosen as the reference data. It can be seen that all the results of power are higher than the reference data, except Case 6. All the results of thrust are lower than the FAST result. The big difference of thrust results are found from Case 5 to Case 8. The main reason for this is that the tip loss correction is not considered, and it has a great influence on aerodynamic performance. Case 4 has unsatisfied power result and there is no data in Case 3 under the rated situation. Besides, Case 1 has the most approximate result to the reference data among all studies. Its power difference is less than 50 kW and thrust difference is not more than 50 kN.   (a) (b) To further examine the accuracy of EALM under different wind speed conditions, power and thrust are computed and compared with Case 2, which are shown in Figure 11. Besides, their relative errors are given in Table 7. The errors are smaller when the wind speed is close to the rated speed, which is 11.4 m/s, than those when the wind speed is below 8 m/s. This phenomenon is caused by the coefficients of lift and drag C l and C d in Equation (1). These two parameters are referenced from airfoil data tables and these tables are achieved by physical experiment which is given by NREL official report [8]. In that report, the airfoil data is only obtained under nearly rated wind speed. Therefore, the lift and drag of a blade element can get more accurate if the wind speed closed to 11.4 m/s. On the contrary, the result of the force in the blade element will not be satisfactory if the wind speed is far away from 11.4 m/s. In case of application that the airfoil data corresponds to the wind speed, the result will be perfect, which will be improved in future research. From results of Picture 11, the predict results show good agreements with NREL's reference data. In addition, the differences between NREL and EALM are getting small with the inlet wind velocity increasing whether the output power or thrust. Especially in the rated situation, the relative errors of both power and thrust are less than 5%.
wind speed, the result will be perfect, which will be improved in future research. From results of Picture 11, the predict results show good agreements with NREL's reference data. In addition, the differences between NREL and EALM are getting small with the inlet wind velocity increasing whether the output power or thrust. Especially in the rated situation, the relative errors of both power and thrust are less than 5%.
According to the validation above, it could be concluded that EALM can predict the power and thrust accurately. This conclusion can be predictable because the aerodynamic calculations of EALM are based on blade element theory. In this theory, the lift and drag of blade element are achieved from two-dimensional airfoil data, which is obtained from physical experiments.    According to the validation above, it could be concluded that EALM can predict the power and thrust accurately. This conclusion can be predictable because the aerodynamic calculations of EALM are based on blade element theory. In this theory, the lift and drag of blade element are achieved from two-dimensional airfoil data, which is obtained from physical experiments.

Comparison of the Tip Displacement
In this section, the structural responses of blade tip are compared in 5 cases. Case 1 is the present research and obtained by EALM. Case 2 is from the official technical report and its data are given by NREL's FAST code [8]. Case 3 is the results of CFD and multi-body dynamics method by Li et al. [17]. Case 4 is set up by Yu et al. [32,33], in which CFD and CSD (Computational Structural Dynamics) are combined to calculate the coupled problem. The last Case 5 is using the actuator line finite-element beam method (ALFBM) by Ma et al. [13]. Detailed information in different cases can refer to Table 5.
Similarly, Table 8 lists the tip displacement results of 5 cases with 8 m/s and 11.4 m/s, and their differences are provided in the form of histograms, which are shown in Figures 12 and 13. Additionally, the NREL's reported results are chosen as the reference data. In these pictures, 0 direction means out of the rotor plane and 1 direction is in the rotor plane. The results of Case 3 and Case 5 are larger than those in Case 2, while the results of Case 4 are smaller. At the same time, the tip displacement of Case 1 is the nearest with the reference data in 0 direction and the farthest away from that reference data in 1 direction among 5 cases. It indicates that the tip displacement calculated by the beam solver, used in actual design, is similar to the modal approach, used to achieve structural dynamic characteristics, in 0 direction. Besides, the means used in this research could only predict approximate results in 1 direction. Even so, this beam solver can be still utilized to provide structural responses because of little deflection in edgewise. calculated by the beam solver, used in actual design, is similar to the modal approach, used to achieve structural dynamic characteristics, in 0 direction. Besides, the means used in this research could only predict approximate results in 1 direction. Even so, this beam solver can be still utilized to provide structural responses because of little deflection in edgewise. The tip displacement and its relative error in varied wind speeds are compared with Case 2 in Figure 14 and Table 9. The predict tendency of the tip displacement in 0 direction is almost the same as the reference data from Figure 13a, and its relative error is less than 5%. That indicates once again that the result of beam solver used in this study is approximate to that of the modal approach in 0 direction. Although the predicted value is not satisfied in high wind velocity in 1 direction, the relative differences are still not more than 15% except in the velocity of 10 m/s.   The tip displacement and its relative error in varied wind speeds are compared with Case 2 in Figure 14 and Table 9. The predict tendency of the tip displacement in 0 direction is almost the same as the reference data from Figure 13a, and its relative error is less than 5%. That indicates once again that the result of beam solver used in this study is approximate to that of the modal approach in 0 direction. Although the predicted value is not satisfied in high wind velocity in 1 direction, the relative differences are still not more than 15% except in the velocity of 10 m/s.  To represent the exactitude of blade tip deformation in different positions when considering the rotating motion, the azimuthal variations of tip displacement are compared with FAST code in Figures 15 and 16. In addition, the Pearson simplified correlation coefficient r is introduced to  To represent the exactitude of blade tip deformation in different positions when considering the rotating motion, the azimuthal variations of tip displacement are compared with FAST code in Figures 15 and 16. In addition, the Pearson simplified correlation coefficient r is introduced to describe the fit degree of two curves. It is a measure of the linear correlation between two variables in statistics. It has a value from −1 to +1, where 1 is total positive linear correlation, 0 is no linear correlation, and −1 is total negative linear correlation. Moreover, the larger the absolute value of the correlation coefficient is, the stronger proximity of two variables becomes. In this study, the coefficient r is rearranged in the formulas like Equation (35).
where n is sample size; x i and y i are the individual sample points indexed with i, and they represent the data from FAST and EALM in this research.
where n is sample size; i x and i y are the individual sample points indexed with i, and they represent the data from FAST and EALM in this research. From the Figure 15, the results agree well with each other. Moreover, it can be seen that all the Pearson correlation coefficients are not less than 99% from Table 10, which means the curves of two methods have strong relativity. Consequently, it could be concluded that EALM can calculate reliable structural responses.  where n is sample size; i x and i y are the individual sample points indexed with i, and they represent the data from FAST and EALM in this research. From the Figure 15, the results agree well with each other. Moreover, it can be seen that all the Pearson correlation coefficients are not less than 99% from Table 10, which means the curves of two methods have strong relativity. Consequently, it could be concluded that EALM can calculate reliable structural responses.  From the Figure 15, the results agree well with each other. Moreover, it can be seen that all the Pearson correlation coefficients are not less than 99% from Table 10, which means the curves of two methods have strong relativity. Consequently, it could be concluded that EALM can calculate reliable structural responses.

Influence on Aerodynamic Performance
This section will discuss the influence of elastic blade on the aerodynamic performance. Figures 17  and 18 depict the power and thrust azimuthal history with traditional actuator line model (ALM) and elastic actuator line model (EALM) at 8 m/s and 11.4 m/s. The ALM results are achieved by hiding the part of rotating beam solver in EALM, and the parameters setting are all the same. The curves fluctuate every 60 degrees, which is caused by the tower shadow effect. When the structural deformation is considered, the mean power reduces about 11.38 kW in 8 m/s and 11.43 kW in 11.4 m/s. Besides, the averaged thrust decreases approximately 1.08 kN in 8 m/s and 0.35 kN in 11.4 m/s. It indicates that the influence of elasticity on mean power and thrust is small in the rated condition. As the blade passes through tower, the reduction of output power turns into 15.33 kW in 8 m/s and 31.45 kW in 11.4 m/s. The decrement of thrust changes into 1.32 kN in 8 m/s and 0.96 kN in 11.4 m/s. The distance between tower and blade gets closer due to the blade deformation. Hence, the tower shadow effect becomes more serious and the decrement of power and thrust is larger. Because the blade displacement in high wind speed is large, the tower effect is the most serious in the rated wind speed and its reduction is more than that in 8 m/s.

Influence on Aerodynamic Performance
This section will discuss the influence of elastic blade on the aerodynamic performance. Figures  17 and 18 depict the power and thrust azimuthal history with traditional actuator line model (ALM) and elastic actuator line model (EALM) at 8 m/s and 11.4 m/s. The ALM results are achieved by hiding the part of rotating beam solver in EALM, and the parameters setting are all the same. The curves fluctuate every 60 degrees, which is caused by the tower shadow effect. When the structural deformation is considered, the mean power reduces about 11.38 kW in 8 m/s and 11.43 kW in 11.4 m/s. Besides, the averaged thrust decreases approximately 1.08 kN in 8 m/s and 0.35 kN in 11.4 m/s. It indicates that the influence of elasticity on mean power and thrust is small in the rated condition. As the blade passes through tower, the reduction of output power turns into 15.33 kW in 8 m/s and 31.45 kW in 11.4 m/s. The decrement of thrust changes into 1.32 kN in 8 m/s and 0.96 kN in 11.4 m/s. The distance between tower and blade gets closer due to the blade deformation. Hence, the tower shadow effect becomes more serious and the decrement of power and thrust is larger. Because the blade displacement in high wind speed is large, the tower effect is the most serious in the rated wind speed and its reduction is more than that in 8 m/s.  The stabilization of the normal and tangential force is present by standard deviation, as shown in Figure 19. The standard deviation s is defined as follows: where n is sample size; i x is the individual sample points indexed with I; and x is expectation, which is the mean value of sample. From the figure, the standard deviation of the loads on blade increases a little when the elasticity of the blades is considered. Figure 20 shows the wake structure comparison with EALM and ALM at 8 m/s. The wake vorticities are described by Q criterion, in which Q is calculated by the second invariant of the velocity gradient tensor and are stained by velocity field. The tip and hub vorticities are clear to be The stabilization of the normal and tangential force is present by standard deviation, as shown in Figure 19. The standard deviation s is defined as follows: where n is sample size; x i is the individual sample points indexed with I; and x is expectation, which is the mean value of sample. From the figure, the standard deviation of the loads on blade increases a little when the elasticity of the blades is considered. Figure 20 shows the wake structure comparison with EALM and ALM at 8 m/s. The wake vorticities are described by Q criterion, in which Q is calculated by the second invariant of the velocity gradient tensor and are stained by velocity field. The tip and hub vorticities are clear to be seen behind the rotor. The velocity in the internal surface of vorticities is smaller than the outside one. That is because the wind turbine extract momentum from the incoming airflow passing through it, and the wind energy is turned into mechanical energy, which causes the wind velocity after the rotor decreased.
The blade deflection is abbreviated and described in Figure 21. The deformation is too small to be discovered in 1 direction from left side view. Besides, it focuses on blade tip in 0 direction from front and vertical side view.
(a) (b) Figure 19. Stabilization of forces on blade: (a) normal force; and (b) tangential force. Figure 19. Stabilization of forces on blade: (a) normal force; and (b) tangential force. Figure 20 shows the wake structure comparison with EALM and ALM at 8 m/s. The wake vorticities are described by Q criterion, in which Q is calculated by the second invariant of the velocity gradient tensor and are stained by velocity field. The tip and hub vorticities are clear to be seen behind the rotor. The velocity in the internal surface of vorticities is smaller than the outside one. That is because the wind turbine extract momentum from the incoming airflow passing through it, and the wind energy is turned into mechanical energy, which causes the wind velocity after the rotor decreased. The blade deflection is abbreviated and described in Figure 21. The deformation is too small to be discovered in 1 direction from left side view. Besides, it focuses on blade tip in 0 direction from front and vertical side view.

Wake Flows Analysis
To analyze the influence of elastic blade in the wake flows, an array of two NREL 5 MW wind turbines is studied (see Figure 7) with EALM. The results are extracted after the upstream wind turbine rotating 40 revolutions, and the corresponding simulation time is around 260 s. According to research [28], the rotor speed of downstream wind turbine is set as 8.11 rpm to maximize its output power.
The comparison of mean output power among EALM, ALM, and Jha et al. [34] is presented in Table 11, in which the power generated by downstream wind turbine reduces a lot and only about 40% of upstream wind turbine because of the wake flow effects. In addition, the elastic blade makes the mean power of downstream wind turbine decrease about 14.6 kW, which is bigger than 11.38 kW of upstream wind turbine. Furthermore, the power and thrust azimuthal history of downstream wind turbine are depicted in Figure 22. When the blade passes through tower, the output power reduces 24.83 kW and the thrust decreases 7.01 kN. In contrast with the reduction of 8 m/s in Section 4.1, which are 15.33 kW and 1.32 kN, respectively, the tower shadow effect becomes more serious. This phenomenon is caused by the instability of the flow filed around the blade. The velocity and turbulence intensity in the wake of upstream wind turbine are changed into unstable. In the flow

Wake Flows Analysis
To analyze the influence of elastic blade in the wake flows, an array of two NREL 5 MW wind turbines is studied (see Figure 7) with EALM. The results are extracted after the upstream wind turbine rotating 40 revolutions, and the corresponding simulation time is around 260 s. According to research [28], the rotor speed of downstream wind turbine is set as 8.11 rpm to maximize its output power.
The comparison of mean output power among EALM, ALM, and Jha et al. [34] is presented in Table 11, in which the power generated by downstream wind turbine reduces a lot and only about 40% of upstream wind turbine because of the wake flow effects. In addition, the elastic blade makes the mean power of downstream wind turbine decrease about 14.6 kW, which is bigger than 11.38 kW of upstream wind turbine. Furthermore, the power and thrust azimuthal history of downstream wind turbine are depicted in Figure 22. When the blade passes through tower, the output power reduces 24.83 kW and the thrust decreases 7.01 kN. In contrast with the reduction of 8 m/s in Section 4.1, which are 15.33 kW and 1.32 kN, respectively, the tower shadow effect becomes more serious. This phenomenon is caused by the instability of the flow filed around the blade. The velocity and turbulence intensity in the wake of upstream wind turbine are changed into unstable. In the flow field of the downstream wind turbine blade, the variety of blade deflection with time will further disturb them. That will impact the relative velocity of the blade element. Thus, it can be concluded that the decrement of power and thrust caused by elastic blade are getting larger in the wake flows according to the comparisons above. The stabilization of the normal and tangential force on downstream wind turbine blade are compared again in Figure 23. In contrast with Figure 19, the standard deviation of the loads on blade in the wake caused by upstream wind turbine getting larger when the elasticity of the blades is considered. It indicates that the influence of elastic blade becomes more serious when it is in the wake of upstream wind turbine. That will aggravate the fatigue loads and should be paid much more attention during corresponding studies. The stabilization of the normal and tangential force on downstream wind turbine blade are compared again in Figure 23. In contrast with Figure 19, the standard deviation of the loads on blade in the wake caused by upstream wind turbine getting larger when the elasticity of the blades is considered. It indicates that the influence of elastic blade becomes more serious when it is in the wake of upstream wind turbine. That will aggravate the fatigue loads and should be paid much more attention during corresponding studies. The stabilization of the normal and tangential force on downstream wind turbine blade are compared again in Figure 23. In contrast with Figure 19, the standard deviation of the loads on blade in the wake caused by upstream wind turbine getting larger when the elasticity of the blades is considered. It indicates that the influence of elastic blade becomes more serious when it is in the wake of upstream wind turbine. That will aggravate the fatigue loads and should be paid much more attention during corresponding studies.  Figure 24 illustrates the process of double wind turbines array vorticity development. The wake expands as expected and the tip vorticity turns into continuum gradually in Figure 24a. Then the vorticity breakdown and the smaller-scale turbulence appears with distance increasing from upstream wind turbine (referring to Figure 24b). When the downstream wind turbine meets with this complicated wake, where the turbulence levels raise because of upstream wind turbine wakes and the momentum has not recovery completely, its tip vorticities break down earlier. That makes the wakes behind downstream wind turbine more complex, which can be seen in Figure 24c. Figure 25 shows the development of wake velocity field correspondingly. The phenomenon of velocity deficit is easy to observe behind the rotor. Before the downstream wind turbine disturbed from the wakes generated by upstream wind turbine, the velocity field is stable relatively. As the interference status appears, the distribution of velocity field in nearby regions of downstream wind turbine is in a muddle condition. That will increase the fatigue loads on blade of downstream wind turbine.   Figure 25 shows the development of wake velocity field correspondingly. The phenomenon of velocity deficit is easy to observe behind the rotor. Before the downstream wind turbine disturbed from the wakes generated by upstream wind turbine, the velocity field is stable relatively. As the interference status appears, the distribution of velocity field in nearby regions of downstream wind turbine is in a muddle condition. That will increase the fatigue loads on blade of downstream wind turbine.

Conclusions and Discussions
A new elastic actuator line model, which combines the traditional actuator line model and a beam solver used in the blade design, is proposed to study the aeroelastic performance of wind turbine efficiently. The validation is set up through comparing the result of the power, the thrust, and the tip deflection with different research. It is found that this new model can obtain acceptable prediction, including aeroelastic performance and wake fields. Besides, the tip displacement in 0 direction is the nearest with the reference data from official technical report given by NREL's FAST code among these studies. Furthermore, the elasticity of the blades can aggravate tower shadow effect and has the possibility of wind turbine instability. The fatigue loads on the blades also become more serious. When the wake effect generated by upstream wind turbine is considered, the impact above will get more obvious.
However, this new elastic actuator line model is incomplete and needs to be improved. For example, the modal shape and natural frequency results are not yet considered. The purpose of this

Conclusions and Discussions
A new elastic actuator line model, which combines the traditional actuator line model and a beam solver used in the blade design, is proposed to study the aeroelastic performance of wind turbine efficiently. The validation is set up through comparing the result of the power, the thrust, and the tip deflection with different research. It is found that this new model can obtain acceptable prediction, including aeroelastic performance and wake fields. Besides, the tip displacement in 0 direction is the nearest with the reference data from official technical report given by NREL's FAST code among these studies. Furthermore, the elasticity of the blades can aggravate tower shadow effect and has the possibility of wind turbine instability. The fatigue loads on the blades also become more serious. When the wake effect generated by upstream wind turbine is considered, the impact above will get more obvious.
However, this new elastic actuator line model is incomplete and needs to be improved. For example, the modal shape and natural frequency results are not yet considered. The purpose of this study is to propose an approach for investigating the aeroelastic performance rapidly. Further study is still required to refine the solver of structural dynamic characteristics and wave interactions [35].