An efficient rotational sampling method of wind fields for wind turbine blade fatigue analysis

Renewable Energy 146 (2020) 2170e2187

Contents lists available at ScienceDirect

Renewable Energy journal homepage: www.elsevier.com/locate/renene

An efficient rotational sampling method of wind fields for wind turbine blade fatigue analysis Jianbing Chen a, Yupeng Song b, Yongbo Peng c, *, Søren R.K. Nielsen d, Zili Zhang e a State Key Laboratory of Disaster Reduction in Civil Engineering & College of Civil Engineering, Tongji University, 1239 Siping Road, Shanghai, 200092, PR China b College of Civil Engineering, Tongji University, 1239 Siping Road, Shanghai, 200092, PR China c State Key Laboratory of Disaster Reduction in Civil Engineering & Shanghai Institute of Disaster Prevention and Relief, Tongji University, 1239 Siping Road, Shanghai, 200092, PR China d Department of Civil Engineering, Aalborg University, 9000, Aalborg, Denmark e Department of Engineering, Aarhus University, 8000, Aarhus C, Denmark

a r t i c l e i n f o

a b s t r a c t

Article history: Received 9 March 2019 Received in revised form 9 July 2019 Accepted 4 August 2019 Available online 7 August 2019

Both the operational and ultimate load conditions should be considered in the structural design and reliability assessment of wind turbine systems. In the operational condition, the fatigue load experienced by wind turbine blades is of great concern in design which highly relies upon the rotor’s rotation. Three kinds of methods have been developed to explore the rotational sampling effect of wind speeds on wind turbine blades, which, however, are somewhat inconvenient in practical applications. In view of the recent developments in wind field simulation, a novel rotational sampling method allowing for the analytical expression of fluctuating wind speeds on rotating blades is proposed in the present paper. In contrast to the existing methods, the proposed method circumvents the decomposition of cross power spectrum density (PSD) matrix and the interpolation in spatial and temporal dimensions. In particular, a closed-form expression of the rotational sampling spectrum is provided, thereby the mechanism of transfer of turbulent kinetic energy in frequency domain is quantitatively revealed. For illustrative purposes, fatigue analysis of the blades of a 5-MW offshore wind turbine is carried out, demonstrating the non-negligible influence of the rotational sampling on the fatigue load of blades and the competitive efficiency of the proposed method. © 2019 Elsevier Ltd. All rights reserved.

Keywords: Rotational sampling Wind turbine blades Fatigue analysis Stochastic harmonic function representation Wavenumber-frequency joint spectrum

1. Introduction As a family of infrastructures supporting green energy harvesting, the wind turbine systems have received increasing attention in recent years [1]. The operational and ultimate loading conditions are both necessary to be considered in the structural design and reliability assessment of wind turbine systems [2,3]. Unlike the standstill in the ultimate condition, the rotor rotates at a specified or adaptive speed/rate in the operational condition. The rotor’s rotation not only results in the centrifugal stiffening effect on the blades [4], but also causes the spectrum of wind speeds experienced by the rotating blade to be quite different from that by

* Corresponding author. E-mail addresses: [email protected] (J. Chen), [email protected] (Y. Song), [email protected] (Y. Peng), [email protected] (S.R.K. Nielsen), [email protected] (Z. Zhang). https://doi.org/10.1016/j.renene.2019.08.015 0960-1481/© 2019 Elsevier Ltd. All rights reserved.

a fixed blade, which is known as the rotational sampling effect [5]. Prediction of the operational fatigue loads in the lifecycle of wind turbine systems is one of the most important concerns for the design of a wind farm [2,6]. However, for a long-term service, the fatigue load experienced by wind turbine systems exhibits a high possibility of exceeding its expected value significantly due to the lack of awareness of the rotational sampling effect [7]. To consider the rotational sampling effect, three well-visited methods were proposed in 1980s [8], i.e., the Purdue method, the PNL (Pacific Northwest Laboratory) method and the Sandia method. In the Purdue method [9] many time series of wind speeds at specified positions in the rotor plane are firstly generated so as to construct a wind field. Then, as the blades rotate, the instantaneous wind speeds applied on the blades are obtained by picking up the wind speeds from the simulated wind field, which is known as the rotational sampling manipulation. In this procedure, interpolations in spatial and temporal dimensions are required since the wind

J. Chen et al. / Renewable Energy 146 (2020) 2170e2187

Nomenclature BEM CLD Cz ; Cy

blade element momentum constant lifetime diagram vertical and horizontal exponential decay coefficients in coherence model D fatigue damage Eð,Þ expectation operator fb rotational speed of rotor in Hz Hhub the height of the hub Jn ð,Þ the n-th order of Bessel function of the first kind kz ; ky vertical and horizontal wavenumbers PðUj Þ assigned probability of mean wind speed Uj r radial distance from the blade element to hub RSS rotational sampling spectrum SHF stochastic harmonic function S-N curve stress-cycle curve SRM spectral representation method SDav ðuÞ Davenport spectrum SKai ðz; uÞ Kaimal spectrum Sðkz ; ky ; uÞ Sðz; kz ; ky ; uÞjoint spectrum of homogeneous and nonhomogeneous wind speed field SR ðuÞ rotational sampling spectrum

field is constructed by discrete time series at finite spatial points. It is noted that in the Purdue method the wind field is modeled as a three-dimensional (3D) temporal-spatial random field. To improve the efficiency of wind field simulation, a 3D fast Fourier transform (FFT) technique is adopted based on the spectral representation method (SRM) [10]. However, the deduced joint spectrum is only capable of generating homogeneous wind fields and the extended coherence model has not yet been verified. Moreover, the 3D FFT technique requires a large amount of computer memory. Therefore, the Purdue method has not been widely-used and recommended in practice. In the PNL method [11], a celebrated rotational sampling spectrum based on the assumption of isotropic homogeneous turbulence was derived. In this method, the rotational effect is firstly represented through the von Karman correlation function and is then passed onto frequency domain by the Fourier transform. Thus, the rotational sampling fluctuating wind speeds can be generated by the SRM [12]. Along this way, He and Li derived a rotational Fourier spectrum based on the stochastic function description of random processes [13]. The influence of wind shear on the rotational Fourier spectrum was further investigated [14]. Similar to the PNL method, a more thorough approach was developed later [15], through which three components of the rotational sampling fluctuating wind speeds can be generated for the variable speed rotor. Although the PNL method is of rigorousness in derivation, it is difficult to obtain the analytical expression of the rotational sampling spectrum, which results in much inconvenience in applications. In addition, the assumption of isotropic homogeneous turbulence is limited in practice as well. Therefore, the PNL method has not been widely-used in practice either in recent years. The Sandia method [16,17], which differs from the Purdue method merely in the wind field simulation, is the most widelyused method at present [18e21]. This method is based on the SRM for multivariate stochastic processes [22]. It has been investigated over decades to simulate homogeneous and nonhomogeneous, stationary and nonstationary fluctuating wind fields [23]. However, the decomposition of the cross PSD matrix still remains at each discretized frequency, resulting in computational inefficiency

2171

Sh ðuÞ Pierson-Moscowitz (P-M) wave spectrum T time length of simulation Tlife lifetime of the blade UðzÞ; Uhub ; U10 ; U19:5 mean wind speed at the heights of z, hub, 10 m and 19.5 m Uðr; tÞ; uðr; tÞ time-varying mean wind and fluctuating wind experienced by the blade element with a radial distance r to hub Uj selected representative mean wind speed at hub for fatigue damage calculation u shear velocity Vj ; VolðVj Þ partitioned wavenumber-frequency subdomain and its volume yðr; tÞ; zðr; tÞ time-varying spatial coordinates of the blade element with a radial distance r to hub z0 roughness length ub rotational speed of rotor in rad/s k von Karman constant sm ; sa mean stress and stress amplitude rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðyÞ ~ðzÞ ; ~ ~ j Þ,VolðVj Þ kj ; u Aj ¼ 4Sðk j rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ðyÞ 2 ~ðzÞ Þ k Þ þ ðk B ¼ r ð~ j

j

j

for the large scale wind field simulation, or even numerical singularity [24]. Besides, the interpolation in spatial and temporal dimensions, which is also involved in the Purdue method, causes the Sandia method to be somewhat inconvenient in practice. Despite the preceding advances in modeling the rotational sampling effect, it is still somewhat cumbersome to apply these methods in practice. Consequently, the rotational sampling of fluctuating wind speeds in some cases is treated in an artificially approximate manner [5] or even not included in the investigation [25,26]. In order to explore a feasible approach to deal with the rotational sampling problem associated with the wind turbine analysis and to pursue an analytical expression of rotational sampling spectrum, a novel method for rotational sampling of fluctuating wind speeds is proposed based on the newly developed method of wind field simulation [27e29]. The remaining sections in this paper are organized as follows. Section 2 revisits the wind field simulation method associated with the evolutionary frequency-wavenumber joint spectrum and the stochastic harmonic function representation. Section 3 addresses the closed-form expression of the rotational sampling fluctuating wind speeds. Section 4 derives the analytical expression of the rotational sampling spectrum, based on which the mechanism of transfer of turbulent kinetic energy in frequency domain is clearly revealed. To evaluate the rotational sampling effect on the fatigue lifetime of wind turbine blades, numerical simulations of a 5-MW offshore wind turbine are carried out in Section 5. Concluding remarks pertaining to this study are provided in Section 6.

2. Simulation of fluctuating wind speed fields According to the IEC standards [2], two turbulence models are often recommended for the determination of design loads, i.e. the Mann uniform shear model [30], and the Kaimal spectrum and exponential coherence model. Due to its simplicity, the latter is used in the present paper for developing the novel rotational sampling method.

