Evaluation of an Uncoupled Method for Analyzing the Seismic Response of Wind Turbines Excited by Wind and Earthquake Loads

: There is a signiﬁcant interaction between wind and earthquakes for large-scaled wind turbines due to an aeroelastic e ﬀ ect. This study evaluates the accuracy of an uncoupled method extensively utilized to analyze the seismic response of wind turbines at the operational state. Initially, the oscillation of the blade for the National Renewable Energy Laboratory (NREL) 5 MW wind turbine excited by wind and wind-earthquake combination, respectively, is compared using the fully coupled method to verify the assumption in this uncoupled method. Subsequently, the inﬂuence of ground motions on the aerodynamic loadings of the rotor is discussed to evaluate the interaction between wind and earthquake loads. In addition, the accuracy of the uncoupled method is assessed by comparing the analysis results of the coupled and uncoupled methods, where di ﬀ erent mean wind speed and equivalent aerodynamic damping ratio are considered. The results indicate that the oscillation velocity of blades and thrust on the rotor are signiﬁcantly inﬂuenced by ground motions. Moreover, the amplitude of thrust variations caused by earthquakes increases monotonously with the oscillation velocity amplitude of blade-root. The errors between the two models are beyond the engineering margins for some earthquakes, such that it is di ﬃ cult to optimize a consistent aerodynamic damping in the uncoupled model to accurately predict the seismic response of wind turbines.


Introduction
As the primary form of wind energy utilization, wind power can fulfill the global energy demand and reduce CO 2 emissions [1]. The global annual installations of wind turbines have been more than 50 GW since 2014, and it will be nearly 71 GW until 2024 [2]. With more and more wind turbines being built in seismically active regions [3], it is necessary to accurately predict the seismic response of wind turbines for assessment and design purposes. The large-scaled wind turbine is a high-slender structure and sensitive to the loadings in the lateral directions. Consequently, the updated specifications [4,5] recommended that earthquake forces should be combined with aerodynamic loadings induced by the normal wind condition, where the analysis methods can be divided into coupled and uncoupled.
The coupled analysis of the seismic response of wind turbines was firstly performed by the GH BLADED software [6]. The aerodynamic loadings were calculated by blade element momentum (BEM) theory, and the simulation results indicated that the aeroelasticity of the wind turbine and influence of the controller could be modeled with this code. Subsequently, the seismic module was added to the FAST (Fatigue, Aerodynamics, Structures, and Turbulence) code, an open-source software, such that where f is the frequency (Hz); Vhub is the average wind speed of hub-height; and are the 122 standard deviation and integral scale parameter for every velocity component, respectively, in which 123 the subscript k refers to the component along the x-, y-, and z-axes. According to IEC 61400-1 [4], the 124 power law is adopted to determine the wind profile, where the exponent is 0.2. The IEC turbulence 125 level B is selected to determine all the parameters.

126
The cross power spectra density function between nodes i and j is defined as: where Si,i and Sj,j denote the wind velocity auto spectra at nodes i and j, respectively; (∆ , )is the 128 coherence function defined in Equation (3).
where Δr is the distance of the two nodes, b represents the coherence decrement, with a value of 12 130 for this study; LC represents the coherence scale parameter, with a value of 340.2 m according to IEC 131 61400-1 [4].
132 Figure 2. Element distribution of the blade.

Wind Field
The fluctuating wind velocity was assumed to be a stationary, random vector field. The Kaimal spectrum is used to describe the turbulent characteristic of wind process: where f is the frequency (Hz); V hub is the average wind speed of hub-height; σ k and L k are the standard deviation and integral scale parameter for every velocity component, respectively, in which the subscript k refers to the component along the x-, y-, and z-axes. According to IEC 61400-1 [4], the power law is adopted to determine the wind profile, where the exponent is 0.2. The IEC turbulence level B is selected to determine all the parameters. The cross power spectra density function between nodes i and j is defined as: Energies 2020, 13, 3833 5 of 27 where S i,i and S j,j denote the wind velocity auto spectra at nodes i and j, respectively; C(∆r, f ) is the coherence function defined in Equation (3).
where ∆r is the distance of the two nodes, b represents the coherence decrement, with a value of 12 for this study; L C represents the coherence scale parameter, with a value of 340.2 m according to IEC 61400-1 [4]. The TurbSim code [40] is utilized to generate the full-field wind samples with random seeds. The duration of the wind speed time-history is 650 s to cover the simulations. As shown in Figure 1, the wind field is discretized into a 160 m × 160 m vertical rectangle grid to encompass the entire rotor. The distance between adjacent nodes is 8 m in both the horizontal and vertical directions.

