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

Abstract

There is a significant interaction between wind and earthquakes for large-scaled wind turbines due to an aeroelastic effect. 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 influence 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 different 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 significantly influenced 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 difficult to optimize a consistent aerodynamic damping in the uncoupled model to accurately predict the seismic response of wind turbines.

1. Introduction

As the primary form of wind energy utilization, wind power can fulfill the global energy demand and reduce CO₂ 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 the coupled analysis of wind and earthquake loads could be performed by this code [7]. Immediately, the reliability of this code was validated by seismic vibration table tests and numerical simulations [8]. As the FAST code could be developed by the users, it was widely used to explore the dynamic behavior of wind turbines experiencing wind and earthquake loads [9,10]. With this code, the dynamic response of the National Renewable Energy Laboratory (NREL) 5 MW wind turbine excited by wind and earthquake loads was investigated, and the results revealed the fundamental role of seismic loads in the reliability and economics of wind turbines [11]. The aerodynamic damping effect of wind on the seismic response of wind turbines was studied using the FAST code [12,13]. The coupled analysis was also performed by user-developed codes and these studies [14,15,16,17] emphasized the necessity of the coupled method predicting the seismic response of wind turbines. Furthermore, the interaction between wind and earthquake loads was also revealed by the model test [18]. However, the computational costs of the coupled method may be significant due to the randomness of wind and earthquakes, and the existing aeroelastic code for wind turbines had very limited ability to model the support structure and foundation. These shortcomings of the coupled method promoted the development of the uncoupled analysis method.

Actually, the uncoupled method was utilized to analyze the seismic response of wind turbines at the operational state earlier than the coupled method. In the preliminary stage, the aerodynamic loadings on the rotor were modeled as a static force acting at the tower top where a concentrated mass replaced the rotor-nacelle-assembly (RNA) [19]. With this method, seismic responses of wind turbines excited by wind and earthquake were analyzed [20,21], and the dominant role of earthquake load in high-seismic hazard regions was emphasized by using the deterministic analysis [22] and probabilistic analysis [23]. Though the interaction between wind and earthquake was not considered in these studies, the attention of researchers to the anti-seismic wind turbine was attracted. Subsequently, according to the BEM theory, the aerodynamic loadings on the rotor were divided into the mean and fluctuating thrust, where the latter was related to aeroelastic damping [24]. The wind–earthquake interaction was simulated by equivalent aerodynamic damping, and the simulation indicated that the response of the tower caused by a weak earthquake was more significant than that induced by wind loading.

As the large-scaled wind turbines are used in onshore and offshore wind farms, the blade is becoming more flexible. According to the reference [25], the maximum oscillation velocity of the blade tip normal to the chord will be as high as 6 m/s when the displacement amplitude is 1 m. This value can directly be compared with the variations of inflow speed, such that the aerodynamic loadings on the rotor should be computed by an aeroelastic tool [26]. For the multi-megawatt wind turbines, the comparison on structural models of the RNA indicated that higher modes were essential for the seismic response of the system [27]. As the lack of details for the blade, the beam element was usually adopted to model the rotor of wind turbines [28,29]. On the basis of these studies, an improved uncoupled method was proposed [30], where the aerodynamic forces of the blade nodes for each time step were computed by the AeroDyn code and the structural model of the wind turbine is established by the commercial finite element software ABAQUS. An equivalent aerodynamic damping ratio of 4% for the first two fore-aft (FA) tower modes was added to the structural damping to contain the interaction between wind and earthquake. This method was initially used to evaluate the seismic fragility of the NREL 5 MW wind turbine for different mean wind speeds. Immediately, this method was implemented in the GH BLADED software and assessed by comparisons between the coupled and uncoupled analysis results [31,32,33]. Besides, the uncoupled analysis was also performed in ABAQUS where aerodynamic loading was computed by the BEM theory and the average pitch angle was used to consider the effect of pitch regulation [34,35]. The uncoupled analysis method proposed by Asareh et al. determined the aerodynamic loading on the rotor using the aeroelastic code and could be directly implemented in the commercial software with high efficiency. Therefore, it has been widely used to analyze the seismic response of wind turbines [36,37].