2172

J. Chen et al. / Renewable Energy 146 (2020) 2170e2187

2.1. Wavenumber-frequency joint spectrum of fluctuating wind speed fields Although the turbulent wind velocity consists of three components, i.e., the longitudinal, the lateral and the vertical components, it is acceptable to only consider the longitudinal component for preliminary analysis of the substructures of tower and foundations [31]. For clarity, the longitudinal turbulent wind velocity in the rotor plane is denoted as uðz;y;tÞ, where z; y denote the vertical and lateral spatial directions, respectively, and t denotes the time. Clearly, the fluctuating wind speed field uðz; y; tÞ is essentially a 3D temporal-spatial random field. Along this line, a wavenumberfrequency joint spectrum for the homogeneous fluctuating wind speed field in two-spatial dimensions was developed [27]. It was then extended to nonhomogeneous cases by introducing the concept of evolutionary spectrum [28]. Since the Kaimal spectrum describing nonhomogeneous fluctuating wind fields is usually recommended for the design of wind turbines [2], only the joint spectrum for nonhomogeneous cases is revisited here. In this case, the Kaimal spectrum and the two-dimensional Davenport’s coherence model are, respectively, given by Ref. [32]

SKai ðz; uÞ ¼

50zu2*  5=3 pUðzÞ 1 þ 2p50 j u jz UðzÞ 





ru1 u2 ðuÞ ¼ r xz ; xy ; u ¼ exp 

(1)

this method. In other words, fluctuating wind spectra such as the von Karman spectrum etc. [32], can also be adopted. Other coherence models, e.g., the exponential coherence model of Eq. (5) suggested in the IEC code [2], can be adopted as well







qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

r xz ; xy ; u ¼ exp  12 ½u=ð2pUhub Þ 2 þ ð0:12=Lc Þ2 ,jDrj  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ¼ exp  12 ½u=ð2pUhub Þ 2 þ ð0:12=Lc Þ2 , x2z þ x2y qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi i h ¼ exp  lðuÞ, x2z þ x2y

(5) where Uhub denotes the mean wind speed at the hub height; Lc denotes the coherence scale parameter that is often assumed to depend on the height z; Dr is the spatial distance between two qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi points in the rotor plane, Dr ¼ x2z þ x2y ; and lðuÞ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 12 ½u=ð2pUhub Þ 2 þ ð0:12=Lc Þ2 . Similar to the derivation procedure in the previous investigation [27], a two-fold Fourier transform of Eq. (5) with respect to xz ; xy yields





rðWFÞ kz ; ky ; u ¼

1



2pl ðuÞ 1 þ 2

1 j uj 2pU10

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  C 2z x2z þ C 2y x2y

(2)

then the wavenumber-frequency joint spectrum can be derived as [28]



1 l2 ðuÞ



k2z

þ

k2y

(6)

3=2

Correspondingly, the joint spectrum becomes

  S z; kz ; ky ; u ¼ SKai ðz; uÞ, 2



2pl ðuÞ 1 þ

1 1 l2 ðuÞ



k2z þ k2y

3=2

(7)

    S z; kz ; ky ; u ¼ SKai ðz; uÞ,rðWFÞ kz ; ky ; u 1 ¼ SKai ðz; uÞ, 2pCz Cy 

1 1 j uj 2pU10

2

" 1þ

1 ky Cy

2



1 2 #,

1 þ kz Cz

where SKai ðz; uÞ denotes the double-sided Kaimal spectrum; u is the circular frequency; u is the shear velocity; UðzÞ is the mean wind speed at the height z and follows the logarithmic law, i.e.,

UðzÞ ¼

u*

k

ln

z z0

(4)

in which k is the von Karman constant and z0 is the roughness length [32]; ru1 u2 ðuÞ denotes the coherence function; xz ; xy are the spatial distances, i.e. xz ¼ z1  z2 ,xy ¼ y1  y2 ; U10 denotes the mean wind speed at 10 m height; Cz ; Cy are the exponential decay coefficients in z; y directions, respectively; Sðz; kz ; ky ; uÞ denotes the joint spectrum; kz ; ky are the wavenumbers in z; y directions, respectively; and rðWFÞ ðkz ; ky ; uÞ is the two-fold Fourier transform of rðxz ; xy ; uÞ with respect to xz ; xy . It is noted that there is no restriction on the wind spectra used in

1 j uj 2pU10

2 !32

(3)

2.2. Stochastic harmonic function representation for fluctuating wind speed fields Once the joint spectrum of the fluctuating wind speed field is derived, the classical SRM for multi-dimensional processes [22] can be employed to generate samples of the fluctuating wind field. To improve the simulation efficiency, uneven discretization techniques by tensor-product and acceptance-rejection methods were proposed [27,28]. To further reduce the random variables as well as the computational burden, the stochastic harmonic function (SHF) representation was applied [29], which is also adopted in the present paper. In the SHF representation, both the phase angles and the discretized frequencies/wavenumbers are treated as random variables, which distinguishes it from the SRM where only the phase angles are random variables. The random fluctuating wind field is

J. Chen et al. / Renewable Energy 146 (2020) 2170e2187

2173

z

represented by

0

1

u@z; y; t A ¼

N X

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

ðzÞ ðyÞ

~ ;u ~ ;k ~ j ,Vol Vj  4S z; k j j

j¼1

ðzÞ

ðzÞ ~ðyÞ y þ u ~ðyÞ y  u ~ zþk ~ zþk ~ j t þ fjð1Þ þ cos k ~ j t þ fjð2Þ cos k j j j j

ðzÞ

i ðzÞ ~ zk ~ zk ~ðyÞ y þ u ~ðyÞ y  u ~ j t þ fjð3Þ þ cos k ~ j t þ fjð4Þ þcos k j j j j h

r

(8) ~ðzÞ ; k ~ðyÞ ; where N is the total number of the harmonic components; ðk j j ~ j Þ’s are independent 3D random vectors uniformly distributed u over the subdomain Vj ðj ¼ 1; 2;:::;NÞ. The subdomains Vj ’s are nonU overlapping such that Vj ∩Vm ¼ ∅; cjsm and ∪N j¼1 Vj ¼ ½0; u   ½0; U U U U U kz   ½0;ky , where u is the upper cut-off frequency, kz ; ky are the upper cut-off wavenumbers in z; y directions, respectively. To determine the subdomains Vj ’s, the Voronoi cell partitioning scheme is adopted after determining the representative points by the acceptance-rejection method [29,33]. Volð ,Þ denotes the volð1Þ ð2Þ ð3Þ ume (i.e., Lebesgue measure) of the subdomain. fj ’s, fj ’s, fj ’s ð4Þ and fj ’s are four different sets of independent random phase angles uniformly distributed over [0, 2p]. It is noticeable that Eq. (8) represents a zero-mean Gaussian field as N/∞. It can be seen that the SHF representation is an extension of the classical SRM. However, the number of the harmonic components involved in Eq. (8) is much less than that in the classical SRM [29]. Therefore, the computational efforts can be reduced, thus making it more feasible in practical applications. ðyÞ ~ðzÞ ; ~ ~ j Þ’s can also be taken as dependent It is noted that ðk kj ; u j random vectors but will not be detailed here to avoid lengthiness of the paper. More details about the SHF representation can be found in Refs. [29,33e35].

According to Section 2, the instantaneous fluctuating wind speeds at an arbitrary spatial point ðz; yÞ can be generated by Eq. (8). In the following, the rotational sampling wind speed of the blade element P with a radial distance of r is to be considered, as shown in Fig. 1. As the blade rotates, the spatial coordinates ðy; zÞ of the element P at the time t can be expressed as

8 <

yðr; tÞ ¼ r cosðqi ðtÞÞ : zðr; tÞ ¼ Hhub þ r sinðqi ðtÞÞ

(9)

where Hhub denotes the height of the hub, qi denotes the instantaneous azimuth between the i  th ði ¼ 1; 2; 3Þ blade and the horizontal axis at the time t. For the variable speed wind turbine system, qi ðtÞ can be determined by

ðt

qi ðtÞ ¼ qið0Þ þ 2pf  b ðtÞdt 0

o

(10)

y

Fig. 1. Schematic diagram of the rotor plane.

ð0Þ

where qi denotes the initial azimuth of the i  th blade as shown in Fig. 1. ~f b ðtÞ is the time-varying rotational speed (in Hertz) of the rotor. Substituting Eq. (9) into Eq. (8), the instantaneous fluctuating wind speed experienced by the element P at the time t can be obtained as,

1

u@r; t A ¼

3.1. Analytical expression of rotational sampling wind speeds

0

i

0

3. Rotational sampling method of wind speeds

P

h

N X

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

ðzÞ ðyÞ

~ ;~ ~ j ,Vol Vj  kj ; u 4S z r; t ; k j

j¼1



ðzÞ

ðyÞ ~ z r; t þ ~ ~ j t þ fð1Þ þ kj y r; t þ u cos k j j

ðzÞ



ðyÞ ð2Þ ~ y r; t  u ~ z r; t þ k ~ j t þ fj þ cos k j j



ðzÞ

ðyÞ ~ y r; t þ u ~ z r; t  k ~ j t þ fjð3Þ þ cos k j j

i ðzÞ



ðyÞ ~ z r; t  k ~ y r; t  u ~ j t þ fjð4Þ cos k j j

(11)

Eq. (11) provides the analytical expression of the rotational sampling fluctuating wind speeds for a wind turbine. Once the instantaneous azimuth of the blade is determined, the fluctuating wind speed experienced by an arbitrary element of the blade can be generated synchronously. Particularly, if the blade rotates in the clockwise direction at an invariant rate fb , Eq. (9) can be further written as