Earthquake Load
The earthquake has three components, and their relationship is stochastic. Therefore, some researchers selected one component to analyze the seismic response of wind turbines [12,41]. The difference between the coupled and uncoupled methods lies in the solving method of the aeroelastic effect. Under the excitation of wind, the aerodynamic loading in the FA direction is more significant than that of the side-to-side (SS) direction [11], such that the aeroelasticity is mainly in the FA direction of wind turbines. To eliminate the influence of earthquake forces in the SS direction, one of the two horizontal components of earthquakes which will cause larger a tower-base bending moment is selected as the input ground motion and applied along the FA direction of the wind turbine. The studies conducted by Santangelo et al. [31][32][33] indicated that the error between the coupled and uncoupled methods was related to earthquakes. Consequently, a set of earthquakes shown in Table 3 are selected from the database of Pacific Earthquake Engineering Research center (PEER) to evaluate the uncoupled model thoroughly, where the far field and near field records recommended by Federal Emergency Management Agency (FEMA) [42] are included. For Table 3, the component row lists the horizontal component of earthquakes selected in this study.  The FAST code [43] can input the ground motion for the tower base through the damped oscillator or large mass method. For the damped oscillator method, spring and damper were added to the tower base, and then the system was considered as a damped oscillator. The actuator frequency should be approximately 10 times the highest frequency of the turbine model when excited, which indicates that the time step was quite small and the numerical simulation might be unstable if the input ground motion includes impulse. Therefore, the large mass method [44] was employed to input the ground motion for the tower base in the present study. An artificially large mass is added to the tower base according to the large mass method, and a concentrated force is applied at the tower base to produce the desirable ground motion: where M is the artificially large mass; m represents the mass of NREL 5MW wind turbine and is 697 460 kg; and, a(t) is the acceleration time history of the input ground motion. Referring to the literature [44,45], the large mass M is set as 7 × 10 9 , which is 10,000 times larger than the mass of NREL 5MW wind turbine.

Coupled and Uncoupled Methods
As shown in Figure 3, the rotor of the wind turbine is subjected to turbulent wind and rotating with an angular velocity ω. When the rotor oscillates in the FA direction of the wind turbine, the blade element with radius r and length dr experiences an aerodynamic loading dT, which can be expressed as: where ρ is the density of air; C L and C D are the lift coefficients and drag coefficients, respectively; ϕ is the local inflow angle; α is the local angle of attack (AOA); β P (t) is the full-span pitch angle; and κ(t) is the aerodynamic twist of blades; a and a are the tangential induction factor and axial induction factor, respectively; V W is the wind speed; . u is the oscillation velocity of the blade in the FA direction of the wind turbine.
According to Equations (5)-(9), the oscillation velocity of the blade influences the aerodynamic loadings on the rotor, such that the seismic response of wind turbines subjected to wind and earthquake loads is a typical aeroelastic problem. Moreover, modern wind turbines widely adopted variable speed and pitch regulation technology [46]. Therefore, the aero-servo-elastic coupling effect should be considered when predicting the dynamic response of wind turbines. In the coupled method, Equqtions (5)-(9) should be solved with the equation of motion of the dynamic system which is derived using the Kane method in the FAST code. The analysis results obtained by the coupled method are taken as the benchmark to evaluate the accuracy of the uncoupled method.
The total oscillation velocity of the blade . u can be divided into: where . u w and . u E are the vibration velocity of the blade induced by wind and earthquake, respectively. For the large-scale wind turbines, the oscillation velocity of the blade caused by wind can directly be compared with the variations of inflow speed [25]. Assuming . u E is small compared to the wind speed, the aerodynamic loadings on the blade can be decoupled through a first-order Taylor expansion: where dT| .
u= . u w represents the aerodynamic loadings on the blade excited by wind; − u E is the increment of the aerodynamic loading caused by earthquakes and indicates damping effect on the wind turbine.
where ρ is the density of air; CL and CD are the lift coefficients and drag coefficients, respectively; φ 172 is the local inflow angle; α is the local angle of attack (AOA); P ( ) is the full-span pitch angle; and ( ) is the aerodynamic twist of blades; ′ and are the tangential induction factor and axial 174 induction factor, respectively; VW is the wind speed; ̇ is the oscillation velocity of the blade in the 175 FA direction of the wind turbine.

176
According to Equations (5)-(9), the oscillation velocity of the blade influences the aerodynamic 177 loadings on the rotor, such that the seismic response of wind turbines subjected to wind and 178 earthquake loads is a typical aeroelastic problem. Moreover, modern wind turbines widely adopted 179 variable speed and pitch regulation technology [46]. Therefore, the aero-servo-elastic coupling effect

184
The total oscillation velocity of the blade ̇ can be divided into: where ̇w and ̇E are the vibration velocity of the blade induced by wind and earthquake,   al. [30] and utilized by other researchers [31][32][33][34][35][36][37] was based on Equation (11). In this uncoupled 199 method, the responses of wind turbines induced by wind and earthquake, respectively, are calculated al. [30], it is assumed that the aerodynamic damping of the first two FA tower modes of wind turbines 212 is equal, which was also adopted by other researchers [31][32][33].    By Equation (11), the dynamic response of wind turbines excited by wind and earthquake can be decoupled. Though not mentioned in their papers, the uncoupled method proposed by Asareh et al. [30] and utilized by other researchers [31][32][33][34][35][36][37] was based on Equation (11). In this uncoupled method, the responses of wind turbines induced by wind and earthquake, respectively, are calculated initially, and they are linearly combined subsequently. The rotor is spinning under the wind excitation, while it is in the parked state for earthquake excitation only. In order to represent the wind-earthquake interaction, an equivalent aerodynamic damping ratio should be explicitly added to the wind turbines in the uncoupled method [47]. However, there has been a debate on the aerodynamic damping of wind turbines until now [47][48][49]. Consequently, the aerodynamic damping ratios in the uncoupled method are selected in a wide range for two purposes. Firstly, it is to cover the real value of aerodynamic damping. Secondly, it is attempted to determine the optimal aerodynamic damping ratio, which can maintain the accuracy of the uncoupled model for all the earthquakes and response quantities. Therefore, the equivalent aerodynamic damping ratio varies within the interval 0-10% at steps equal to 1.0% in this study. Moreover, the existing aerodynamic damping theory only considers the first FA tower mode. According to the suggestions of Asareh et al. [30], it is assumed that the aerodynamic damping of the first two FA tower modes of wind turbines is equal, which was also adopted by other researchers [31][32][33].

