Inﬂuence of the Flexible Tower on Aeroelastic Loads of the Wind Turbine

: The ﬁnite element discretization of a tower system based on the two-node Euler-Bernoulli beam is carried out by taking the cubic Hermite polynomial as the form function of the beam unit, calculating the structural characteristic matrix of the tower system, and establishing the wind turbine-nacelle-tower multi-degree-of-freedom ﬁnite element numerical model. The equation for calculating the aerodynamic load for any nacelle attitude angle is derived. The effect of the ﬂexible tower vibration feedback on the aerodynamic load of the wind turbine is studied. The results show that, when the stiffness of the tower is large, the effect of having tower vibration feedback or not on the aeroelastic load of the wind turbine is small. For the more ﬂexible tower system, wind-induced vibration time-varying feedback will cause larger aeroelastic load variations, especially the top of the tower overturning moment, thus causing a larger impact on the dynamic behavior of the tower downwind and crosswind. As the ﬂexibility of the tower system increases, the interaction between tower vibration and pneumatic load is also gradually increasing. Taking into account the inﬂuence of ﬂexible towers on the aeroelastic load of a wind turbine can help predict the pneumatic load of a wind turbine more accurately and improve the efﬁciency of wind energy utilization on the one hand and analyze the dynamic behavior of the ﬂexible structure of a wind turbine more accurately on the other hand, which is extremely beneﬁcial to the structural optimization of wind turbine.


Introduction
As the supporting structure of the wind turbine, the safety and reliability of the tower are directly related to the operation state of the wind turbine [1]. During the operation of the wind turbine, the dynamic behavior of the blades and the tower influence each other. To obtain more accurate dynamic response results, the blade-tower coupling effect must be considered in the dynamic analysis of wind turbines [2,3]. For wind turbines with complex structures [4,5], to effectively analyze their dynamic behavior, it is necessary to establish a simplified model of their dynamics [6,7]. The accurate calculation of tower wind load has an important influence on the accurate assessment of tower dynamic response. The nacelle is one of the main components of a wind turbine, and different nacelle attitude angles correspond to different aerodynamic loads of the wind turbine.
In this research, a simplified multi-degree-of-freedom numerical model of the horizontal axis wind turbine tower system, including the wind turbine-nacelle-tower foundation [8], was established and used to study the dynamic behavior of offshore HAWT tower system [9]. At the same time, based on the blade element momentum theory, a calculation method for the aerodynamic load of rotating blades including time-varying nacelle attitude feedback is proposed. This method uses a set of Euler angles to describe the time-varying attitude of the horizontal-axis wind turbine nacelle. This paper takes four basic cases as the research object and calculates the load on the top of the tower under the action of the flexibility of the foundation. These loads are used to analyze the dynamic response of the Appl. Sci. 2021, 11, 8876 2 of 19 tower system and discuss the influence of the flexibility characteristics of the foundation on the dynamic characteristics of the tower.

Equations of Motion for Vibrating Systems
To effectively analyze the inherent characteristics of the wind turbine tower system, a simplified model of the tower system was developed, as shown in Figure 1. The model simplifies the wind turbine nacelle and the wind wheel as concentrated masses M s , which are attached to the top of the tower. Additionally, the flexibility of the HAWT foundation is simulated using rotating springs and horizontal tension springs.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 2 of 19 the time-varying attitude of the horizontal-axis wind turbine nacelle. This paper takes four basic cases as the research object and calculates the load on the top of the tower under the action of the flexibility of the foundation. These loads are used to analyze the dynamic response of the tower system and discuss the influence of the flexibility characteristics of the foundation on the dynamic characteristics of the tower.

Equations of Motion for Vibrating Systems
To effectively analyze the inherent characteristics of the wind turbine tower system, a simplified model of the tower system was developed, as shown in Figure 1. The model simplifies the wind turbine nacelle and the wind wheel as concentrated masses , which are attached to the top of the tower. Additionally, the flexibility of the HAWT foundation is simulated using rotating springs and horizontal tension springs. Based on the geometric properties of the structure, material properties, and external excitation properties, the multi-degree-of-freedom equation of motion is generally obtained as: (1) where M is the mass matrix; C is the structural damping matrix, and K is the stiffness matrix, which is the external load acting on the tower.