8

> < yðr; tÞ ¼ r cos  2pfb t þ qð0Þ i

ð0Þ > : zðr; tÞ ¼ Hhub þ r sin  2pfb t þ qi Correspondingly, Eq. (11) becomes

(12)

2174

J. Chen et al. / Renewable Energy 146 (2020) 2170e2187

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ! u N u

X   ð0Þ t  ðzÞ  ðyÞ  ; k j ; k j ; uj ,Vol Vj  uðr; tÞ ¼ 4S Hhub þ r sin  2pfb t þ qi (

"

j¼1

#



ð0Þ ð0Þ ð1Þ  ðzÞ  ðyÞ  k j þ r cos  2pfb t þ qi k j þ uj t þ fj cos Hhub þ r sin  2pfb t þ qi " #



ðzÞ ðyÞ ð0Þ ð0Þ ð2Þ k j þ r cos  2pfb t þ qi k j  u t þ f þcos Hhub þ r sin  2pfb t þ qi j j " #



ð0Þ ð0Þ ð3Þ  ðzÞ  ðyÞ  k j  r cos  2pfb t þ qi k j þ uj t þ f j þcos Hhub þ r sin  2pfb t þ qi " #)



ð0Þ ð0Þ ð4Þ  ðzÞ  ðyÞ  k j  r cos  2pfb t þ qi k j  uj t þ f j þcos Hhub þ r sin  2pfb t þ qi

It is noted that the mean wind speed UðzÞ involved in the Kaimal spectrum as well as in the joint spectrum also varies with the rotor rotation. Since UðzÞ follows the logarithmic law of wind profile, substituting Eq. (9) or (12) into Eq. (4) yields

Uðr; tÞ ¼

u*

k

ln

zðr; tÞ z0

(14)

which denotes the instantaneous “mean” wind speed, or more rigorously the “deterministic” part of wind speed, experienced by that element. It is noted that, the wind speed experienced by the element is no longer stationary due to rotational effects. However, it can be still decomposed into a deterministic part and a random part. That is, the rotational sampling of the total wind speed can be expressed as

~ tÞ ¼ Uðr; tÞ þ uðr; tÞ Uðr;

(15)

which can be directly utilized to calculate the aerodynamic loads on the rotating blades. Interestingly, the fluctuating wind field as a 3D temporal-spatial random field is viewed by both the proposed method and the Purdue method. However, the joint spectrum adopted in this paper

Vr ðr; tÞ ¼

derivation processes. Besides, in the Purdue method, the interpolation in spatial and temporal dimensions is still necessary to implement the rotational sampling on the basis of the provided time histories of wind speeds at a collection of spatial points. In contrast, the proposed method facilitates an easier way of generating the wind speeds experienced by the blades as the rotor rotates, which is of much convenience in practice. In addition, the computational efforts are acceptable for applications due to the adoption of the SHF representation instead of the classical SRM. 3.2. Applications in calculation of aerodynamic loads on blades To illustrate how the rotational sampling wind speed is applied in the analysis of wind turbines, the calculation of aerodynamic loads on blades is briefly introduced in this subsection. In the determination of the aerodynamic loads on wind turbine blades, the blade element momentum (BEM) method is the most widely-applied method [6], and is also adopted in the present paper. In the BEM method, each blade is divided into several elements. It is assumed that all the element sections are independent so that wind loads on each element can be calculated separately. Considering an air-foil of blades as shown in Fig. 2, the relative wind velocity applied on the air-foil at the time t is given by

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi _ tÞ 2 ð1  aðr; tÞ Þ2 þ ½2pfb rð1 þ a’ ðr; tÞ Þ 2 ½Uðr; tÞ þ uðr; tÞ  xðr;

is more rigorous in terms of theoretical foundations and the

2 fbr (1

(16)

_ tÞ denotes the velocity of the blade element due to where xðr; structural vibration in the direction of mean wind; aðr; tÞ and a’ðr; tÞ denote the axial and tangential induction factors, respectively, which can be determined by iterative algorithms. Thus, the angle of attack aðr; tÞ is obtained from

a) Rotor plane

U u xi 1 a

(13)

Vr

Fig. 2. Blade air-foil and the applied wind velocities.

aðr; tÞ ¼ 4ðr; tÞ  bðtÞ  gðrÞ

(17)

_ tÞ =½2pfb rð1 þ a’ Þ  g is where 4ðr; tÞ ¼ arctanf½Uðr; tÞ þ uðr; tÞ  xðr; the flow angle, bðtÞ is the pitch angle of the blade, and gðrÞ is the pre-twist angle at the cross section. Then, the local lift and drag forces on the element section can be determined, respectively, by

J. Chen et al. / Renewable Energy 146 (2020) 2170e2187

pN

2175

One can observe that the rotational sampling wind speeds in the procedure are just involved in Eq. (16). The flowchart of dynamic analysis of the blade system by means of the proposed method is shown in Fig. 4. Instead of generating wind speed field first and then carrying out interpolations in temporal and spatial dimensions, the rotational sampling wind speeds at any given time t can be obtained straightforwardly by the proposed method.

pD

pL

3.3. Variance of rotational sampling wind speeds

pT

The variance of fluctuating wind speed is one of the most critical factors that dominate the wind induced response of structures. In this subsection, the variance of rotational sampling wind speeds is investigated. For the simplicity of notation, define the following operational rules with respect to the 3D vectors a ¼ ða1 ; a2 ; a3 Þ and b ¼ ðb1 ;b2 ;b3 Þ

Rotor plane

V

8 > > a,þþ bba1 b1 þ a2 b2 þ a3 b3 > > < þ a, bba1 b1 þ a2 b2  a3 b3 þ > > a, bba1 b1  a2 b2 þ a3 b3 > > : a, bba1 b1  a2 b2  a3 b3

Fig. 3. Lift and drag forces applied on a blade air-foil.

8 > > > 1 > 2 > < pL ðr; tÞ ¼ ra V r ðr; tÞcðrÞCL ðr; aÞ 2 > 1 > > p ðr; tÞ ¼ ra V 2r ðr; tÞcðrÞCD ðr; aÞ > > : D 2

(18) in this way, Eq. (11) can be rewritten as

0

in which, ra denotes the density of air; c denotes the chord length; CL and CD are the lift and drag coefficients of the air-foil, respectively. Further, the wind loads on the cross section along the out-ofplane and in-plane directions can be obtained by the projection of pL and pD , as shown in Fig. 3,

8 <

pN ðr; tÞ ¼ pL cos 4 þ pD sin 4 : pT ðr; tÞ ¼ pL sin 4  pD cos 4

h



s2uðr;tÞ ¼ E u2 r; t

(20)

(19)

1

u@r; t A ¼

N qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

X    ffi h ð1Þ 4S zðr; tÞ; K j Vol Vj  cos K j ,þþ x þ fj

j¼1



ð2Þ ð3Þ þcos K j ,þ x þ fj þ cos K j ,þ x þ fj

i ð4Þ þcos K j , x þ fj (21) ðzÞ ðyÞ ~ ;u ~ j Þ and x ¼ ðzðr; tÞ; yðr; tÞ; tÞ are vectors. where K j ¼ ð~ kj ; k j ð1Þ ð2Þ ð3Þ ð4Þ Due to the independence among fj ’s, fj ’s, fj ’s, fj ’s and K j ’s, it is easy to obtain the variance of uðr; tÞ by

i

8 N qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi


   ffi h ð1Þ ð2Þ 4S zðr; tÞ; K j Vol Vj , cos K j ,þþ x þ fj þ cos K j ,þ x þ fj ¼E : j¼1



io2 ð3Þ ð4Þ þ cos K j , x þ fj þcos K j ,þ x þ fj N



X      h ð1Þ ð2Þ þ cos2 K j ,þ x þ fj ¼ E 4S zðr; tÞ; K j Vol Vj , cos2 K j ,þþ x þ fj j¼1





io ð3Þ ð4Þ þ cos2 K j , x þ fj þcos2 K j ,þ x þ fj N N ð



X     X E 8S zðr; tÞ; K j Vol Vj ¼ 8S zðr; tÞ; kj Vol Vj pK j kj dkj ¼ j¼1

¼

j¼1

N ð X j¼1



 8S zðr; tÞ; kj dkj ¼

Vj Integral domain/∞

ð

8Sðzðr; tÞ; kÞdk V0 þ∞ ð þ∞ ð þ∞ ð

! þ∞ ð þ∞ ð ð þ∞

¼ ∞ ∞ ∞ þ∞ ð ∞

  S zðr; tÞ; kz ; ky ; u dkz dky du

∞ ∞ ∞





SKai zðr; tÞ; u rðWFÞ kz ; ky ; u dkz dky du



SKai zðr; tÞ; u du

¼

(22)

Vj

2176

J. Chen et al. / Renewable Energy 146 (2020) 2170e2187

Initial azimuth of blades

uðr; tÞ ¼

N X

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðzÞ ðyÞ

  ~ j ,Vol Vj  kj ; u 4S ~ kj ; ~

j¼1

i n h ðzÞ ~ðyÞ cosðu tÞ þ u ~ sinðu tÞ þ r k ~ðzÞ ~ j t þ fjð1Þ þ Hhub k cos r k b b j j j i h ðzÞ ~ sinðu tÞ þ r k ~ðyÞ cosðu tÞ  u ~ðzÞ ~ j t þ fjð2Þ þ Hhub k þcos r k b b j j j i h ðzÞ ~ðzÞ ~ sinðu tÞ  r k ~ðyÞ cosðu tÞ þ u ~ j t þ fjð3Þ þ Hhub k þcos r k b b j j j io h ðzÞ ~ sinðu tÞ  r k ~ðyÞ cosðu tÞ  u ~ðzÞ ~ j t þ fjð4Þ þ Hhub k þcos r k b b j j j

Generation of rotational sampling wind speeds at the present moment: Eqs. (11)(14)(15)

Calculation of wind loads: Eqs. (16)-(19)

(24) For clarity, denote

Dynamic equations

i h ðzÞ ðyÞ ðzÞ ~ sinðu tÞ þ r ~ ~ j t þ fð1Þ f1;j ðtÞ ¼ cos r k kj cosðub tÞ þ u þ Hhub ~ kj b j j i h ðzÞ ðyÞ ðzÞ ~ sinðu tÞ þ r ~ ~ j t þ fð2Þ kj f2;j ðtÞ ¼ cos r k kj cosðub tÞ  u þ Hhub ~ b j j i h ðzÞ ðyÞ ðzÞ ~ sinðu tÞ  r ~ ~ j t þ fð3Þ kj f3;j ðtÞ ¼ cos r k kj cosðub tÞ þ u þ Hhub ~ b j j i h ðzÞ ðyÞ ðzÞ ~ sinðu tÞ  r ~ ~ j t þ fð4Þ kj kj cosðub tÞ  u þ Hhub ~ f4;j ðtÞ ¼ cos r k b j j

Responses of blades; Pitch angle of blades; Rotational speed of rotor; Azimuth of blades in the next step Fig. 4. Flow chart of dynamic analysis of blades.

(25) Then, by employing the Jacobi-Anger expansion [36].

where Eð,Þ denotes the expectation operator. It is noted that the integration of the Kaimal spectrum over the frequency domain is a constant, i.e.

8 ∞ X > > > J2n ðzÞcosð2nqÞ > cosðzsinqÞ ¼ J0 ðzÞ þ 2 < n¼1

þ∞ ð

∞ X > > > J2n1 ðzÞsin½ð2n  1Þq  > : sinðzsinqÞ ¼ 2

SKai ðz; uÞdu

n¼1

∞ þ∞ ð

¼2

 0

¼

(26)

50zu2*

pUðzÞ 1 þ

50 uz 2pUðzÞ

5=3 du

(23)

6u2*

which indicates that the turbulent kinetic energy described by the Kaimal spectrum does not vary with spatial positions. Therefore, the turbulent kinetic energy s2uðr;tÞ also maintains invariant after the disturbing of the rotational sampling.

4. Analytical expression of rotational sampling spectrum (RSS) To derive the analytical expression of the RSS, the rotational speed of the blade is assumed to be a constant. Then, the rotational sampling wind speed can be represented by Eq. (13).

where Jn ð,Þ is the n-th order of Bessel function of the first kind, n is a integer, and it has the following characteristics

Jn ðzÞ ¼ ð1Þn Jn ðzÞ;

cz

(27)

One can derive

i h ðzÞ ðyÞ ðzÞ ~ sinðu tÞ þ r ~ ~ j t þ fð1Þ kj f1;j ðtÞ ¼ cos r k kj cosðub tÞ þ u þ Hhub ~ b j j i h   ~ j t þ fjð1Þ þ 4j ¼ cos Bj sin ub t þ aj þ u i h   ~ j t þ fjð1Þ þ 4j ¼ cos Bj sin ub t þ aj cos u i   h ~ j t þ fjð1Þ þ 4j sin Bj sin ub t þ aj sin u ( ) ∞ i h X       ~ j t þ fjð1Þ þ 4j ¼ J0 Bj þ 2 Jn Bj cos n ub t þ aj cos u n¼even ) ( ∞ i h X     ~ j t þ fjð1Þ þ 4j Jn Bj sin n ub t þ aj sin u  2 n¼odd þ∞ X

¼

n¼∞ þ∞ X n¼∞

For simplicity, the Davenport spectrum SDav ðuÞ (double-sided spectrum) which is not dependent on height is employed instead of ð0Þ the Kaimal spectrum. The initial azimuth angle qi is set to be zero, and denotes ub ¼  2pfb . Thus, Eq. (13) is rewritten as

J 2n ðzÞ ¼ 1

n¼∞

¼ 4.1. Expression of the RSS

∞ X

i h     ~ j t þ fjð1Þ þ 4j Jn Bj cos n ub t þ aj þ u i h    ~ j t þ naj þ fjð1Þ þ 4j Jn Bj cos nub þ u

(28) 0 1 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðyÞ 2 2 ~ ~ðzÞ ~ðyÞ Þ þ ðk ~ðzÞ Þ ; a ¼ arctan@kj A; 4 ¼ H k where Bj ¼ r ðk ðzÞ j j hub j . j j ~ kj Similarly,

J. Chen et al. / Renewable Energy 146 (2020) 2170e2187

f2;j ðtÞ ¼ f3;j ðtÞ ¼ f4;j ðtÞ ¼

þ∞ X n¼∞ þ∞ X n¼∞ þ∞ X n¼∞

i h    ~ j t þ naj þ fjð2Þ þ 4j Jn Bj cos nub  u

FðuÞ ¼

i h    ~ j t  naj þ fjð3Þ þ 4j Jn Bj cos nub þ u

(29)

þ

i h    ~ j t  naj þ fjð4Þ þ 4j Jn Bj cos nub  u

2177

N n  h   i 1X ð0Þ  ~ j þ bjð0Þ d u  u ~j Aj  J0 Bj aj d u þ u 2 j¼1

þ∞  h X     i ð1Þ  ~ j þ bð1Þ ~ Jn Bj aj d u þ nub þ u d u  n u þ u j b j n¼1

þ

þ∞ X

   h ð2Þ    i ~ j þ bð2Þ ~j Jn Bj aj d u þ nub  u d u  n ub  u j

)