Vibration and Aerodynamic Loadings on the Rotor
The aeroelasticity of wind turbines is associated with the oscillation velocity of the blade in the FA direction. First, the oscillation velocity of blades under the excitation of wind only and wind-earthquake combination was analyzed using the fully coupled method where the average wind speed of hub-height was set to 11.4 m/s. Next, the influence of earthquakes on the aerodynamic loadings on the rotor were evaluated. Considering the symmetry of blades, the kinematic analysis was performed for the blade 1. In all the following simulations, the duration and time step were set as 600 s and 0.002 s, respectively.

Excited by Wind
As shown in Figure 3, two coordinate systems were introduced to analyze the vibration of the blade. The x, y, z coordinate system fixed at the ground is a static frame of reference, with the x-axis increasing downwind along the FA direction of the wind turbine. The ξ, η, γ coordinate system, fixed on the nacelle, translates and rotates concerning the x, y, z system. The coordinate of the origin of this moving frame is (0, 0, 90) in the x, y, z system. According to the kinematic theory [50], the absolute velocity of the blade along the x-axis, . u, can be expressed as: Energies 2020, 13, 3833 where V rx is the relative velocity of the blade along the x-axis concerning the ξ, η, γ frame; V NA,x is the absolute velocity of the origin of ξ, η, γ frame along the x-axis, ω NA,y is the angular velocity of ξ, η, γ frame about the y-axis, r is the radius of the section from the hub, θ is the azimuth angle of the blade. According to their definition, V NA,x and ω NA,y are the translational velocity in the x-axis and angular velocity about the y-axis of the nacelle, respectively. In Equation (12), the term V NA,x + ω NA,y × rcosθ represents the motion of ξ, η, γ frame observed from the fixed x, y, z frame and is referred to as transport velocity. The relative and transport velocity are caused by the deformation of the blade and support structure, respectively. Consequently, the absolute velocity of the blade is associated with the deformation of the blade and support structure. In the existing aerodynamic damping models of wind turbine towers [47][48][49], the blades were modeled as rigid and it is assumed that the motion of the entire rotor is consistent with that of the tower top but not rotation. According to Equation (12), the absolute velocity of the rigid blade along the x-axis is identical to the transport velocity. Therefore, only if the angular velocity of the nacelle is small will the oscillation velocity of the blade meet the assumption in aerodynamic damping models. Figure 4 is the time history of wind speed at the hub height. When the NREL 5 MW wind turbine was excited by this wind sample, the translational velocity V NA,x and angular velocity ω NA,y of the nacelle were recorded and are shown in Figure 5, and their amplitudes were 0.15 m/s and 0.003 rad/s, respectively. Obviously, the oscillation velocity of the nacelle excited by this wind sample was small enough to guarantee the serviceability of the wind turbine.

9
where rx is the relative velocity of the blade along the x-axis concerning the ξ, η, γ frame; NA,x is 230 the absolute velocity of the origin of ξ, η, γ frame along the x-axis, NA,y is the angular velocity of ξ, 231 η, γ frame about the y-axis, r is the radius of the section from the hub, is the azimuth angle of the 232 blade. According to their definition, NA,x and NA,y are the translational velocity in the x-axis and 233 angular velocity about the y-axis of the nacelle, respectively. In Equation (12), the term NA,x + 234 NA,y × represents the motion of ξ, η, γ frame observed from the fixed x, y, z frame and is 235 referred to as transport velocity.

236
The relative and transport velocity are caused by the deformation of the blade and support 237 structure, respectively. Consequently, the absolute velocity of the blade is associated with the 238 deformation of the blade and support structure. In the existing aerodynamic damping models of 239 wind turbine towers [47][48][49], the blades were modeled as rigid and it is assumed that the motion of 240 the entire rotor is consistent with that of the tower top but not rotation. According to Equation (12)    where rx is the relative velocity of the blade along the x-axis concerning the ξ, η, γ frame; NA,x is 230 the absolute velocity of the origin of ξ, η, γ frame along the x-axis, NA,y is the angular velocity of ξ, η, γ frame about the y-axis, r is the radius of the section from the hub, is the azimuth angle of the 232 blade. According to their definition, NA,x and NA,y are the translational velocity in the x-axis and 233 angular velocity about the y-axis of the nacelle, respectively. In Equation (12), the term NA,x + 234 NA,y × represents the motion of ξ, η, γ frame observed from the fixed x, y, z frame and is 235 referred to as transport velocity.

236
The relative and transport velocity are caused by the deformation of the blade and support 237 structure, respectively. Consequently, the absolute velocity of the blade is associated with the 238 deformation of the blade and support structure. In the existing aerodynamic damping models of 239 wind turbine towers [47][48][49], the blades were modeled as rigid and it is assumed that the motion of 240 the entire rotor is consistent with that of the tower top but not rotation. According to Equation (12)