Though comparisons between the coupled and uncoupled method were performed by Santangelo et al., the assumptions of the uncoupled method have not been illustrated and evaluated, and the input ground motions were limited. This study aims to further evaluate the uncoupled analysis method for the seismic response of wind turbines at operational condition. First, the numerical model of simulations and assumptions of the uncoupled method are illustrated in detail. Next, the oscillation velocity of blades and aerodynamic loading on the rotor are investigated to evaluate the interaction between wind and earthquake load. Finally, the seismic response of the NREL 5 MW wind turbine analyzed by the coupled and uncoupled methods are compared for different mean wind speed at hub height.

2. Numerical Modeling

The reference 5 MW wind turbine created by the NREL was taken as the research object, and its dimensions are shown in Figure 1. The main specifications of this turbine are listed in

Table 1. Properties of National Renewable Energy Laboratory (NREL) 5-MW Wind Turbine.

Part/Item	Property	Value
Blade	Rotor diameter	126 m
	Hub height	90 m
	Cut-in, rated and cut-out wind speed	3 m/s, 11.4 m/s, 25 m/s
	Cut-in and rated rotor speed	6.9 rpm, 12.1 rpm
	Length	61.5 m
	Overall mass	17,740 kg
	Structural damping ratio	0.5%
Hub and nacelle	Hub diameter	3 m
	Hub mass	56,780 kg
	Nacelle mass	240,000 kg
Tower	Bottom and top outer diameter	6 m, 3.87 m
	Bottom and top wall thickness	0.027 m, 0.019 m
	Overall mass	347,460 kg
	Structural damping ratio	1%

, and the details can be found in an existing report [38]. The natural frequencies of the NREL 5 MW wind turbine obtained from this design report are listed in

Table 2. Frequencies of the NREL 5 MW wind turbine.

Modes	Description	Natural Frequency (Hz)
1	1st Tower Fore-Aft	0.3240
2	1st Tower Side-to-Side	0.3120
3	1st Drivetrain Torsion	0.6205
4	1st Blade Asymmetric Flapwise Yaw	0.6664
5	1st Blade Asymmetric Flapwise Pitch	0.6675
6	1st Blade Collective Flap	0.6993
7	1st Blade Asymmetric Edgewise Pitch	1.0793
8	1st Blade Asymmetric Edgewise Yaw	1.0898
9	2nd Blade Asymmetric Flapwise Yaw	1.9337
10	2nd Blade Asymmetric Flapwise Pitch	1.9223
11	2nd Blade Collective Flap	2.0205
12	2nd Tower Fore-Aft	2.9003
13	2nd Tower Side-to-Side	2.9361

. The tower is divided into 50 beam elements of the same length. The blade shown in Figure 2 is discretized into 17 beam elements with the following length distribution: 2.73 m for the first three elements, 4.1 m for the next eleven elements, and 2.73 m for the remaining three elements. The FAST code [39] is used to analyze the seismic response of the NREL 5 MW wind turbine.

2.1. 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:

$$S_k\left(f\right)=\frac{4{\sigma }_k{L}_k/V_{\mathrm{hub}}}{{\left(1+6f{L}_k/V_{\mathrm{hub}}\right)}^{5/3}}$$

(1)

where f is the frequency (Hz); V_hub is the average wind speed of hub-height; ${\sigma }_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:

$$S_{i,j}\left(f\right)=C\left(\Delta r,f\right)\sqrt{S_{i,i}\left(f\right)·S_{j,j}\left(f\right)}$$

(2)

where S_i,i and S_j,j denote the wind velocity auto spectra at nodes i and j, respectively; $C\left(\Delta r,f\right)$ is the coherence function defined in Equation (3).

$$C\left(\Delta r,f\right)=\exp \left[-b\sqrt{\left(\frac{f\Delta r}{V_{\mathrm{hub}}}\right)+\left(\frac{0.12\Delta r}{L_{\mathrm{C}}}\right)}\right]$$