n¼1

Therefore, uðr; tÞ can be rewritten as

(33)

uðr; tÞ ¼

N X

(

n¼∞

j¼1

where Aj ¼ Since

1 2p

þ∞ ð

þ∞ X

Aj

i h   n  ~ j t þ naj þ fjð1Þ þ 4j Jn Bj  cos nub þ u i h  ~ j t þ naj þ fjð2Þ þ 4j þcos nub  u i h  ~ j t  naj þ fjð3Þ þ 4j þcos nub þ u ioo h  ~ j t  naj þ fjð4Þ þ 4j þcos nub  u

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~ðzÞ ; k ~ðyÞ ; u ~ Þ,VolðV Þ. 4Sðk j

j

j

~ j  0, and where u

j

i 1 h if e dðu þ u0 Þ þ eif dðu  u0 Þ 2

cosðu0 t þ fÞ,eiut dt ¼

∞

(30)

" ð0Þ aj

(31)





ð1Þ

i fj þ4j

¼ e





ð3Þ

i fj þ4j

þe



ð2Þ



i fj þ4j

þe



ð4Þ

i fj þ4j

#

þe

(34)

the Fourier spectrum of uðr; tÞ is derived as

1 FðuÞ ¼ 2p 1 ¼ 2

N X j¼1

(

þ∞ ð

uðr; tÞ,eiut dt

∞

Aj

þ∞ X n¼∞



#  ("   ð1Þ ð3Þ   i naj þfj þ4j i naj þfj þ4j þe Jn Bj d u þ n ub þ u e j "

ð1Þ



i naj þfj þ4j

þ e "



ð2Þ

þe

i naj þfj þ4j

þ e

"

ð2Þ

i naj þfj þ4j

þ e



ð3Þ

i naj þfj þ4j



ð4Þ

# 

i naj þfj þ4j

# 

þe

ð4Þ

i naj þfj þ4j

where dð,Þ is the Dirac function. Further, by adopting Eq. (27), FðuÞ can be written as





d u þ n ub  u j

þe





d u  nub þ u j

# 



d u  nub  u j

" ð0Þ bj



 ))

ð1Þ

i fj þ4j

¼ e

(32)





ð3Þ

i fj þ4j

þe





ð2Þ

i fj þ4j

þe





ð4Þ

i fj þ4j

#

þe

(35)

2178

J. Chen et al. / Renewable Energy 146 (2020) 2170e2187

" ð1Þ aj





ð1Þ

i naj þfj þ4j

¼ e

" þ ð1Þ

ð1Þ



i naj þfj þ4j

SR ðuÞ ¼ E jFðuÞj2 N n    X    ~j þ d u  u ~j A2j  J 20 Bj d u þ u ¼



ð4Þ

#

þe

ð3Þ

i naj þfj þ4j

(36)

j¼1



ð2Þ

i naj þfj þ4j



ð4Þ

i naj þfj þ4j

#

þe

e

þ∞ X

        ~ j þ d u  n ub þ u ~j J 2n Bj d u þ nub þ u n¼1 ) þ∞    X      ~ j þ d u  n ub  u ~j þ J 2n Bj d u þ nub  u þ

#

þe



" n

#

i naj þfj þ4j

e

¼ e

þ ð1Þ



ð2Þ

ð1Þ

ð3Þ

þe

i naj þfj þ4j

n

" bj



i naj þfj þ4j

n¼1

(40)

(37) Accordingly, the continuous form of the RSS is given by

þ∞ ð þ∞ ð ð þ∞

S R ð uÞ ¼

∞ 0

þ

þ∞ X

  n    1 ~  J 20 B kz ; ky ½dðu þ u ~ Þ þ dðu  u ~ Þ  4S kz ; ky ; u 2

0

   ~ ÞÞ þ dðu  ðnub þ u ~ ÞÞ J 2n B kz ; ky ½dðu þ ðnub þ u

n¼1

þ

)    ~ ÞÞ þ dðu  ðnub  u ~ ÞÞ dkz dky du ~ J 2n B kz ; ky ½dðu þ ðnub  u