253
The oscillation velocity of the blade in the x-axis excited by this wind sample is displayed in 254 Figure 6. As shown in Figure 6a, the transport velocity of the blade tip was less than 0.2 m/s, which The oscillation velocity of the blade in the x-axis excited by this wind sample is displayed in Figure 6. As shown in Figure 6a, the transport velocity of the blade tip was less than 0.2 m/s, which indicates the support structure has little contribution to the blade-tip velocity. From Figure 6b, the amplitude of relative velocity for the blade tip was 3.5 m/s, which was much larger than the transport velocity. Consequently, the absolute velocity of the blade tip depended on the deformation of the blade, and it was similar to the relative velocity of blade tip. The amplitudes of oscillation velocity for the blade are shown in Figure 6d, where V ax and V ex are the amplitude of the absolute velocity and transport velocity, respectively. As the angular velocity of the nacelle was small, the transport velocity amplitude of the blade was almost constant, which meets the assumption of the vibration velocity of the blade in the existing aerodynamic damping theory for wind turbine towers.
The version 13 of AeroDyn_code (NREL, Golden, CO, USA) [51], the aerodynamic module of FAST, output the aerodynamic loadings on blade 1. However, the local wind speed on the rotor disk was different because of the wind shear, which led to the aerodynamic loadings on the three blades not being identical. The necessary modifications were made for AeroDyn to output the aerodynamic loadings on the three blades, and the aerodynamic thrust, F, was the resultant force of aerodynamic loadings on the three blades in the x-axis: where l is the index of the blade. As shown in Figure 7, the amplitude of thrust in this example was 0.8 MN and the mean was 0.6 MN, where the aeroelastic effect was considered.
where l is the index of the blade. As shown in Figure 7, the amplitude of thrust in this example was   subjected to a single wind-earthquake event was analyzed to obtain preliminary insight into its 280 dynamic behavior. Considering that the spectral characteristics of seismic records 3 and 4 in Table 3 the period of the first and second FA tower modes, respectively. The seismic record 3 had abundant

Excited by Wind and Earthquake
The dynamic response of the wind turbine experiencing wind and earthquake loads is discussed using the coupled method in this section. To eliminate the initial transient behavior, the earthquake began at the 400th second in the following sections. The seismic response of the wind turbine subjected to a single wind-earthquake event was analyzed to obtain preliminary insight into its dynamic behavior. Considering that the spectral characteristics of seismic records 3 and 4 in Table 3 are typical, they were taken as the input ground motions, respectively. Their acceleration time history, response spectrum, and power spectrum are shown in Figure 8, where T 1 and T 2 represent the period of the first and second FA tower modes, respectively. The seismic record 3 had abundant high-frequency content, where the peak period of acceleration response spectrum was approximately identical to the second FA tower mode period, i.e., T 2 . Therefore, the high-order modes of the wind turbine may be significantly stimulated under the excitation of seismic record 3. For seismic record 4, its peak period of the power spectrum was close to the fundamental period of tower. Combining the two seismic records with the wind field in Section 3.1.1, respectively, the seismic response of the wind turbine was analyzed.  subjected to a single wind-earthquake event was analyzed to obtain preliminary insight into its 280 dynamic behavior. Considering that the spectral characteristics of seismic records 3 and 4 in Table 3 281 are typical, they were taken as the input ground motions, respectively. Their acceleration time  First, the oscillation velocity of the nacelle is shown in Figure 9 when the input ground motion is seismic record 3. The translational velocity amplitude of the nacelle along the x-axis was 0.5 m/s which was 2.5 times larger than that induced by the wind only. The angular velocity amplitude of the nacelle about the y-axis was 0.07 rad/s, which was 20 times larger than that of wind turbines under wind excitation. Therefore, the translational velocity and angular velocity of the nacelle were significantly increased by the earthquake load due to the deformation of the support structure.
The oscillation velocity of the blade in the x-axis is shown in Figure 10. From Figure 10a, the relative velocity amplitude of the blade tip was 4 m/s. According to Figure 10b, the absolute velocity amplitude of the blade tip was 6 m/s, which was larger than that excited by wind only. For the wind turbine at operational state, the increment of oscillation velocity of the blade caused by the earthquake . u E can be computed by u w based on Equation (10). The oscillation velocity increment of the blade tip is shown in the Figure 10c and its amplitude was 4.6 m/s, which was even larger than the velocity amplitude induced by wind. The oscillation velocity amplitudes of the blade are compared in Figure 10d, where V Wx represents the velocity amplitude of the blade excited by wind and ∆V represents the amplitude of blade velocity increment. The amplitude of blade velocity increment in the FA direction was non-monotonic, which was evidence that higher modes of blades are stimulated by seismic record 3. Moreover, the distribution of the transport velocity amplitude of the blade was in conflict with the assumption in the aerodynamic damping theory of wind turbine towers.

293
First, the oscillation velocity of the nacelle is shown in Figure 9 when the input ground motion 294 is seismic record 3. The translational velocity amplitude of the nacelle along the x-axis was 0.5 m/s 295 which was 2.5 times larger than that induced by the wind only. The angular velocity amplitude of 296 the nacelle about the y-axis was 0.07 rad/s, which was 20 times larger than that of wind turbines under 297 wind excitation. Therefore, the translational velocity and angular velocity of the nacelle were 298 significantly increased by the earthquake load due to the deformation of the support structure.