Tower Finite Element Modeling
In this paper, the top mass of the tower is simplified to a concentrated mass block and attached to the top of the tower. As for the thin-walled tower with a variable cross-section, a two-node Euler-Bernoulli beam is used to discretize the finite element of the wind turbine tower, and a cubic Hermite polynomial is chosen as the form function of each beam unit. Figure 2 shows a beam unit of the wind turbine tower with length l, mass per unit length m(x), and bending stiffness EI(x). The two nodes of this unit are located at the two endpoints of the beam; through these two nodes, the finite unit can be assembled into a single structure. Based on the geometric properties of the structure, material properties, and external excitation properties, the multi-degree-of-freedom equation of motion is generally obtained as: M ..
where M is the mass matrix; C is the structural damping matrix, and K is the stiffness matrix, which is the external load acting on the tower.

Tower Finite Element Modeling
In this paper, the top mass of the tower is simplified to a concentrated mass block and attached to the top of the tower. As for the thin-walled tower with a variable crosssection, a two-node Euler-Bernoulli beam is used to discretize the finite element of the wind turbine tower, and a cubic Hermite polynomial is chosen as the form function of each beam unit. Figure 2 shows a beam unit of the wind turbine tower with length l, mass per unit length m(x), and bending stiffness EI(x). The two nodes of this unit are located at the two endpoints of the beam; through these two nodes, the finite unit can be assembled into a single structure.
If only plane displacement is considered, each node of the beam unit has only two degrees of freedom, which are lateral displacement and rotation. The relationship between the displacement of the beam unit and the four degrees of freedom can be expressed as: where ψ i (x) is defined as the unit displacement due to the occurrence of unit displacement u i (x) while keeping the other degrees of freedom at zero. ψ i (x) satisfies the following boundary conditions:  If only plane displacement is considered, each node of the beam unit has only two degrees of freedom, which are lateral displacement and rotation. The relationship between the displacement of the beam unit and the four degrees of freedom can be expressed as: , where is defined as the unit displacement due to the occurrence of unit displacement while keeping the other degrees of freedom at zero. satisfies the following boundary conditions: The equilibrium equation for a beam of equal cross-section considering only the load acting at both ends without taking into account the shear deformation is: Its general solution is: . The constant can be determined by substituting the boundary conditions, which gives: The equilibrium equation for a beam of equal cross-section considering only the load acting at both ends without taking into account the shear deformation is: Its general solution is: u(x) = a 1 + a 2 x l + a 3 x l 2 + a 4 x l 3 . The constant a i can be determined by substituting the boundary conditions, which gives: where ψ 1 (x), ψ 2 (x), ψ 3 (x), and ψ 4 (x) are the cubic Hermite interpolation functions, which ensure the continuity of deflection and cornering between beam units at the boundary. Using the principle of imaginary displacement for a beam unit with bending stiffness EI(x) and length l, the general expressions for the stiffness influence coefficient k ij and the mass influence coefficient m ij of the beam unit can be derived as follows: For cells with EI(x) = EI and uniformly distributed masses (m(x) ≡ m), when i, j = 1, 2, 3, and 4, the analytical solutions of the cell stiffness matrix and the consistent mass matrix can be calculated as: The wind turbine tower is always subjected to the gravitational force of the top mass of the tower, and this axial pressure also has a certain degree of influence on the stiffness of the tower. Using Hermite cubic polynomials as an interpolation function yields the general form of the unit geometric stiffness influence coefficient: where N(x) denotes the axial force acting on the unit. When N(x) = N G is constant, the consistent geometric stiffness matrix of the unit is obtained as: The axial force along the full length of the unit is assumed to be constant. By superimposing the above two effects, elastic and geometric, the total stiffness matrix of the structural unit is obtained, expressed as: where the negative sign indicates that the presence of axial load N(x) increases the deflection of the tower. If units m, n, and p are connected to node i of the structure, the stiffness coefficient of this node can be calculated using the direct stiffness method as: where the notation 'ˆ' indicates the coefficient in the overall coordinate system. The stiffness matrix of the whole structure can be obtained by Equation (13). Similarly, the consistent mass matrix of the whole structure can be calculated.
The damping characteristics of structures generally need to be determined by experimental methods, but it would be quite difficult and impractical to obtain the damping of each structure by experimental tests. Therefore, the dynamic properties of a structure are often analyzed by using a classical damping model. In this paper, the Rayleigh damping model is used, which integrates the effects of mass-proportional damping and stiffness-proportional damping, and its calculation equation is as follows: where η 0 and η 1 are scaling factors.
where ω m and ω n are the m-order and n-order frequencies of the vibrating system, and ξ m and ξ n are the damping ratios of the corresponding modal frequencies [8].
When ξ m = ξ n = ξ, it can be simplified as follows: In practical engineering applications, ω m is usually taken as the fundamental frequency of the vibrating system, and ω n is taken as the higher order frequency that has a significant effect on the dynamic behavior of the system. When the modal damping and modal frequency of the tower system are known, its damping matrix can be determined by Equations (15) and (16). After calculating the mass matrix, stiffness matrix and damping matrix of the whole tower, the multi-degree-of-freedom finite element numerical model of the tower system can be established by treating the top mass of the tower as a boundary condition [10].