þ∞ X n¼1

¼

þ∞ ð þ∞ ð

1 2

0

1 þ 2

þ∞ ð þ∞ ð

0

1 þ 2

ð2Þ

aj

þ∞  X       4 S kz ; ky ; u  nub þ S kz ; ky ; u  nub J 2n B kz ; ky dkz dky

þ∞  X       4 S kz ; ky ; u þ nub þ S kz ; ky ; u þ nub J 2n B kz ; ky dkz dky n¼1

0

"

0

n¼1

0

þ∞ ð þ∞ ð

0

i na

¼ e

ð2Þ j þfj þ4j

" þ ð1Þ

i na

þe



ð2Þ

"







#

ð3Þ

where



ð4Þ

i naj þfj þ4j





ð3Þ

þe

     ðkÞ 2  ðkÞ 2 Since Eðaj  Þ ¼ Eðbj  Þ ¼ 4;

r

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðkz Þ2 þ ðky Þ2 .Since

Sðkz ;

uÞ ¼

ky ;

i 1 h SR ðuÞ ¼ h0 SDav ðuÞ þ SDav ð  uÞ þ 2 þ∞ h i X 1 hn ðu þ nub Þ SDav ðu þ nub Þ þ SDav ð  u  nub Þ þ 2 n¼1

#

i naj þfj þ4j

ky Þ ¼

(38)

þe ð1Þ

Bðkz ;

SDav ðuÞrðWFÞ ðkz ;ky ; uÞ, SDav ðuÞ and rðWFÞ ðkz ; ky ; uÞ are symmetric, Eq. (41) can be further rewritten as

#

þe

i naj þfj þ4j

e

ð4Þ j þfj þ4j

i naj þfj þ4j

e

¼ e

þ ð1Þ



ð1Þ

i naj þfj þ4j

n



i naj þfj þ4j

n

" ð2Þ bj

(41)

       4 S kz ; ky ; u þ S kz ; ky ; u J 20 B kz ; ky dkz dky

#

h

þ∞ X 1

(39)

k ¼ 0; 1; 2, the RSS is finally

i

hn ðu  nub Þ SDav ðu  nub Þ þ SDav ð  u þ nub Þ

n¼1

2

¼ h0 SDav ðuÞ þ

þ∞ X

hn ðu þ nub ÞSDav ðu þ nub Þþ

n¼1