299
The oscillation velocity of the blade in the x-axis is shown in Figure 10. From Figure 10a

293
First, the oscillation velocity of the nacelle is shown in Figure 9 when the input ground motion 294 is seismic record 3. The translational velocity amplitude of the nacelle along the x-axis was 0.5 m/s 295 which was 2.5 times larger than that induced by the wind only. The angular velocity amplitude of 296 the nacelle about the y-axis was 0.07 rad/s, which was 20 times larger than that of wind turbines under 297 wind excitation. Therefore, the translational velocity and angular velocity of the nacelle were 298 significantly increased by the earthquake load due to the deformation of the support structure.

299
The oscillation velocity of the blade in the x-axis is shown in Figure 10. From Figure 10a, the  x-axis is shown in Figure 12. From Figure 12a, the relative velocity amplitude of the blade tip was 323 5m/s, which was larger than that caused by the wind. According Figure 12b, the absolute velocity    Next, the oscillation velocity of the nacelle is shown in Figure 11, where the input ground motion was seismic record 4. According to Figure 11a, the translational velocity amplitude of the nacelle along the x-axis was 1.5 m/s, which was 10 times larger than that induced by the wind only. From Figure 11b, the angular velocity amplitude of the nacelle about the y-axis was 0.03 rad/s which was also 10 times larger than that induced by the wind only. The vibration velocity of the blade along the x-axis is shown in Figure 12. From Figure 12a, the relative velocity amplitude of the blade tip was 5m/s, which was larger than that caused by the wind. According Figure 12b, the absolute velocity amplitude of the blade tip was 5.9 m/s. The increment of oscillation velocity of the blade tip is shown in Figure 12c and the amplitude was 4.5 m/s. According to Figure 12d, the amplitude of the blade velocity increment along the x-axis increased monotonously from the root to the tip, and it was larger than the blade velocity induced by the wind. Similar to the former example, the transport velocity of the blade induced by the deformation of the support structure was in conflict with the assumption in the aerodynamic damping theory of wind turbine towers.

317
Next, the oscillation velocity of the nacelle is shown in Figure 11, where the input ground motion  x-axis is shown in Figure 12. From Figure 12a, the relative velocity amplitude of the blade tip was 323 5m/s, which was larger than that caused by the wind. According Figure 12b, the absolute velocity

335
Considering the randomness of wind, five wind samples were generated by using the TurbSim 336 code with different random seeds and then combined with all the seismic records in Table 3. To   respectively. The velocity increment at the middle section of the blade was even larger than the 345 oscillation velocity induced by the wind for almost all the earthquake records, and the blade-tip 346 velocity increment factor was larger than 50% for most earthquakes. Therefore, the oscillation velocity 347 of the blade for the wind turbine at the operational state induced by some earthquakes did not meet 348 the requirement in Equation (11), which is the base of the uncoupled method. Considering the randomness of wind, five wind samples were generated by using the TurbSim code with different random seeds and then combined with all the seismic records in Table 3. To evaluate the influence of earthquakes on the vibration of the blade, the velocity increment factor ξ is defined as.
where j represents the index of wind sample; . u − . u W j,max is the amplitude of blade velocity increment along the x-axis caused by earthquakes when the wind field is sample j; . u W j,max represents the amplitude of blade velocity about the x-axis induced by the wind sample j, and m = 5 denotes the sample size of the wind field. Figure 13a,b represent the velocity increment factor for the middle and tip of the blade, respectively. The velocity increment at the middle section of the blade was even larger than the oscillation velocity induced by the wind for almost all the earthquake records, and the blade-tip velocity increment factor was larger than 50% for most earthquakes. Therefore, the oscillation velocity of the blade for the wind turbine at the operational state induced by some earthquakes did not meet the requirement in Equation (11), which is the base of the uncoupled method.

360
For the earthquakes in Table 3

Aerodynamic Loadings on the Rotor
For the wind turbine at the operational state, the oscillation velocity of blades in the FA direction is significantly changed by earthquakes, which will directly influence the aerodynamic loadings on the rotor. To assess this influence, the thrust variation is defined as: where F W and F W,E are the aerodynamic thrust on the rotor excited by wind only and combined wind-earthquake load. When the earthquake was seismic record 3, the thrust increment on the rotor of wind turbines is shown in Figure 14a, and the amplitude was 37 kN. When the earthquake was seismic record 4, the thrust increment is shown in Figure 14b, and the amplitude was 270 kN.   is significantly changed by earthquakes, which will directly influence the aerodynamic loadings on 354 the rotor. To assess this influence, the thrust variation is defined as: where W and W,E are the aerodynamic thrust on the rotor excited by wind only and combined 356 wind-earthquake load. When the earthquake was seismic record 3, the thrust increment on the rotor 357 of wind turbines is shown in Figure 14a, and the amplitude was 37 kN. When the earthquake was 358 seismic record 4, the thrust increment is shown in Figure 14b, and the amplitude was 270 kN.