(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.

2.2. 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. Details of input ground motions obtained from Pacific Earthquake Engineering Research (PEER) database.

ID.	Event	Year	Station	Manitude	Component	PGA(g)
1	Kocaeli, Turkey	1999	Arcelik	7.5	ARC000	0.22
2	Duzce, Turkey	1999	Bolu	7.1	BOL000	0.73
3	Loma Prieta	1989	Capitola	6.9	CAP000	0.53
4	Chi-Chi, Taiwan	1999	CHY101	7.6	CHY101-E	0.35
5	Imperial Valley	1979	Delta	6.5	DLT262	0.24
6	Kocaeli, Turkey	1999	Duzce	7.5	DZC180	0.31
7	Imperial Valley	1979	El Centro Array #11	6.5	H-E11140	0.36
8	Loma Prieta	1989	Gilroy Array #3	6.9	G03090	0.37
9	Hector Mine	1999	Hector	7.1	HEC090	0.34
10	Superstition Hills	1987	El Centro Imp. Co.	6.5	B-ICC090	0.26
11	Northridge	1994	Canyon Country-WLC	6.7	LOS000	0.41
12	Northridge	1994	Beverly Hills-Mulhol	6.7	MUL009	0.42
13	Kobe, Japan	1995	Nishi-Akashi	6.9	NIS000	0.51
14	San Fernando	1971	LA-Hollywood Stor	6.6	PEL090	0.21
15	Superstition Hills	1987	Poe Road (temp)	6.5	POE360	0.30
16	Cape Mendocino	1992	Rio Dell Overpass	7.0	RIO270	0.39
17	Kobe, Japan	1995	Shin-Osaka	6.9	SHI000	0.24
18	Friuli, Italy	1976	Tolmezzo	6.5	TMZ000	0.31
19	Landers	1992	Yermo Fire Station	7.3	YER270	0.24
20	Manjil, Iran	1990	Abbar	7.4	ABBAR-T	0.21
21	Darfield	2010	Christchurch Cathedral College	7.0	CN26W	0.19
22	Chi-Chi, Taiwan	1999	CHY104	7.6	CHY104-N	0.17
23	Mexico	2010	Calexico Fire Station	7.2	CXO090	0.27
24	Mexico	2010	Cerro Prieto Geothermal	7.2	GEO000	0.26
25	Darfield	2010	Christchurch Hospital	7.0	HCS89W	0.15
26	Chi-Chi, Taiwan	1999	TCU070	7.6	TCU070-N	0.16
27	Chi-Chi, Taiwan	1999	TCU109	7.6	TCU109-N	0.16
28	Duzce, Turkey	1999	Duzce	7.1	DZC-180	0.30
29	Imperial Valley	1979	El Centro Array-6	6.5	H-E06230	0.44
30	Imperial Valley	1979	El Centro Array-7	6.5	H-E07140	0.34
31	Erzican, Turkey	1992	Erzincan	6.7	ERZ-NS	0.51
32	Kocaeli, Turkey	1999	Izmit	7.5	LZT090	0.22
33	Landers	1992	Lucerne	7.3	LCN260	0.73
34	Cape Mendocino	1992	Petrolia	7.0	PET090	0.66
35	Superstition Hills	1987	Parachute Test Site	6.5	PTS225	0.45
36	Northridge-01	1994	Rinaldi Receiving Sta	6.7	RRS228	0.83
37	Loma Prieta	1989	Saratoga-Aloha	6.9	STG090	0.32
38	Irpinia, Italy	1980	Sturno	6.9	STU270	0.36
39	Northridge	1994	Sylmar-Olive View	6.7	SYL360	1.68
40	Chi-Chi, Taiwan	1999	TCU065	7.6	TCU065-E	0.81
41	Chi-Chi, Taiwan	1999	TCU102	7.6	TCU102-E	0.30
42	Northridge-01	1994	LA-Sepulveda VA	6.7	0637-360	0.93
43	Imperial Valley-06	1979	Bonds Corner	6.5	BCR140	0.59
44	Loma Prieta	1989	BRAN	6.9	BRN000	0.48
45	Imperial Valley-06	1979	Chihuahua	6.5	CHI282	0.25
46	Loma Prieta	1989	Corralitos	6.9	CLS000	0.68
47	Gazli	1976	Karakyr	6.8	GAZ0	0.61
48	Nahanni	1985	Site 2	6.8	S2240	0.49
49	Nahanni	1985	Site 1	6.8	S1010	0.91
50	Northridge	1994	Northridge-Saticoy	6.7	STC090	0.37
51	Chi-Chi, Taiwan	1999	TCU067	7.6	TCU067-E	0.50
52	Chi-Chi, Taiwan	1999	TCU084	7.6	TCU084-E	1.16
53	Kocaeli, Turkey	1999	Yarimca	7.5	YPT330	0.35

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:

$$F_{\mathrm{a}}\left(t\right)=\left(M+m\right)·a\left(t\right)$$

(4)

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⁹, which is 10,000 times larger than the mass of NREL 5MW wind turbine.

2.3. 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 $\omega$. 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:

$$\mathrm{d}T=\frac{1}{2}\rho V_{\mathrm{rel}}^2\left[C_{\mathrm{L}}\left(\alpha \right)\cos \phi +C_{\mathrm{D}}\left(\alpha \right)\sin \phi \right]dr$$

(5)

$$V_{\mathrm{rel}}^2=V_{\mathrm{x}}^2+{\left[\omega r\left(1+a^{\prime }\right)\right]}^2$$

(6)

$$V_{\mathrm{x}}=V_{\mathrm{w}}\left(1-a\right)-\dot{u}$$

(7)

$$\phi ={\tan }^{-1}\left(\frac{V_{\mathrm{w}}\left(1-a\right)-\dot{u}}{\omega r\left(1+a^{\prime }\right)}\right)$$

(8)

$$\alpha =\phi -{\beta }_{\mathrm{P}}\left(t\right)-\kappa \left(t\right)$$

(9)

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); ${\beta }_{\mathrm{P}}\left(t\right)$ is the full-span pitch angle; and $\kappa \left(t\right)$ is the aerodynamic twist of blades; $a^{\prime }$ and $a$ are the tangential induction factor and axial induction factor, respectively; V_W is the wind speed; $\dot{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 $\dot{u}$ can be divided into:

$$\dot{u}={\dot{u}}_{\mathrm{w}}+{\dot{u}}_{\mathrm{E}}$$

(10)

where ${\dot{u}}_{\mathrm{w}}$ and ${\dot{u}}_{\mathrm{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 ${\dot{u}}_{\mathrm{E}}$ is small compared to the wind speed, the aerodynamic loadings on the blade can be decoupled through a first-order Taylor expansion:

$$\mathrm{d}T={\mathrm{d}T|}_{\dot{u}={\dot{u}}_{\mathrm{w}}}-\frac{\delta \left(\mathrm{d}T\right)}{\delta \dot{u}}·{\dot{u}}_{\mathrm{E}}$$

(11)

where ${\mathrm{d}T|}_{\dot{u}={\dot{u}}_{\mathrm{w}}}$ represents the aerodynamic loadings on the blade excited by wind; $-\frac{\delta \left(\mathrm{d}T\right)}{\delta \dot{u}}·{\dot{u}}_{\mathrm{E}}$ is the increment of the aerodynamic loading caused by earthquakes and indicates damping effect on the wind turbine.

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].

3. 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.

3.1. Oscillation Velocity of the Blade in the FA Direction

3.1.1. 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, $\dot{u}$, can be expressed as:

$$\dot{u}=V_{\mathrm{rx}}+\left(V_{\mathrm{NA},\mathrm{x}}+{\omega }_{\mathrm{NA},\mathrm{y}}\times rcos\theta \right)$$

(12)

where $V_{\mathrm{rx}}$ is the relative velocity of the blade along the x-axis concerning the ξ, η, γ frame; $V_{\mathrm{NA},\mathrm{x}}$ is the absolute velocity of the origin of ξ, η, γ frame along the x-axis, ${\omega }_{\mathrm{NA},\mathrm{y}}$ is the angular velocity of ξ, η, γ frame about the y-axis, r is the radius of the section from the hub, $\theta$ is the azimuth angle of the blade. According to their definition, $V_{\mathrm{NA},\mathrm{x}}$ and ${\omega }_{\mathrm{NA},\mathrm{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_{\mathrm{NA},\mathrm{x}}+{\omega }_{\mathrm{NA},\mathrm{y}}\times rcos\theta$ 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_{\mathrm{NA},\mathrm{x}}$ and angular velocity ${\omega }_{\mathrm{NA},\mathrm{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.

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:

$$F={\displaystyle \sum _{l=1}^3{\displaystyle {\int }_0^R\mathrm{d}T}}$$

(13)

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.

3.1.2. 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₁ and T₂ 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₂. 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.

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 ${\dot{u}}_{\mathrm{E}}$ can be computed by ${\dot{u}}_{\mathrm{E}}=\dot{u}-{\dot{u}}_{\mathrm{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.

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.

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.

$$\xi =\frac{\frac{1}{m}\times {\displaystyle \sum _{j=1}^m{\left|\dot{u}-{\dot{u}}_{\mathrm{W}}\right|}_{\mathrm{j},\max }}}{\frac{1}{m}\times {\displaystyle \sum _{j=1}^m{\left|{\dot{u}}_{\mathrm{W}}\right|}_{\mathrm{j},\max }}}$$

(14)

where j represents the index of wind sample; ${\left|\dot{u}-{\dot{u}}_{\mathrm{W}}\right|}_{\mathrm{j},\max }$ is the amplitude of blade velocity increment along the x-axis caused by earthquakes when the wind field is sample j; ${\left|{\dot{u}}_{\mathrm{W}}\right|}_{\mathrm{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.

3.2. 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:

$$\Delta F=F^{\mathrm{W},\mathrm{E}}-F^{\mathrm{W}}$$

(15)

where $F^{\mathrm{W}}$ and $F^{\mathrm{W},\mathrm{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.

For the earthquakes in Table 3, the thrust variation factor of the wind turbine at the operational state is defined as:

$$\eta =\frac{\frac{1}{m}\times {\displaystyle \sum _{j=1}^m{\left|F^{W,E}\left(j\right)-F^W\left(j\right)\right|}_{\max }}}{\frac{1}{m}\times {\displaystyle \sum _{j=1}^m{\left|F^W\left(j\right)\right|}_{\max }}}$$

(16)

where j, and m have the same meaning and value as those of Equation (14); ${\left|F^{W,E}\left(j\right)-F^W\left(j\right)\right|}_{\max }$ is the amplitude of thrust variation caused by earthquakes when the wind field is sample j; ${\left|F^W\left(j\right)\right|}_{\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.

4. 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.

4.1. 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. Response amplitude of the wind turbine as the earthquake is seismic record 3.

Model	Tower-Top Displacement	Tower-Top Acceleration	Tower-Base Shear-Force	Tower-Base Bending-Moment
	m	m/s2	MN	MN·m
coupled	0.53	1.61	3.04	120
uncouple	0.51	0.71	2.14	90

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.

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. Response amplitude of the wind turbine as the earthquake is seismic record 4.

Model	Tower-Top Displacement	Tower-Top Acceleration	Tower-Base Shear-Force	Tower-Base Bending-Moment
	m	m/s2	MN	MN·m
coupled	1.05	3.15	2.18	172
uncouple	1.25	4.03	2.49	211

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.

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.

4.2. 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:

$${\zeta }_{\mathrm{d}}=\frac{\tilde{d}-d}{d}$$

(17)

$${\zeta }_{\mathrm{a}}=\frac{\tilde{a}-a}{a}$$

(18)

$${\zeta }_{\mathrm{F}}=\frac{\tilde{F}-F}{F}$$

(19)

$${\zeta }_{\mathrm{M}}=\frac{\tilde{M}-M}{M}$$

(20)

where d and a are the tower-top displacement and acceleration amplitudes from the coupled model, while $\tilde{d}$ and $\tilde{a}$ 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, while $\tilde{F}$ and $\tilde{M}$ 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. Parameters of wind field.

Mean Wind Speed at Hub, Vhub	Wind Samples
5 m/s	5
11.4 m/s	5
18 m/s	5

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.

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 Figure 20a, Figure 22a and Figure 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 Figure 20a, Figure 22a and Figure 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.

5. 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.