~ j  0), given by (u

þ∞ X

hn ðu  nub ÞSDav ðu  nub Þ

n¼1

(42) where u2ð∞; þ∞Þ and

J. Chen et al. / Renewable Energy 146 (2020) 2170e2187

h0 ðuÞ ¼

þ∞ ð ð þ∞

0

     4rðWFÞ kz ; ky ; u J 20 B kz ; ky dkz dky

(43)

þ∞ ð

0

SR ðuÞdu ¼

hn ðuÞ ¼

0

     4rðWFÞ kz ; ky ; u J 2n B kz ; ky dkz dky

(44)

" # þ∞   X   A2j  2J 20 Bj þ 4 J 2n Bj n¼1

j¼1

∞ þ∞ ð ð þ∞

N X

2179

¼

N X

"

A2j  2

n¼∞

j¼1

0

¼

þ∞ X

N X

# N X 2  J n Bj ¼ 2A2j j¼1

ðzÞ ðyÞ

  ~ ;k ~ ;u ~ j Vol Vj 8S k j j

j¼1 þ∞ ð þ∞ ð ð þ∞

¼

(45)

  S kz ; ky ; u dkz dky du

∞ ∞ ∞

4.2. Total energy of the RSS Based on Eq. (40), the total energy of the RSS over the frequency domain is derived as

in which Eq. (27) is adopted. Eq. (45) demonstrates that the kinetic energy of the fluctuating wind speeds keeps unchanged after the rotational sampling, which is consistent with the result in Sections 3.3.

Fig. 5. Energy transfer mechanism of the RSS.

2180

J. Chen et al. / Renewable Energy 146 (2020) 2170e2187

4.3. Energy transfer mechanism of the RSS The kinetic energy of fluctuating wind speeds over frequency domain exhibits a redistribution phenomenon due to the rotational sampling [37]. The energy transfer mechanism and its magnitude can be directly quantitively revealed from Eq. (40). For clarity, SR ðuÞ is divided into three parts, i.e., ð1Þ

ð2Þ

ð3Þ

SR ðuÞ ¼ SR ðuÞ þ SR ðuÞ þ SR ðuÞ

(46)

where

ð1Þ

SR ðuÞ ¼ ð2Þ SR ðuÞ ð3Þ

Another interesting phenomenon is that the values of the RSS at the exact multiples of the rotor’s rotational speed fall suddenly. This is because the value of the Davenport spectrum at zero frequency is zero. If the Kaimal spectrum is used, the result is different, which can be seen in Section 4.4. The comparison between Kaimal spectrum and Davenport spectrum is shown in Fig. 6. However, one should note that the zero frequency is often not included in the simulation of wind field, then the values of the RSS at the multiples of the rotor’s rotational speed will still fall suddenly, therefore, a smoothing technique is usually employed, just as seen in Fig. 9 and Fig. 10. To further illustrate the energy transfer mechanism of RSS, the

N n X       o ~j þ d u  u ~j A2j J 20 Bj d u þ u

j¼1 ( ) N þ∞    X X      ~ j þ d u  n ub þ u ~j ¼ J 2n Bj d u þ nub þ u A2j

SR ðuÞ ¼

j¼1 ( N X

n¼1 þ∞ X

j¼1

n¼1

A2j

(47)

)         ~ j þ d u  n ub  u ~j J 2n Bj d u þ nub  u

It is noted that J 20 ðBj Þ < 1 and the value of J 2n ðBj Þ decreases with the increasing of n. Therefore, one can observe that the energy at ~ j ðu ~ j  0Þ of the source spectrum, i.e., 2A2j , is the frequency u decomposed into three parts: (1) a ratio of J 20 ðBj Þ of the energy ~ j and  u ~ j ; (2) a ratio of J 2n ðBj Þ of the remains at the frequency u ~ j  nub , and ~ j þ nub and  u energy is transferred to the frequency u the transferred energy decreases with n; and (3) a ratio of J 2n ðBj Þ of ~ j þ n ub , ~ j  nub and  u the energy is transferred to the frequency u and the transferred energy also decreases with n. ð1Þ ð2Þ ð3Þ The physical meanings indicated in SR ðuÞ, SR ðuÞ and SR ðuÞ of Eq. (47) can be easily understood in view of Fig. 5. It is clearly revealed in Fig. 5 that the reason for the occurrence of peaks in the RSS around the multiples of the rotor’s rotational speed is that the vast majority of the energy of source spectrum is concentrated in the region close to zero frequency.

Fig. 6. Comparison between Kaimal spectrum and Davenport spectrum.

source wind spectrum is replaced by an assumed virtual spectrum Sðf Þ ¼ aebjf fp j , in which the size parameter a ¼ 30 m2 =s, the shape parameter b ¼ 5 s, and the position parameter fp ¼ 0:5 Hz. The rotor’s rotational speed fb ¼ 0:2 Hz. The comparison between the source spectrum and its corresponding RSS is shown in Fig. 7. In this case, there are no peaks around the multiples of the rotor’s rotational speed, and the energy is transferred from a narrow band to its both sides. Therefore, the influence of rotational sampling on wind spectrum is actually a spreading frequency effect. 4.4. Characteristics of the RSS In order to understand the characteristics of the RSS more comprehensively, the influences of the radial distance r and the rotational speed fb on the RSS are quantitatively investigated in this subsection. Consider a standard model of 5-MW offshore wind turbine with

Fig. 7. Comparison between an assumed spectrum and its RSS.

Speed(m/s)

J. Chen et al. / Renewable Energy 146 (2020) 2170e2187

5

r=0m

0 -5 0

Speed(m/s)

2181

100

200

300

400

5

500

600

r=30m

0 -5

Speed(m/s)

0

100

200

300

400

5

500

600

r=60m

0 -5 0

100

200

300

400

500

600

Time(s) Fig. 8. Rotational sampling fluctuating wind speeds at different positions of blade.

three blades [38]. The height of the hub is Hhub ¼ 90m, and the radius of the blades is R ¼ 63m. The cut-in and cut-out wind speeds are 3 m/s and 25 m/s, respectively. The cut-in and rated rotor’s speeds fb are 0.12 Hz and 0.20 Hz, respectively. The parameters of the joint spectrum representing the normal wind field, i.e. Eq. (3), take the following values: the mean wind speed at the hub height is Uhub ¼ 15m=s; the turbulence intensity is Iref ¼ 0:14 [2]; the standard deviation of turbulence is su ¼ Iref Uhub ¼ 2:1m=s. The shear velocity u* can be determined by Eq. (23): u* ¼ 0:857m=s; z0 and U10 are then specified by Eq. (4): z0 ¼ 0:082m, U10 ¼ 10:3m=s, and the decay coefficients Cz ¼ Cy ¼ 7 [39]. The other parameters for the numerical analysis are: the upU per cut-off wavenumbers and frequency are kU z ¼ ky ¼ p rad=m U and u ¼ 4p rad=s, the time history length of wind speeds is T ¼ 600s, and the time step of the simulation is Dt ¼ 0.25s. Case 1. Influences of radial distance r on RSS To investigate the influences of radial distance r on the RSS, the radial distance r is chosen to be 0 m, 30 m and 60 m, respectively, with a fixed rotational speed of blades fb ¼ 0:20Hz. A population of 100 samples of rotational sampling fluctuating wind speeds are simulated for each case based on Eq. (13), of which one typical sample at different positions of the blade is shown in Fig. 8. Based on the simulation samples, the RSS can be estimated by Ref. [27].

0 1 Su @uA ¼

E jFu ðuÞj2 T

(48)

where Fu ðuÞ denotes the finite Fourier transform of uðr; tÞ. The estimated RSS of rotational sampling fluctuating wind speeds at different positions of the blade is shown in Fig. 9. It is easily observed that the rotational sampling effect is more pronounced with the increase of r. To verify the conclusions of Sections 3.3 and 4.2, the variances of the rotational sampling fluctuating wind speeds at the three radial positions are calculated by statistics to be 4.20 m2/s2, 4.19 m2/s2 and

4.20 m2/s2 for r ¼ 0 m, r ¼ 30 m and r ¼ 60 m, respectively. Therefore, the total turbulent kinetic energy keeps unchanged but arises to redistribution over the frequency domain. Case 2. Influences of rotational speed fb on RSS To this end, the rational speed of blades fb takes the value of 0.1 Hz, 0.2 Hz and 0.3 Hz, respectively, with a fixed radial distance r ¼ 63m. Likewise, 100 samples of rotational sampling fluctuating wind speeds are simulated for each case. The estimated RSS for different rotational speeds is shown in Fig. 10. One might recognize that all the peaks occur at the multiples of fb , the turbulent kinetic energy is correspondingly transferred to a higher frequency region in the case of a larger fb. Therefore, the rotational sampling may have different influences on the fatigue loads of wind turbines with different rotational speeds.

5. Fatigue analysis of offshore wind turbine blades To investigate the rotational sampling effect on the fatigue lifetime of wind turbine blades, numerical analysis of a 5-MW offshore wind turbine is performed in this section with and without considering the rotational sampling, respectively. 5.1. Integrated numerical model for a 5-MW offshore wind turbine The stochastic dynamic response analysis of offshore wind turbine (StoDRAOWT) model [40,41] is employed here for the numerical analysis. This blade-hub-tower-pile integrated model is based on the definition of a 5-MW offshore wind turbine by the NREL [38]. The model is capable of considering the rotational effect of blades by combining Kane’s equation with the finite element method. With the Euler-Bernoulli beam model as the basic element, a total of 390 degrees of freedom are included to describe the wind turbine system. Besides, the geometric nonlinearity of blades, the aero-elastic effect, the blade-pitch control and the nonlinear generator torque control are all taken into account. Considering the soil-pile interaction, a nonlinear soil spring model

2182

J. Chen et al. / Renewable Energy 146 (2020) 2170e2187

101

0

PSD(m2 /s)

10-1

10-2

r=60m r=30m r=0 Kaimal Spectrum at Hub

3

2

PSD(m /s)

10

4

r=60m r=30m r=0 Kaimal Spectrum at Hub

10-3

2

1

0 10-1 Frequency(Hz)

100

0.2

(a) in logarithmic coordinates

0.4 0.6 Frequency(Hz)

0.8

1

(b) in normal coordinates

Fig. 9. Estimated RSS at different positions of blades.

101

4 f b=0.3Hz

10-1 f b=0.3Hz 10-2

f b=0.2Hz

3 PSD(m 2/s)

PSD(m 2/s)

100

f b=0.2Hz

f b=0.1Hz Kaimal Spectrum at Hub 2

1

f b=0.1Hz Kaimal Spectrum at Hub

10-3

0 10-1 Frequency(Hz)

100

0.1

(a) in logarithmic coordinates

0.2

0.3 0.4 0.5 Frequency(Hz)

0.6

0.7

(b) in normal coordinates

Fig. 10. Estimated RSS at the blade tip for different rotational speeds.

is employed. The StoDRAOWT model has been verified to be accurate and efficient in the structural dynamic response analysis of offshore wind turbines [42]. For details of the model, one can refer to Ref. [41].

The P-M spectrum has the following formulation [6].

5.2. Determination of environmental loads

in which, a ¼ 8:1  103 ; b ¼ 0:74 are the dimensionless constants, g denotes the gravitational acceleration and U19:5 denotes the mean wind speed at 19.5 m height. Then, the sea surface elevation hðtÞ can be simulated as

In the analysis of offshore wind turbine systems, the aerodynamic loads on the blades and tower and the wave loads on the monopile need to be taken into account. For the calculation of aerodynamic loads on the blades, the BEM method is employed as introduced in Section 3.2. The wind load per unit length applied on the tower is then given by Ref. [6].

"



4 #

ag2 g Sh ðuÞ ¼ exp  b 5 U 19:5 u juj

hðtÞ ¼

N qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X   2Sh ðuÞDucos uj t þ fj

(50)

(51)

j¼1

1 ft ðtÞ ¼ Cdt dt ra ½UðzÞ þ uðtÞ2 2

(49)

where Cdt and dt are the drag coefficient and diameter of the tower, respectively. It is noted that uðtÞ is directly generated by Eq. (8) since there is no rotational sampling involved in the wind field simulation of the tower. To calculate the wave loads, the Pierson-Moscowitz (P-M) spectrum is adopted to simulate the sea surface elevation firstly.

where uj ¼ jDu; j ¼ 1; 2; /; N; Du ¼ uu =N, uu denotes the upper cut-off frequency, N denotes the number of discretized frequencies, and fj denotes the random phase angle uniformly distributed over [0, 2p]. Further, based on the two-dimensional linear wave theory [6], the horizontal velocity and acceleration of a water particle at the position ðxw ; zw Þ can be expressed as

J. Chen et al. / Renewable Energy 146 (2020) 2170e2187

8

N > X >     uj cosh kj ðzw þ hÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > >   ðx ; z ; tÞ ¼ 2Sh uj Ducos kj xw  uj t þ fj v > x w w > > sinh k h < j j¼1

> N u2 cosh k ðz þ hÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > X     > j w j > a ðx ; z ; tÞ ¼  >   2Sh uj Dusin kj xw  uj t þ fj x w w > > h sinh k : j j¼1

in which xw denotes the horizontal coordinate, and is consistent with the mean wind direction in this paper; zw denotes the vertical coordinate, with the sea level as the original elevation and the upward direction as the positive direction; h denotes the water depth, which is chosen to be 18 m in the present study; kj denotes the wavenumber corresponding to the frequency uj and is determined by the following dispersion relationship [6]

u2 ¼ gk tanhðkhÞ

(53)

Thus, the wave load per unit length on the monopile can be calculated by Morison’s equation [6]

pdp 1 ax ðx; z; tÞ fp ðx; z; tÞ ¼ rw Cdp dp vx ðx; z; tÞjvx ðx; z; tÞj þ rw CI 2 4 2

2183

(52)

by

D¼

Nc X

1   N s m;i ; sa;i i¼1

(55)

where Nc denotes the number of full stress cycles experienced by a specified blade element under a certain wind condition within a certain time period. During the lifetime of a wind turbine, the wind conditions vary within a broad range. Since the wind conditions can be mainly characterized by the mean wind speed and the fluctuating wind speed for the wind turbine with yaw capability, the expected fatigue damage of a blade during its entire lifetime can be estimated by

(54) where Cdp , CI denote the drag and inertia coefficients of the pile, respectively, and Cdp ¼ 1:2, CI ¼ 2:0 are adopted in the present study; rw denotes the density of sea water; and dp denotes the diameter of the pile.

EðDÞ ¼

Tlife T

Uðout

EðDjUhub ; TÞpðUhub ÞdUhub Uin

¼

Tlife T

(56) NW   X    E DUj ; T P Uj j¼1

5.3. Prediction of the fatigue lifetime of wind turbine blades According to IEC code [2], the fatigue failure of wind turbine blades results from the accumulation of damage due to fluctuating loads. To estimate the fatigue damage, the number of full cycles of stress should be firstly counted by the Rainflow counting algorithm [43] after the structural dynamic analysis. The mean stress sm and the stress amplitude sa can be used to represent the characteristic of the stress cycle. Then, the fatigue damage increment caused by each stress cycle can be determined by means of a specified constant lifetime diagram (CLD). Further, the total fatigue damage can be obtained in terms of Miner’s rule with the assumption that the damage accumulates linearly and independently for each cycle. An optimal CLD for the fiberglass composite materials used in wind turbine blades was suggested by Sutherland & Mandell [44], and was further developed by Toft & Sørensen [3], which is adopted in this paper. To allow the usage of the CLD for lifetime prediction of the materials with similar fatigue behaviors, the CLD is usually presented in a normalized form [45]. In this case, both the axes of the CLD are normalized to the larger one between the static tensile and compressive strength. The normalized formulation of the CLD similar to Ref. [3] is presented in Fig. 11. In the CLD, N ¼ Nðsm ; sa Þ denotes the number of cycles to failure under the constant stress amplitude sa and the mean stress sm ; R ¼ ðsm  sa Þ=ðsm þ sa Þ, denotes the stress ratio. As seen in Fig. 11, the CLD is constructed by the log-log stress-cycle (S-N) curves at 7 R-values. The parameters for the S-N model were provided in previous investigations [3]. Thus, the number of cycles to fatigue for an arbitrary stress cycle can be obtained by interpolations in the CLD. The detailed interpolation procedure can be found in Ref. [3]. Then, the total damage D after a number of full stress cycles is given

where Tlife denotes the lifetime of the blade; Uin and Uout denote the cut-in and cut-out mean wind speeds of the wind turbine; pðUhub Þ denotes the probability density function of the mean wind speed at the hub; T is a specific time period, usually taken as the duration of one time of structural numerical analysis; NW denotes the number of selected representative mean wind speeds used for numerical analysis; PðUj Þ is the assigned probability of the j  th  representative mean wind speed Uj ; EðDUj ; TÞ denotes the expected damage under the mean wind speed Uj within the time period T, be estimated by Ref. [46] which involves the randomness of the fluctuating wind speeds and can NS      1 X E DUj ; T ¼ D Uj ; u l ; T NS

(57)

l¼1

where NS denotes the number of fluctuating wind field samples on the condition of Uj ; ul denotes the l  th fluctuating wind field sample; and DðUj ; ul ; TÞ can be calculated by Eq. (55). The mean wind speed at the hub Uhub is assumed to follow a Rayleigh distribution given by Ref. [2].

i h FðUhub Þ ¼ 1  exp  pðUhub =2Uave Þ2

(58)

where FðUhub Þ denotes the cumulative density function. The mean wind velocity is Uave ¼ 0:2Uref with the reference velocity defined as Uref ¼ 50m=s in wind turbine class I [2], which will be adopted in the numerical analysis. Thus, the lifetime can be estimated by setting EðDÞ ¼ 1 in Eq. (56).

2184

J. Chen et al. / Renewable Energy 146 (2020) 2170e2187

Normalised stress amplitude ( a /UTS)

0.50 N=10 3

R=-2.5

R=-1.0 R=-0.4

N=10 4 N=10 5

0.40

N=10 6 N=10 7

R=0.1 0.30

N=10 8

R=10.0

N=10 9 N=10 10

R=2.0

R=0.5

0.20

0.10

0.00 -1.0

-0.8

-0.6

-0.4

(UCS/UTS)

-0.2

0.0

0.2

Normalised mean stress (

0.4 m

0.6

0.8

1.0

(UTS/UTS)

/UTS)

Fig. 11. Normalized CLD for the materials of wind turbine blades. (UCS: ultimate compressive strength; UTS: ultimate tensile strength)

Further, the expected lifetime can be obtained from Eq. (59)

5.4. Numerical results The representative mean wind speeds Uj as well as their assigned probabilities are listed in Table 1. Under a given mean wind speed condition, a number of 30 deterministic structural analyses are performed. Each structural analysis is associated with a random fluctuating wind speed field. The duration of a deterministic structural analysis is T ¼ 500s. The stress of the blade root due to the moment in flapwise direction is considered since it is pertinent to fatigue lifetime prediction of blades [47]. The ultimate tensile strength of the blade is 255 MPa [48]. Using the above-mentioned procedure, the fatigue damage with and without considering the rotational sampling effect within the time period T is estimated and included in Table 1. For illustration, a representative time history of stress at the blade root in the case of Uhub ¼ 22m=s is shown in Fig. 12. The histograms of the full stress cycles identified from the stress time histories by the Rainflow counting algorithm are shown in Fig. 13. It is obvious that the stress amplitudes become larger with consideration of the rotational sampling effect. Within a time period of 500s, the number of full stress cycles are around 1200 and 200 in average, respectively, in the cases of considering and not considering the rotational sampling effect. Therefore, the estimated damage is larger when the rotational sampling effect is considered, as can be seen in Table 1.

Tlife ¼

N W P

T 1 ,T ¼      EðDjTÞ E DUj ; T P Uj

(59)

j¼1

Thus, the expected fatigue lifetimes with and without considering the rotational sampling effect are calculated as shown in Table 1. It is seen that the predicted fatigue lifetime of blades without consideration of rotational sampling effect increases by more than 50% compared to that with consideration of rotational sampling effect. Therefore, the rotational sampling of wind speeds has a significant influence on the fatigue lifetime of blades, and should be taken into account in the structural design of wind turbines. It is noted that the present analyses are carried out based on a variable speed pitch-regulated wind turbine model. However, the invariant speed wind turbines are widely used in early ages. Therefore, the influence of rotational speed on the fatigue life needs to be further investigated, in which the rotational sampling is considered, as shown in Table 1. It is clearly shown that the fatigue life of blades is very sensitive to the rotational speed and considerably decreases with increasing rotational speed of the rotor.

Table 1 Estimated fatigue damage and lifetime of the blade. Uj (m/s)

4

P(Uj) without Rotational Sampling with Rotational Sampling a fixed fb ¼ 0.19 Hz a fixed fb ¼ 0.20 Hz a fixed fb ¼ 0.21 Hz

0.1432 0.2216 E(D|T) ¼ 2.8747  107, Tlife ¼ 55 yrs E(D|T) ¼ 4.3839  107, Tlife ¼ 36 yrs Tlife ¼ 90 yrs Tlife ¼ 29 yrs Tlife ¼ 6 yrs

7

10

13

16

19

22

0.2131

0.1621

0.1016

0.0534

0.0295

J. Chen et al. / Renewable Energy 146 (2020) 2170e2187

2185

60 Stress at blade root (MPa)

With Rotational Sampling 40 20 0 -20 0

50

100

150

200

60

250

300

350

400

450

500

400

450

500

Without Rotational Sampling

40 20 0 -20 0

50

100

150

200

250

300

350

Time(s) Fig. 12. A representative time history of stress at the blade root.

5.5. Efficiency assessment To evaluate the efficiency of the proposed method in considering rotational sampling, the consumed time of the proposed method and the classical method is compared in this subsection. In the StoDRAOWT model, each blade is discretized into a number of 16 finite elements based on the properties of blades [38]. In the structural dynamic analysis, the aerodynamic load on each blade element is calculated based on the rotational sampling wind speed of that element. Therefore, only a total of 48 series of rotational sampling wind speeds need to be generated by means of the proposed method. On the other hand, in the traditional method a large series of normal fluctuating wind speeds should be firstly simulated at prescribed spatial positions. As shown in Fig. 14, the spatial range of wind field is denoted by [-64 m, 64 m]  [26 m, 154 m], and the coordinate separation is 8 m between two adjacent discretized positions [49]. Therefore, the wind speeds are firstly simulated at 289 spatial positions by the conventional Cholesky decomposition and FFT based method [22]. Then, the rotational sampling wind speeds of blades are obtained by the interpolations in spatial and temporal dimensions as the blades rotate. The duration and time

step of the simulated wind field are 600s and 0.25s, respectively, and the rotational sampling frequency takes 4 Hz. Routines are written in MATLAB platform and run on a PC with Intel (R) Core (TM) i7-4790 K CPU @ 4.00 GHz and 16 GB main memory under the WIN7 operating environment. Assuming that the blades rotate at a fixed speed fb ¼ 0.2 Hz, and only the rotational sampling fluctuating wind speeds for the 48 blade elements are simulated, the consumed time are 244s and 148s for the traditional method and the proposed method, respectively. Therefore, besides the merits of theoretical rigorousness and accuracy of avoiding interpolation, the proposed method is more competitive in efficiency as well. 6. Concluding remarks In conjunction with the wavenumber-frequency joint spectrum of fluctuating wind field, a novel method for rotational sampling of wind speeds has been proposed. Using this method, the instantaneous wind speeds experienced by wind turbine blades can be readily simulated as the rotor rotates. The analytical explicit expression of the rotational sampling spectrum has been derived and its characteristics have been investigated based on the

20

Number of cycles

Number of cycles

60

40

20

Str e

0 0

ss

20

plit

10 5 0 0

Str

10

am

15

ud

e(M

30

Pa

)

-20

0

20

40

(MPa)

Mean stress

(a) with rotational sampling

es

sa

10

mp

litu

de

20

(M

Pa

)

0

10

20

40

30

Pa)