360
For the earthquakes in Table 3 For the earthquakes in Table 3, the thrust variation factor of the wind turbine at the operational state is defined as: (16) where j, and m have the same meaning and value as those of Equation (14); F W,E ( j) − F W ( j) max is the amplitude of thrust variation caused by earthquakes when the wind field is sample j; F W ( j) max represents the amplitude of thrust induced by wind sample j.
The thrust variation factor shown in Figure 15a was larger than 0.15 for more than half of the earthquakes, and its maximum was 0.54. The relationship of thrust variation factor with respect to the amplitude of blade root velocity is examined in Figure 15b which indicates that the effect of earthquakes on the thrust increases with the amplitude of blade-root velocity. The aerodynamic loadings on the rotor were significantly changed by earthquakes. Consequently, the interaction of wind and earthquake loads must be considered in the coupled and uncoupled analysis. For the uncoupled model, the thrust variation induced by the earthquake was replaced by the equivalent aerodynamic damping.

16
The thrust variation factor shown in Figure 15a was larger than 0.15 for more than half of the 366 earthquakes, and its maximum was 0.54. The relationship of thrust variation factor with respect to 367 the amplitude of blade root velocity is examined in Figure 15b

384
The equivalent aerodynamic damping ratio for the uncoupled method is set to 4%, referring to 385 the literature [31][32][33], and the wind field is the same as that of Section 3.1. For the first example, the 386 input ground motion was the seismic record 3, and the responses were compared in Figure 16. The 387 differences between the two methods were significant for tower-top acceleration and tower-base 388 shear force and bending moment. Their amplitudes are listed in Table 4 where the tower-base

Comparisons of the Coupled and Uncoupled Methods
The analysis results of the coupled and uncoupled methods were compared to evaluate the accuracy of the uncoupled analysis method predicting the seismic response of wind turbines at the operational state. First, the response time-history of the two methods was compared for typical seismic records; next, the response amplitudes using the two methods were compared for different mean wind speed at hub height. Both of the two methods were implemented in the time domain by using the FAST code, where the damping ratios of blade and support structure modes were set as 0.4775% and 1% respectively.

Response Time History
The equivalent aerodynamic damping ratio for the uncoupled method is set to 4%, referring to the literature [31][32][33], and the wind field is the same as that of Section 3.1. For the first example, the input ground motion was the seismic record 3, and the responses were compared in Figure 16. The differences between the two methods were significant for tower-top acceleration and tower-base shear force and bending moment. Their amplitudes are listed in Table 4 where the tower-base bending-moment amplitudes were 120 MN·m and 90 MN·m for the coupled and uncoupled models, respectively. On the whole, the uncoupled model underestimates the amplitudes of the tower-top motion and the tower-base internal forces in this example.

395
Next, the input ground motion is changed to seismic record 4 in Table 3, and the responses of 396 the wind turbine are shown in Figure 17. The response amplitudes are listed in Table 5 where the 397 response amplitudes for the coupled method were less than those for the uncoupled method.

398
Therefore, the uncoupled method overestimates the tower-top motion and the tower-base internal 399 forces in this example. Table 4. Response amplitude of the wind turbine as the earthquake is seismic record 3. Next, the input ground motion is changed to seismic record 4 in Table 3, and the responses of the wind turbine are shown in Figure 17. The response amplitudes are listed in Table 5 where the response amplitudes for the coupled method were less than those for the uncoupled method. Therefore, the uncoupled method overestimates the tower-top motion and the tower-base internal forces in this example. Table 4. Response amplitude of the wind turbine as the earthquake is seismic record 3.

395
Next, the input ground motion is changed to seismic record 4 in Table 3, and the responses of 396 the wind turbine are shown in Figure 17. The response amplitudes are listed in Table 5 where the 397 response amplitudes for the coupled method were less than those for the uncoupled method.

398
Therefore, the uncoupled method overestimates the tower-top motion and the tower-base internal 399 forces in this example. 400 Table 4. Response amplitude of the wind turbine as the earthquake is seismic record 3. 404 Table 5. Response amplitude of the wind turbine as the earthquake is seismic record 4.     Figure 18 compares the response amplitudes of the tower using the coupled and uncoupled methods, where E3 and E4 denote the responses induced by seismic records 3 and 4, respectively. It was observed that the differences between them were significant besides the two sections discussed above. Therefore, there was remarkable discrepancy between the results of these two methods for seismic records 3 and 4. Moreover, when the equivalent aerodynamic damping ratio was 4%, the uncoupled model respectively underestimated and overestimated the seismic response of the wind turbine in the two examples, which indicates the conflict tendency for optimizing the aerodynamic damping ratio to improve the accuracy of the uncoupled method.     In order to evaluate the accuracy of the uncoupled method further, all the seismic records in 417 Table 3 were taken as input ground motions and the simulation results using the two methods were  Table 6 lists the mean wind speed at the hub and its sample size in the wind-earthquake load 429 combinations for this section.