Eigenvalue Analysis
The eigenvalue problem of the structure is to analyze its free vibration. To calculate the natural frequency and modal vibration pattern of the tower system, the effect of damping is generally neglected. Removing the damping matrix and external excitation of Equation (1), the free vibration equation of the tower can be obtained as follows: The characteristic equation of the tower can be solved as: where ω denotes the frequency of free vibration, and φ denotes the time-independent nth-order vector. φ = φ 1 φ 2 · · · φ n , φ n denote the nth-order vibration column vector of the structure. Wind turbines operating in different seas have different foundation characteristics. Adhikari et al. [11] studied and obtained classical values of foundation stiffness for three actual operating wind turbines. Based on this study, 12 foundation cases are selected in this paper to study the effect of flexible foundations on the towers' inherent frequency. The wind turbine numerical model parameters used for simulation in this paper were selected from the NREL 5 MW offshore horizontal axis wind turbine [12]. The first three orders of tower inherent frequencies are calculated as shown in Table 1.
From Table 1, it can be concluded that the rotational stiffness has a greater effect on the first two orders of the towers' natural frequency and a smaller effect on the third order of the towers' natural frequency. The horizontal stiffness of a tower on the natural frequency of the tower is opposite to the effect of rotational stiffness. For a tower with a rigid foundation, the first three orders of inherent frequency are also calculated, as shown in Table 2.
From Tables 1 and 2, it can be concluded that the flexible foundation of a tower has an important effect on the free vibration of the tower. Therefore, the flexibility of the tower needs to be taken into account in the dynamic analysis of offshore wind turbines.

Solving the Nacelle Attitude Angle Using Euler Angles
Attitude solutions are relatively common in computer graphics and aerospace, and Euler angles [13] have been widely used in the field of attitude description due to their simplicity and ease of use. A set of Euler angles, which are the yaw, pitch, and roll angles obtained from the rotation of the nacelle about the X, Y, and Z coordinate axes, are used to describe the attitude of the HAWT nacelle, as shown in Figure 3. From Table 1, it can be concluded that the rotational stiffness has a greater effect on the first two orders of the towersʹ natural frequency and a smaller effect on the third order of the towersʹ natural frequency. The horizontal stiffness of a tower on the natural frequency of the tower is opposite to the effect of rotational stiffness. For a tower with a rigid foundation, the first three orders of inherent frequency are also calculated, as shown in Table 2. From Tables 1 and 2, it can be concluded that the flexible foundation of a tower has an important effect on the free vibration of the tower. Therefore, the flexibility of the tower needs to be taken into account in the dynamic analysis of offshore wind turbines.

Solving the Nacelle Attitude Angle Using Euler Angles
Attitude solutions are relatively common in computer graphics and aerospace, and Euler angles [13] have been widely used in the field of attitude description due to their simplicity and ease of use. A set of Euler angles, which are the yaw, pitch, and roll angles obtained from the rotation of the nacelle about the X, Y, and Z coordinate axes, are used to describe the attitude of the HAWT nacelle, as shown in Figure 3. Assume that the nacelle transitions from one attitude to another in the sequence of transverse roll-pitch-yaw, wherein the transition matrix around the axis is: Assume that the nacelle transitions from one attitude to another in the sequence of transverse roll-pitch-yaw, wherein the transition matrix around the axis is: The attitude matrix of the nacelle can be obtained as follows:

Calculation of Aerodynamic Loads under Consideration of NAF
For any nacelle attitude angle, the horizontal axis wind turbine coordinate system shown in Figure 4 is established to facilitate the calculation of the spatial position and wind speed components of each wing section of the blade [14].