tress(M

Mean s

(b) without rotational sampling

Fig. 13. Histograms of full stress cycles at the blade root within a time period of 500s.

2186

J. Chen et al. / Renewable Energy 146 (2020) 2170e2187

z

Hub

o

y

Fig. 14. Schematic diagram of wind field simulated by classical method.

proposed method. The fatigue lifetime of wind turbine blades has been estimated with and without considering the rotational sampling effect. The efficiency of the proposed method is addressed as well. Some concluding remarks are made as follows: (1) The proposed method provides an analytical expression of rotational sampling wind speeds, which circumvents the decomposition of the cross PSD matrix and the interpolation in spatial and temporal dimensions involved in the widelyused Sandia method. Besides, the simulation efficiency is acceptable for practical applications by integrating the stochastic harmonic function representation. (2) The rotational sampling spectrum (RSS) is given analytically, and the mechanism of transfer of turbulent kinetic energy in frequency domain is clearly revealed. The turbulent kinetic energy is redistributed in the RSS compared to the undisturbed wind spectrum. The RSS is significantly affected by both the rotational speed of the rotor and the radial distance from the hub center. (3) The rotational sampling effect significantly increases the fatigue load of wind turbine blades. Therefore, it is necessary to consider the rotational sampling in the prediction of fatigue lifetime of wind turbines, especially of the large-size wind turbines. (4) The proposed method in consideration of rotational sampling is more competitive in efficiency compared to the classical method. The wind field provided in the present paper is simulated on the basis of the empirical spectra. Future work may focus on extending the applicability of the approach to wind fields that are simulated on the basis of physical mechanism of atmospheric turbulence. Acknowledgments The supports of the National Natural Science Foundation of China (Grant Nos. 11672209, 51725804 and 51878505), the NSFCDFG joint project (Grant No. 11761131014) and the Committee of Science and Technology of Shanghai China (Grant No. 18160712800) are highly appreciated. References [1] Z.L. Zhang, A. Staino, B. Basu, S.R.K. Nielsen, Performance evaluation of fullscale tuned liquid dampers (TLDs) for vibration control of large wind turbines using real-time hybrid testing, Eng. Struct. 126 (2016) 417e431.

[2] IEC 61400-1, International Standard. Wind Turbines - Part 1: Design Requirements, third ed., 2005. Reference number IEC 61400-1: 2005(E). [3] H.S. Toft, J.D. Sørensen, Reliability-based design of wind turbine blades, Struct. Saf. 33 (6) (2011) 333e342. [4] J.W. Larsen, S.R.K. Nielsen, Nonlinear parametric instability of wind turbine wings, J. Sound Vib. 299 (2007) 64e82. [5] P.J. Murtagh, B. Basu, B.M. Broderick, Along-wind response of a wind turbine tower with blade coupling subjected to rotationally sampled wind loading, Eng. Struct. 27 (8) (2005) 1209e1219. [6] T. Burton, N. Jenkins, D. Sharpe, E. Bossanyi, Wind Energy Handbook, second ed., John Wiley & Sons, 2011. [7] J.R. Connell, Turbulence Spectrum Observed by a Fast-Rotating Wind Turbine Blade, Tech. Rep. PNL-3426, Pacific Northwest Laboratory, Richland, Washington, 1980. [8] D.C. Powell, J.R. Connell, Review of Wind Simulation Methods for HorizontalAxis Wind Turbine Analysis, Tech. Rep. PNL-5903, Pacific Northwest Laboratory, Richland, Washington, 1986. [9] R.M. Sundar, J.P. Sullivan, Performance of wind turbines in a turbulent atmosphere, Sol. Energy 31 (6) (1983) 567e575. [10] M. Shinozuka, Simulation of multivariate and multidimensional random processes, J. Acoust. Soc. Am. 49 (1) (1971) 357e368. [11] J.R. Connell, The Spectrum of Wind Speed Fluctuations Encountered by a Rotating Blade of a Wind Energy Conversion System, Tech. Rep. PNL-4083, Pacific Northwest Laboratory, Richland, Washington, 1981. [12] A. Spagnoli, L. Montanari, Along-wind simplified analysis of wind turbines through a coupled blade-tower model, Wind Struct. 17 (6) (2013) 589e608. [13] G.L. He, J. Li, Wind field simulation of wind turbine systems, J. Tongji Univ. Nat. Sci. 38 (7) (2010) 976e981 ([in Chinese]). [14] W. He, D. Tian, W.L. Wang, Effect of wind shear on rotational Fourier spectrum of wind turbine, Appl. Mech. Mater. 271 (2013) 872e876. [15] W. Bierbooms, J.B. Dragt, A stochastic 3D wind field generator for design calculations, in: Proceedings 1996 European Wind Energy Conference, Goteborg, 1996, pp. 942e945. [16] P.S. Veers, Modeling Stochastic Wind Loads on Vertical Axis Wind Turbines, Tech. Rep. SAND83-1909, Sandia National Laboratories, Albuquerque, New Mexico, 1984. [17] P.S. Veers, Three-Dimensional Wind Simulation, Tech. Rep. SAND88-0152, Sandia National Laboratories, Albuquerque, New Mexico, 1988. [18] P.J. Moriarty, A.C. Hansen, AeroDyn Theory Manual, Tech. Rep. NREL/TP-50036881, National Renewable Energy Laboratory, Golden, Colorado, 2005. [19] B.J. Jonkman, L. Kilcher, TurbSim User’s Guide: Version 1.06.00, Tech. Rep. NREL/TP-500-39797, National Renewable Energy Laboratory, Golden, Colorado, 2006. [20] W.F. Tao, B. Basu, J. Li, Reliability analysis of active tendon-controlled wind turbines by a computationally efficient wavelet-based probability density evolution method, Struct. Control Health Monit. 25 (2018). https://doi.org/10. 1002/stc.2078. [21] K.A. Abhinav, N. Saha, Nonlinear dynamical behaviour of jacket supported offshore wind turbines in loose sand, Mar. Struct. 57 (2018) 133e151. [22] M. Shinozuka, C.M. Jan, Digital simulation of random processes and its applications, J. Sound Vib. 25 (1) (1972) 111e128. [23] Y.L. Li, K. Togbenou, H.Y. Xiang, N. Chen, Simulation of non-stationary wind velocity field on bridges based on Taylor series, J. Wind Eng. Ind. Aerodyn. 169 (2017) 117e127. [24] T.Y. Tao, H. Wang, A. Kareem, Reduced-Hermite bifold-interpolation assisted schemes for the simulation of random wind field, Probab. Eng. Mech. 53 (2018) 126e142. [25] S.T. Ke, T.G. Wang, Y.J. Ge, H. Wang, Wind-induced fatigue of large HAWT coupled towereblade structures considering aeroelastic and yaw effects, Struct. Design Tall Spec. Build. 27 (2018). https://doi.org/10.1002/tal.1467. [26] H.R. Zuo, K.M. Bi, H. Hao, Dynamic analyses of operating offshore wind turbines including soil-structure interaction, Eng. Struct. 157 (2018) 42e62. [27] J.B. Chen, Y.P. Song, Y.B. Peng, P.D. Spanos, Simulation of homogeneous fluctuating wind field in two spatial dimensions via a wavenumber-frequency joint power spectrum, J. Eng. Mech. 144 (11) (2018), 04018100. [28] Y.P. Song, J.B. Chen, Y.B. Peng, P.D. Spanos, J. Li, Simulation of nonhomogeneous fluctuating wind speed field in two-spatial dimensions via an evolutionary wavenumber-frequency joint power spectrum, J. Wind Eng. Ind. Aerodyn. 179 (2018) 250e259. [29] Y.P. Song, J.B. Chen, M. Beer, L. Comerford, Wind speed field simulation via stochastic harmonic function representation based on wavenumberfrequency spectrum, J. Eng. Mech. (2019) (accepted). [30] J. Mann, The spatial structure of neutral atmospheric surface-layer turbulence, J. Fluid Mech. 273 (1994) 141e168. [31] M. Harte, B. Basu, S.R.K. Nielsen, Dynamic analysis of wind turbines including soil-structure interaction, Eng. Struct. 45 (2012) 509e518. [32] E. Simiu, R.H. Scanlan, Wind Effects on Structures, third ed., John Wiley & Sons, 1996. [33] J.B. Chen, J.R. He, X.D. Ren, J. Li, Stochastic harmonic function representation of random fields for material properties of structures, J. Eng. Mech. 144 (7) (2018), 040118049. [34] J.B. Chen, W.L. Sun, J. Li, J. Xu, Stochastic harmonic function representation of stochastic processes, J. Appl. Mech. 80 (1) (2013), 011001. [35] J.B. Chen, F. Kong, Y.B. Peng, A stochastic harmonic function representation for non-stationary stochastic processes, Mech. Syst. Signal Process. 96 (2017)

J. Chen et al. / Renewable Energy 146 (2020) 2170e2187 31e44. [36] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover Publications, 1983. [37] J.B. Dragt, The spectra of wind speed fluctuations met by a rotating blade, and resulting load fluctuations, in: Proceedings 1984 European Wind Energy Conference, Hamburg, 1984, pp. 453e458. [38] J. Jonkman, S. Butterfield, W. Musial, G. Scott, Definition of a 5-MW Reference Wind Turbine for Offshore System Development, Tech. Rep. NREL/TP-50038060, National Renewable Energy Laboratory, Golden, Colorado, 2009. [39] Y.B. Peng, S.F. Wang, J. Li, Field measurement and investigation of spatial coherence for near-ground strong winds in southeast China, J. Wind Eng. Ind. Aerodyn. 172 (2018) 423e440. [40] J.B. Chen, Y.K. Liu, X.Y. Bai, Shaking table test and numerical analysis of offshore wind turbine tower systems controlled by TLCD, Earthq. Eng. Eng. Vib. 14 (1) (2015) 55e75. [41] J.B. Chen, T. Sun, K. Huang, J. Li, Study on integrated numerical modeling of offshore wind turbine tower systems, J. Dyn. Control 15 (3) (2017) 268e278 ([in Chinese]). [42] T. Sun, J.B. Chen, Comparison of the integrated numerical model for wind turbine systems with FAST, in: Symposium on Reliability of Engineering System, Hangzhou, 2015.

2187

[43] S.D. Downing, D.F. Socie, Simple Rainflow counting algorithms, Int. J. Fatigue 4 (1) (1982) 31e40. [44] H.J. Sutherland, J.F. Mandell, Optimized constant life diagram for the analysis of fiberglass composites used in wind turbine blades, J. Sol. Energy Eng. 127 (4) (2005) 563e569. [45] D.U. Shah, P.J. Schubel, M.J. Clifford, P. Licence, Fatigue life evaluation of aligned plant fibre composites through S-N curves and constant-life diagrams, Compos. Sci. Technol. 74 (2013) 139e149. [46] R.R. Pedersen, S.R.K. Nielsen, P. Thoft-Christensen, Stochastic analysis of the influence of tower shadow on fatigue life of wind turbine blade, Struct. Saf. 35 (1) (2012) 63e71. [47] P.J. Moriarty, W.E. Holley, S.P. Butterfield, Extrapolation of Extreme and Fatigue Loads Using Probabilistic Methods, Tech. Rep. NREL/TP-500-34421, National Renewable Energy Laboratory, Golden, Colorado, 2004. [48] W.Y. Tang, W.G. Lyu, D.Y. Li, Fatigue life prediction of wind turbine blade considering effect of mean stress, Acta Energiae Solaris Sin. 37 (10) (2016) 2703e2709 ([in Chinese]). [49] Y. Chen, J.Y. Zhang, N. Wang, Z.Q. Ye, Wind turbine wind field models study and numerical simulation of turbulence wind field with MATLAB, Acta Energiae Solaris Sin. 27 (9) (2006) 955e960 ([in Chinese]).