Response Amplitude
In order to evaluate the accuracy of the uncoupled method further, all the seismic records in Table 3 were taken as input ground motions and the simulation results using the two methods were compared. The displacement error ζ d , acceleration error ζ a , shear-force error ζ F , and bending-moment error ζ d are defined as: where d and a are the tower-top displacement and acceleration amplitudes from the coupled model, whiled andã are the corresponding quantities from the uncoupled model; F and M are the tower-base shear-force and bending-moment amplitudes from the coupled model, whileF andM are the corresponding quantities for the uncoupled model. Therefore, if the error is larger than 0, the uncoupled model overestimates the seismic response of wind turbines. Conversely, the uncoupled model underestimates the seismic response. According to Chinese seismic standard [52], it is accepted that the uncoupled method has sufficient accuracy if the errors between the coupled and uncoupled methods are within the range of ±0.15, which is filled with grey in the following figures. Table 6 lists the mean wind speed at the hub and its sample size in the wind-earthquake load combinations for this section. When the mean wind speed was 11.4 m/s and the equivalent aerodynamic damping ratio was 4%, the errors between the coupled and uncoupled analysis were calculated and are shown in Figure 19. For some earthquakes, such as seismic records 31-36 in Table 3, the relative errors of the four response quantities were within the range of ±0.15. However, the relative errors were beyond the range of ±0.15 significantly for other earthquakes, which indicated this uncoupled method was not universal to predict the seismic response of wind turbines. The prominent conflict was that the remarkable underestimation and overestimation of the seismic responses emerged simultaneously for different earthquakes, which cannot be solved by just optimizing the equivalent aerodynamic damping ratios.
As shown in Figure 20, all the errors were less than 0.15 when the aerodynamic damping was increased to 7%. Compared with Figure 19, the overestimation of seismic responses was eliminated by increasing the aerodynamic damping, but some errors were much less than −0.15. The error of tower-top acceleration for seismic record 3 was −0.52, while it was −0.42 as the aerodynamic damping ratio was 4%. Consequently, the uncoupled model with a 7% aerodynamic damping was also not suitable to predict the seismic response of wind turbines as it indicated sufficient accuracy only for some earthquakes, e.g., seismic records 22-26, without universality. Figure 21 shows the relative errors between the coupled and uncoupled analysis, where the legend is the ID of the seismic event listed in Table 3. The relative errors decreased monotonously, and they were less than 0 for all the earthquakes when the equivalent aerodynamic damping ratio was 10%. For seismic record 3, the errors were less than 0 as the aerodynamic damping ratio was 1%, and the errors of acceleration, shear force, and bending moment were less than −0.15 once the aerodynamic damping ratio increased to 2%. However, the errors of all response quantities for seismic record 20 were larger than 0.4 while the aerodynamic damping ratio was 2%, and it should reach 7% to reduce the errors to 0.15. Therefore, the uncoupled model was not sufficient to accurately predict the seismic response of wind turbines when the mean wind speed at the hub height was 11.4 m/s.

432
When the mean wind speed was 11.4 m/s and the equivalent aerodynamic damping ratio was 433 4%, the errors between the coupled and uncoupled analysis were calculated and are shown in Figure   434 19. For some earthquakes, such as seismic records 31-36 in Table 3, the relative errors of the four 435 response quantities were within the range of ±0.15. However, the relative errors were beyond the 436 range of ±0.15 significantly for other earthquakes, which indicated this uncoupled method was not 437 universal to predict the seismic response of wind turbines. The prominent conflict was that the 438 remarkable underestimation and overestimation of the seismic responses emerged simultaneously 439 for different earthquakes, which cannot be solved by just optimizing the equivalent aerodynamic 440 damping ratios.