Calculation of Aerodynamic Loads Under Consideration of NAF
For any nacelle attitude angle, the horizontal axis wind turbine coordinate system shown in Figure 4 is established to facilitate the calculation of the spatial position and wind speed components of each wing section of the blade [14].
1. The inertial coordinate system S coincides with the tower bottom coordinate system , and the nacelle coordinate system coincides with the tower top coordinate system . 2. When the azimuth angle θ = 0, the wind turbine rotation plane coordinate system coincides with the hub coordinate system , and rotates with the azimuth angle θ. 3. The tower height is L, and the horizontal distance from the top of the tower to the hub is . 4. The spindle inclination angle and blade taper angle are and . The transformation relationship between the blade coordinate system and the inertial coordinate system exists both in terms of rotation and translation. The transformation between the coordinate systems is described using the chi-square coordinates [15] as 1 . The transformation relationship between the coordinate systems is as follows:

1.
The inertial coordinate system S coincides with the tower bottom coordinate system S BOT , and the nacelle coordinate system S N coincides with the tower top coordinate system S top .

2.
When the azimuth angle θ = 0, the wind turbine rotation plane coordinate system S R coincides with the hub coordinate system L H , and S R rotates with the azimuth angle θ.

3.
The tower height is L, and the horizontal distance from the top of the tower to the hub is L H . 4.
The spindle inclination angle and blade taper angle are β T and β C .
The transformation relationship between the blade coordinate system and the inertial coordinate system exists both in terms of rotation and translation. The transformation between the coordinate systems is described using the chi-square coordinates [15] as S = [ X Y Z 1 ] T . The transformation relationship between the coordinate systems is as follows: where: K L is the transformation matrix from the inertial coordinate system to the nacelle coordinate system; K Lh and sub K βT are the translation and rotation transformation matrices from the nacelle to the hub coordinate system, respectively; R N and K DN are the attitude matrix and deflection matrix of the nacelle, respectively; K θ is the transformation matrix from the hub to the rotating coordinate system of the wind wheel, and K βC is the transformation matrix from the wind wheel to the rotating blade.
From Equations (23) to (32), the transformation relation between the blade coordinate system and the inertial coordinate system can be obtained as: Therefore, the coordinate component X Y Z T in the inertial coordinate system for the leaf element coordinates X B Y B Z B T = 0 0 r T is calculated by Equation (34).
The coordinate component of the wind speed U S = U 0 0 T in the inertial coordinate system in the blade coordinate system can be calculated by the following equation.  The load on the wind turbine tower is mainly from the inertial load on the top mass of the tower and the aerodynamic load on the rotating blades. The blade load and tower The load on the wind turbine tower is mainly from the inertial load on the top mass of the tower and the aerodynamic load on the rotating blades. The blade load and tower top load can be decomposed into component forces and bending moments along the three coordinate axes, as shown in Figure 5. When the wind component of turbulent wind in the lobe element coordinate system is obtained, the angle of entry and relative wind speed of the incoming wind on the lobe element can be calculated by the following equation.
In this paper, the dynamic response of the tower in the forward-backward direction and the left-right direction is investigated. The aerodynamic forces and aerodynamic moments on the blades were calculated using BEM theory [16] with the following equations.
where and denote the lift and drag coefficients of the cross-sectional airfoil; ρ is the air density, and c is the chord length of the blade cross-sectional airfoil.

Tower Load Calculation
The aerodynamic forces F and aerodynamic moments M on the wind turbine are transformed by coordinates to obtain the aerodynamic loads acting on the tower. When the wind component of turbulent wind in the lobe element coordinate system is obtained, the angle of entry and relative wind speed of the incoming wind on the lobe element can be calculated by the following equation.
In this paper, the dynamic response of the tower in the forward-backward direction and the left-right direction is investigated. The aerodynamic forces and aerodynamic moments on the blades were calculated using BEM theory [16] with the following equations.
where C L and C D denote the lift and drag coefficients of the cross-sectional airfoil; ρ is the air density, and c is the chord length of the blade cross-sectional airfoil.

Tower Load Calculation
The aerodynamic forces F and aerodynamic moments M on the wind turbine are transformed by coordinates to obtain the aerodynamic loads acting on the tower. where The eccentricity of the top mass center of the wind turbine tower with the tower center will produce the overturning moment and pitching moment acting on the tower. In the tower top coordinate system, assuming the coordinates of the tower top mass center as X TE Y TE Z TE , the tower top bending moment caused by the eccentricity of the tower top mass can be calculated by the following equation.
where M s is the tower top mass, and g is the acceleration of gravity. The final total load of the tower can be obtained by combining the wind turbine pneumatic load and the eccentric moment of the top mass of the tower on the tower.

Solving Structural Equations of Motion
For a general external excitation P(t), the Newmark method is used in this paper to solve the multi-degree-of-freedom equations of motion.
The Newmark method approximates the velocity and displacement of the system at the moment (t + ∆t) by two assumptions. . x t+∆t = .
x t+∆t When β ≥ 0.5 and α ≥ 0.25(0.5 + β) 2 , the Newmark method is the unconditionally stable format. The equations of motion can be solved by the following steps.

1.
Firstly, the overall characteristic matrices K, M and C of the vibrating system are calculated.

2.
From the initial conditions x 0 and .
Step ∆t is selected with parameters α, β and calculate the integration constants.
The equivalent stiffness matrix K is calculated. 2.
For each time step 1. The equivalent load vector F at moment t + ∆t is determined. .. ..
x t C (52) 2. The displacement at moment t + ∆t is calculated.
3. The velocity at time t + ∆t and the acceleration are solved. .. .

Dynamic Response Analysis Process
The dynamic response of the tower system is calculate as shown in Figure 6, taking into account the effects of nacelle attitude feedback and foundation flexibility of the foundation.
3. The velocity at time Δ and the acceleration are solved.

Dynamic Response Analysis Process
The dynamic response of the tower system is calculate as shown in Figure 6, taking into account the effects of nacelle attitude feedback and foundation flexibility of the foundation.

Tower Load Error Analysis
Under the action of wind load, the tower will vibrate, which will cause the pitch and roll angle of the nacelle to change, and the change of attitude angle will in turn affect the load calculation of the rotating blade, which is a cyclic process. In this paper, we design a method to analyze the dynamics under the action of dynamic structural parameter feedback, and use MATLAB to prepare the corresponding program to explore the laws through numerical examples [9]. In analyzing the effect of flexible pylons on aeroelastic loads, the following assumptions are made first.

1.
The nacelle and the tower have only yaw motion.

2.
The center of mass of the nacelle coincides with the center of the tower top.
The deflection of the tower top in downwind direction and crosswind direction will be equal to the fore-and-aft displacement and left-right displacement of the nacelle respectively, and the angle of rotation of the tower top in fore-and-aft direction and leftright direction will be equal to the pitch angle and cross-roll angle of the nacelle, respectively. The relationship between the tower top turning angle and the nacelle attitude angle is shown in Figure 7.