441
As shown in Figure 20, all the errors were less than 0.15 when the aerodynamic damping was 442 increased to 7%. Compared with Figure 19, the overestimation of seismic responses was eliminated 443 by increasing the aerodynamic damping, but some errors were much less than −0.15. The error of 444 tower-top acceleration for seismic record 3 was −0.52, while it was −0.42 as the aerodynamic damping 445 ratio was 4%. Consequently, the uncoupled model with a 7% aerodynamic damping was also not 446 suitable to predict the seismic response of wind turbines as it indicated sufficient accuracy only for 447 some earthquakes, e.g., seismic records 22-26, without universality. 448 Figure 21 shows the relative errors between the coupled and uncoupled analysis, where the 449 legend is the ID of the seismic event listed in Table 3. The relative errors decreased monotonously, 450 and they were less than 0 for all the earthquakes when the equivalent aerodynamic damping ratio 451 was 10%. For seismic record 3, the errors were less than 0 as the aerodynamic damping ratio was 1%,    When the average wind speed of hub-height was set to 5 m/s and the aerodynamic damping was increased to 3%, Figure 22 shows that the errors between the two models were less than 0.15. The errors for some earthquakes may be significantly less than −0.15, where the errors of tower-base shear force for seismic record 46 were the most remarkable and the minimum was −0.42. The errors between the coupled and uncoupled methods for different aerodynamic damping ratios were compared in Figure 23, where the overestimation of seismic responses was less than that illustrated When the average wind speed of hub-height was set to 5 m/s and the aerodynamic damping was increased to 3%, Figure 22 shows that the errors between the two models were less than 0.15. The errors for some earthquakes may be significantly less than −0.15, where the errors of tower-base shear force for seismic record 46 were the most remarkable and the minimum was −0.42. The errors between the coupled and uncoupled methods for different aerodynamic damping ratios were compared in Figure 23, where the overestimation of seismic responses was less than that illustrated in Figure 21 on the whole. Taking the seismic record 20, for example, when the aerodynamic damping was 4%, all the errors were nearly 0 as the average wind speed was 5 m/s; however, the errors were larger than 0.4 as the mean wind speed was 11.4 m/s. Consequently, the errors between the two methods were associated with the mean wind speed of hub height. When the mean wind speed was 5 m/s, a consistent aerodynamic damping ratio did not exist for the uncoupled model to accurately predict the seismic response of wind turbines. When the average wind speed of hub height was set to 18 m/s and the aerodynamic damping was 7%, Figure 24 shows that the relative errors between the coupled and uncoupled methods were not larger than 0.15. The errors for some earthquakes may be substantially less than −0.15, where the error of tower-base shear force for seismic record 46 was the most prominent and the minimum was −0.51. The errors between the coupled and uncoupled models for different aerodynamic damping ratios were compared in Figure 25. For a specified earthquake and response quantity, there was an equivalent aerodynamic damping ratio to maintain the high accuracy of the uncoupled model. However, for all the earthquakes in Table 3, a consistent aerodynamic damping ratio could not be determined to maintain the accuracy of the uncoupled model.
From Figures 20a, 22a and 24a, the errors of tower-top displacement were within the range of ±15% when the mean wind speed was 11.4 m/s, 5 m/s and 18 m/s. Therefore, just investigating the tower-top displacement was not sufficient to evaluate the uncoupled method. The aerodynamic damping ratio corresponding to the uncoupled model of Figures 20a, 22a and 24a was 7%, 3% and 7%, respectively. This was consistent with the conclusion that the aerodynamic ratio was associated with the mean wind speed. Consequently, a consistent aerodynamic damping ratio cannot be determined for the uncoupled method with different mean wind speeds.
The comparisons between the coupled and uncoupled methods indicate that this uncoupled method is not universal to analyze the seismic response of wind turbines. For the uncoupled method, the thrust variation caused by the ground motion is replaced by the equivalent aerodynamic damping ratio. Nevertheless, the distribution of the transport velocity of the blade in the FA direction excited by wind and earthquake is inconsistent with the assumptions in the existing aerodynamic damping models of wind turbine. Moreover, the existing aerodynamic damping model was only established for the first tower mode. As a result, the update of the aerodynamic damping model for wind turbines may be a feasible way to improve the accuracy of the uncoupled method. not larger than 0.15. The errors for some earthquakes may be substantially less than −0.15, where the error of tower-base shear force for seismic record 46 was the most prominent and the minimum was −0.51. The errors between the coupled and uncoupled models for different aerodynamic damping ratios were compared in Figure 25. For a specified earthquake and response quantity, there was an equivalent aerodynamic damping ratio to maintain the high accuracy of the uncoupled model. However, for all the earthquakes in Table 3, a consistent aerodynamic damping ratio could not be determined to maintain the accuracy of the uncoupled model.  From Figures. 20a, 22a and 24a, the errors of tower-top displacement were within the range of ±15% when the mean wind speed was 11.4 m/s, 5 m/s and 18 m/s. Therefore, just investigating the tower-top displacement was not sufficient to evaluate the uncoupled method. The aerodynamic damping ratio corresponding to the uncoupled model of Figures. 20a, 22a and 24a was 7%, 3% and 7%, respectively. This was consistent with the conclusion that the aerodynamic ratio was associated with the mean wind speed. Consequently, a consistent aerodynamic damping ratio cannot be determined for the uncoupled method with different mean wind speeds.
The comparisons between the coupled and uncoupled methods indicate that this uncoupled

Conclusions
The accuracy of the uncoupled method predicting the seismic response of wind turbines was investigated in this study. Firstly, the vibration of blades and aerodynamic loading on the rotor were analyzed to evaluate the assumptions in the uncoupled method. Subsequently, the simulation results of coupled and uncoupled models were compared to assess the accuracy of the uncoupled method. Based on the results, some conclusions are summarized in the following points.
(1) The oscillation velocity of the blade along the FA direction of wind turbines may be greatly influenced by the ground motions. The angular velocity of the nacelle induced by wind and earthquake load may be significantly larger than that induced by wind only, which illustrates that the transport velocity of the blade does not meet the assumption of aerodynamic damping models for wind turbine towers.
(2) The resultant forces of the aerodynamic loaings on the rotor may be significantly impacted by the ground motions selected in the present study, such that the interaction of wind and earthquake load is substantial. The influence of ground motions should be taken into account when computing the aerodynamic loadings on the rotor of wind turbines excited by the wind and earthquake combination.
(3) The errors between the coupled and uncoupled methods are related to both the mean wind speed at the hub height and input ground motions. The consistent aerodynamic damping ratios cannot be determined to maintain the accuracy of the uncoupled method for different wind speed and earthquakes. Therefore, this uncoupled method cannot be utilized to analyze the seismic response of wind turbines at its current state.
It should be noted that an aerodynamic damping model consistent with the vibration characteristic of the blade induced by an earthquake should be established to improve the uncoupled analysis method. This aerodynamic damping model should include the modal aerodynamic damping of higher modes of wind turbines and contribution of the rotational degree of freedom of the tower top. The examples in Section 4.1 show that the uncoupled method has different performance for ground motions with different spectral characteristics. However, the numerical simulation is carried out in the time domain due to the nonlinearity of the coupled system. Therefore, the interaction between wind and earthquake should be also discussed in the frequency domain.
Funding: This research was funded by the National Natural Science Foundation of China, grant number 51808061, 51722801 and 51678014. The APC was funded by the National Natural Science Foundation of China, grant number 51808061.

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