Tower Load Error Analysis
Under the action of wind load, the tower will vibrate, which will cause the pitch and roll angle of the nacelle to change, and the change of attitude angle will in turn affect the load calculation of the rotating blade, which is a cyclic process. In this paper, we design a method to analyze the dynamics under the action of dynamic structural parameter feedback, and use MATLAB to prepare the corresponding program to explore the laws through numerical examples [9]. In analyzing the effect of flexible pylons on aeroelastic loads, the following assumptions are made first. The deflection of the tower top in downwind direction and crosswind direction will be equal to the fore-and-aft displacement and left-right displacement of the nacelle respectively, and the angle of rotation of the tower top in fore-and-aft direction and leftright direction will be equal to the pitch angle and cross-roll angle of the nacelle, respectively. The relationship between the tower top turning angle and the nacelle attitude angle is shown in Figure 7. and denote the angle of rotation of the top of the tower in the front-to-back direction and the left-to-right direction, respectively. Based on the assumptions of this paper, there are , holds at the same time. The load calculation formula for any nacelle attitude angle derived in this paper is applied to the study of the effect of the flexible tower on the aeroelastic load of the wind turbine. In the steady-state yaw condition, when considering the effect of the flexible tower vibration feedback on the wind turbine aeroelastic load, it is necessary to calculate the large feedback effect of the tower top front-to-back and left-right directional turning angle and displacement at the same time, and it is also necessary to consider the two-way feedback effect of the nacelle's cross-roll angle, pitch angle, downwind tower top displacement, and cross-wind tower top displacement at the same time, and its calculation process is shown in Figure 6. At the same time, the calculation flow without considering the feedback effect is also incorporated into this procedure, only when the feedback is not considered, the influence of tower vibration on the structural parameters needs to be compensated. ψ ty and ψ tx denote the angle of rotation of the top of the tower in the front-to-back direction and the left-to-right direction, respectively. Based on the assumptions of this paper, there are ψ y = ψ ty , ψ x = ψ tx holds at the same time. The load calculation formula for any nacelle attitude angle derived in this paper is applied to the study of the effect of the flexible tower on the aeroelastic load of the wind turbine. In the steady-state yaw condition, when considering the effect of the flexible tower vibration feedback on the wind turbine aeroelastic load, it is necessary to calculate the large feedback effect of the tower top front-to-back and left-right directional turning angle and displacement at the same time, and it is also necessary to consider the two-way feedback effect of the nacelle's cross-roll angle, pitch angle, downwind tower top displacement, and cross-wind tower top displacement at the same time, and its calculation process is shown in Figure 6. At the same time, the calculation flow without considering the feedback effect is also incorporated into this procedure, only when the feedback is not considered, the influence of tower vibration on the structural parameters needs to be compensated.
To simulate the roll of tower vibration feedback in the dynamic response of flexible towers, simulation analysis is performed for towers with different flexibility. First, the other parameters of the towers are maintained consistent with the NREL 5 MW wind turbine; the geometric parameters of the wind turbine tower and the basic operating parameters of the wind turbine are shown in Table 3. The increase of its flexibility is simulated by increasing the height of the tower, and the parameters of the tower used for numerical simulation are shown in Table 4. The relative errors of the tower top load are calculated separately for the two cases of considering feedback and not considering feedback; the comparative analysis results of longitudinal thrust and transverse thrust are shown in Figure 8.
As can be seen from Figure 8, when the tower stiffness is relatively large, the effect on the wind turbine pneumatic load is small when considering the tower bending-bending coupling vibration feedback. When the tower flexibility increases, its load relative error increases rapidly, especially the transverse thrust; for each flexible tower, its relative error amplitude is almost maintained at about two times the relative error of the longitudinal thrust, which will have a significant increase in impact for the more flexible tower, and cause a large impact on the power of the wind turbine.
To investigate the law of mutual influence between tower vibration and pneumatic load, we compared the relative errors of tower top angle and displacement under three tower heights, and the simulation results are shown in Figures 9 and 10.
From Figures 9 and 10, it can be seen that, when the flexibility of the support structure is large, the maximum errors of displacement and the turning angle of the tower top downwind are both around 8%, while the maximum deviations of the crosswind direction are both over 12%. It can be seen that, when the flexibility of the support structure is large, ignoring the effect of the time-varying nacelle attitude feedback will make the analysis of the dynamic less accurate. It can also be seen that the relative errors of displacement and the turning angle at the top of the tower remain the same in terms of trend and magnitude. It can be seen that when the other structural parameters (blade, nacelle mass, etc.) are constant, the degree of influence of the tower bending-bending coupling bidirectional feedback on the dynamical behavior of the tower is positively correlated with the flexibility of the support structure. The relative errors of the tower top load are calculated separately for the two cases of considering feedback and not considering feedback; the comparative analysis results of longitudinal thrust and transverse thrust are shown in Figure 8. As can be seen from Figure 8, when the tower stiffness is relatively large, the effect on the wind turbine pneumatic load is small when considering the tower bending-bending coupling vibration feedback. When the tower flexibility increases, its load relative error increases rapidly, especially the transverse thrust; for each flexible tower, its relative error amplitude is almost maintained at about two times the relative error of the longitudinal thrust, which will have a significant increase in impact for the more flexible tower, and cause a large impact on the power of the wind turbine To investigate the law of mutual influence between tower vibration and pneumatic load, we compared the relative errors of tower top angle and displacement under three tower heights, and the simulation results are shown in Figures 9 and 10. As can be seen from Figure 8, when the tower stiffness is relatively large, the effect on the wind turbine pneumatic load is small when considering the tower bending-bending coupling vibration feedback. When the tower flexibility increases, its load relative error increases rapidly, especially the transverse thrust; for each flexible tower, its relative error amplitude is almost maintained at about two times the relative error of the longitudinal thrust, which will have a significant increase in impact for the more flexible tower, and cause a large impact on the power of the wind turbine To investigate the law of mutual influence between tower vibration and pneumatic load, we compared the relative errors of tower top angle and displacement under three tower heights, and the simulation results are shown in Figures 9 and 10. (a) From Figures 9 and 10, it can be seen that, when the flexibility of the support structure is large, the maximum errors of displacement and the turning angle of the tower top From Figures 9 and 10, it can be seen that, when the flexibility of the support structure is large, the maximum errors of displacement and the turning angle of the tower top

Dynamic Response Analysis of Tower System
Four foundation cases were selected to study the effect of flexible foundations on tower loads. The wind condition used is a steady state wind with 11.4 m/s and a turbulence intensity of 16.4946%.
For each of the four foundation cases in Table 5, the wind load response of the tower system is solved, and the effect of the foundation rotational stiffness on the tower vibration is explored. When considering the flexibility, the wind turbine tower top response time equations are shown in Figures 11 and 12. Rigid foundation placement and the turning angle at the top of the tower remain the same in terms of trend and magnitude. It can be seen that when the other structural parameters (blade, nacelle mass, etc.) are constant, the degree of influence of the tower bending-bending coupling bidirectional feedback on the dynamical behavior of the tower is positively correlated with the flexibility of the support structure.

Dynamic Response Analysis of Tower System
Four foundation cases were selected to study the effect of flexible foundations on tower loads. The wind condition used is a steady state wind with 11.4 m/s and a turbulence intensity of 16.4946%.
For each of the four foundation cases in Table 5, the wind load response of the tower system is solved, and the effect of the foundation rotational stiffness on the tower vibration is explored. When considering the flexibility, the wind turbine tower top response time equations are shown in Figures 11 and 12.    As shown in Figure 11, the maximum displacements of the tower top in the downwind direction are 1.2 m, 0.82 m, 0.61 m, and 0.46 m. The maximum lateral displacements of the tower top are 0.14 m, 0.1 m, 0.08 m, and 0.07 m. It can be seen that the deflection of the tower top increases greatly when considering flexibility. Compared with the rigid foundation, the maximum tower top displacement of Case 1, Case 2, and Case 3 increases nearly 161%, 78%, and 32% in the downwind direction and nearly 100%, 42%, and 14% in the sidewind direction, respectively.
As shown in Figure 12, the tower top velocity response is greatly increased when considering flexibility. Compared with the rigid foundation, the tower top velocities of Case 1, Case 2, and Case 3 increase by nearly 88%, 41%, and 14% in the downwind direction, and by nearly 55%, 23%, and 4% in the sidewind direction, respectively.
From Figures 11 and 12, it can be seen that the tower top displacement and velocity response amplitude increase with decreasing base stiffness when considering base stiffness. Moreover, when the rotational stiffness gradually increases, the closer to the rigid foundation, the smaller the decrease of the top displacement and velocity response amplitude.

Conclusions
(1) The finite element discretization of the pylon system is carried out with a two-node Euler-Bernoulli beam; the cubic Hermite polynomial is taken as the form function of the beam unit; the structural characteristic matrix of the pylon system is calculated, and its multi-degree-of-freedom finite element model is established. (2) The calculation formula of the effect of the nacelle attitude angle on the aerodynamic load of the blade is obtained.
(3) When the tower stiffness is large, the effect of having or not having tower vibration feedback on the calculation of the wind turbine aeroelastic load is small. (4) For a more flexible tower system, wind-induced vibration time-varying feedback will cause a larger aeroelastic load variation, especially the lateral tilting moment at the top of the tower, thus causing a larger impact on the dynamical behavior of the tower downwind and crosswind. It can be seen that, as the flexibility of the tower system increases, the interaction between the tower vibration and the pneumatic load also increases gradually. Taking the influence of the flexible tower on the aeroelastic load of the wind turbine into account can help predict the wind turbine pneumatic load more accurately and improve the efficiency of wind energy utilization; on the other hand, it can analyze the dynamic behavior of the flexible structure of the wind turbine more accurately, which is extremely beneficial to the structural optimization design of the wind turbine. On the other hand, it can analyze the dynamic behavior of the flexible structure of wind turbines more accurately, which is extremely beneficial to the structural optimization design of wind turbines.

Conflicts of Interest:
The authors declare no conflict of interest.