Dynamic analysis of offshore wind turbine blades under the action of wind shear and fluctuating wind

2.1. Physical model under FSI

In order to effectively avoid the deficiency that wind tunnel test can only be based on scale model, and really reflect the interaction between the blade structure of wind turbine and the surrounding flow field, a FSI physical model at the same proportion with the actual blades was established in this work. When considering the rotating state of blades, a rotating internal flow field was set in the blade periphery. Meanwhile, due to the influence of upstream wind turbine on the downstream one, the length of inflow region in external flow field is twice that of the blade, and the length of wake region is six times that of the blade, as shown in Fig. 1.

Fig. 1Physical model under FSI

2.2. Governing equations of fluid domain and structure domain

In the numerical calculation of FSI, Lagrange coordinate system focuses on the motion of material points in solid mechanics and Euler coordinate system concentrates on the state of space points in fluid mechanics. In order to solve the disunity problem of coordinate system at motion interface between fluid domain and solid domain, $\varphi$ is represented as the universal variable, which stands for the velocity components $u^{*}$, $v^{*}$, $w^{*}$ in three-dimensional coordinates, so the fluid governing equations built by ALE method can be written as [12]:

where $\rho$ is the air density. For the continuity equation, $\varphi =$ 1. For the momentum equation, $\varphi =u_j^{*}$, where $u_i^{*}$ and $u_j^{*}$ stand for speed components in each direction, and the indexes of $i$, $j$ range from 1 to 3. For the two-equation $k$-$\epsilon$ turbulence model, $\varphi =\left\{k,\epsilon \right\}$. $\mathrm{\Gamma }$ and $S$, which are ascertained according to $\varphi$, are respectively the generalized diffusion coefficient and generalized source term.

With FSI and wind shear effect taken into account, Eq. (1) is solved by iteration computation to obtain the velocity field which satisfies with the continuity equation, so that the wind pressure distribution on the blade surface is calculated.

In consideration of FSI, the blade motion equations with geometrical nonlinearity can be expressed as [13]:

where ${\sigma }_{ij}$, $U_i$, $\stackrel{-}{\rho }$ are the stress tensor, displacement vector and material density respectively, and $f_i$ is a vector related to the wind pressure distribution.

The constitutive equation and geometric equation for a geometrically nonlinear elastomer can be respectively written as:

where $G$, $\lambda$ are respectively represented as the shear modulus and Lame coefficient, and ${\delta }_{ij}$ is the unit tensor.

2.3. Calculation of wind shear and fluctuating wind speed

Wind shear refers to the wind variation rule and characteristics in the vertical direction in the field of wind power application. Generally, the exponential distribution of experience is used to describe the changing law of the average wind speed near surface layer with the height, and distribution function is defined as [14]:

where $\stackrel{-}{V}\left(z\right)$ is the wind speed at a vertical height $z$ above ground, and 120 m $\le z\le$ 183 m in this work is obtained by assuming that the height from ground to wheel hub is 120 m and the blade length is 63 m. ${\stackrel{-}{V}}_{z_0}$ is the average wind speed at the standard reference height $z_0$, and $z_0$ is generally set as 10 m. $h$ is the empirical shear coefficient associated with atmospheric stability and surface roughness, with a general value of 0.12 for offshore seas, islands, coasts, lakeshore and desert regions [15].

The instantaneous wind speed replaced by random fluctuating one is employed to simulate atmospheric turbulence motion. The wind speed $V$($x$, $y$, $z$, $t$) at any point ($x$, $y$, $z$) on the structure can be expressed by the sum of fluctuating wind speed $v$($x$, $y$, $z$, $t$) and average one $\stackrel{-}{V}\left(z\right)$, namely:

Davenport spectrum adopted widely in engineering is applied to reveal the characteristics of fluctuating wind, and the expressions of power spectrum density are given as [16]:

where $S_v\left(n\right)$ is the power spectrum density of fluctuating wind at the pulsating frequency $n$, ${\stackrel{-}{V}}_{10}$ is the average wind speed at the standard reference height of 10 m, and $k$ is the ground roughness coefficient which is selected as 0.0042 according to empirical shear coefficient [17].

The speed Fourier spectrum is obtained by using the speed power spectrum and the coherent function of a point, and then the values of fluctuating wind speed are eventually computed by inverse transformation. According to the calculation formula of Davenport spectrum, fluctuating wind speed is simulated by the programming of user defined function (UDF) in MATLAB software. When the frequencies of fluctuating wind with ${\stackrel{-}{V}}_{10}=$ 11.4 m/s and ${\stackrel{-}{V}}_{10}=$ 19.8 m/s are taken as 0.01 Hz and 0.015 Hz respectively, the time-history responses of instantaneous wind speed $v\left(t\right)$ are shown in Fig. 2.

Fig. 2Time-history responses of instantaneous wind speed

a)${\stackrel{-}{V}}_{10}=$ 11.4 m/s

b)${\stackrel{-}{V}}_{10}=$ 19.8 m/s

3.1. Entity modeling and parameter setting

Considering the wind field condition of large-scale offshore wind turbines, external flow field model at the same proportion with the actual wind field was established in this work. The three-dimensional model of the blade and fluid domain was built by using UG and Geometry module, as displayed in Fig. 3.

3.2. Grid division and independent validation

Grid generation of entity model is the process of discretization of the solution area, and lays the foundation for numerical solution. The shape, size and method of grid generation have a direct impact on the grid quality, thus affecting the accuracy of the numerical calculation process. Based on the actual computational requirements, the tetrahedral mesh is adopted, and the finite element grids of fluid domain were obtained by dynamic mesh technique with ICEM software, as shown in Fig. 4.

Fig. 3Three-dimensional model of fluid domain and structure domain

Fig. 4Grid generation of fluid domain

a) External flow field

b) Internal flow field

Here to validate the mesh independence of 3D FSI model under rated conditions, the results are listed in Table 1. It can be seen that accurate results can be obtained by using 862496 meshes in fluid domain and 61522 meshes in structure domain, thereby this kind of mesh division is used here to save computational resources.

5.1. Wind shear effect

For the blade under the rated rotating speed and the average wind speed of 11.4 m/s or 19.8 m/s, the comparisons of displacement response and stress response with the wind shear considered or not are respectively carried out, and the corresponding four kinds of working conditions are recorded as Condition I, II, III and IV. Fig. 6 depicts comparative response curves of maximum displacement and Mises stress along the blade wingspan under four conditions. The following can be easily concluded.

1) Response curves of displacement and stress under four conditions all present a nonlinearly attenuation trend. The decay rates of the response curves of Condition I and Condition II are similar, and Condition III and Condition IV are also the case, indicating that the influence of wind shear on the attenuation of the vibration curve is small.

2) From Condition I to IV, the maximum displacement and stress occur at the time of 2.0 s, 1.9 s, 1.8 s and 1.7 s respectively, and the corresponding maximum values of displacement are respectively 1.74 m, 2.00 m, 2.19 m and 2.57 m, while the ones of stress are respectively 16.11 MPa, 18.56 MPa, 20.32 MPa and 23.75 MPa, showing that wind shear effect causes the maximum values of displacement and stress to advance by 0.1 s and the increase of average wind speed makes them occur ahead of 0.2 s.

3) The maximum displacement and Mises stress under Condition II increase by 14.94 % and 15.21 % respectively compared with those under Condition I. Meanwhile, the maximum displacement and Mises stress under Condition IV increase by 17.35 % and 16.88 % respectively compared with those under Condition III. The two sets of data both show that wind shear contributes a lot to the vibration characteristics of wind turbine blades, so wind shear effect should be paid more attention to during the design of large-scale offshore wind turbines.

4) By comparing the first and second sets of data in 3), it can be seen that as the wind speed increases from 11.4 m/s to 19.8 m/s, the contribution of wind shear to maximum displacement increases from 14.94 % to 17.35 %, and the contribution to maximum Mises stress increases from 15.21 % to 16.88 %, implying that the wind shear effect becomes more remarkable with increasing average wind speed.

At the corresponding moments when the displacement and Mises stress under four conditions reach the maximum respectively, the maximum displacement and Mises stress of the blade section along the wingspan are displayed in Fig. 7. It is clear that:

1) The maximum displacement along the wingspan increases nonlinearly and reaches the maximum at the blade tip, and the maximum Mises stress at the relative wingspan of 0.55 runs up to the maximum and reduces basically towards each end of the blade.

2) On account of the existence of stress concentration on the joint between wheel hub and blade, the extreme stress value occurs in the blade root. From the relative wingspan of 0.15 to 0.7, the blade stress is comparatively large and multiple stress extremes appear, owing to the fact that DU airfoils with different thickness result in the change of flexural section coefficients. While from the relative wingspan of 0.7 to 1.0, Mises stress along wingspan decreases smoothly, which is attributed to using airfoil named NACA64_A17 in this region.

3) After comparing Condition I with II and comparing Condition III with IV, it is identified that displacement difference increases in the direction of the wingspan and reaches the maximum at the blade tip, while stress difference is up to the maximum near the blade central section and basically decreases along both sides of the blade, meaning that the contribution of wind shear to displacement is further greater along wingspan and the contribution to stress from the middle to each end of the blade generally decreases gradually. Furthermore, both the displacement difference and the stress difference with wind shear effect or not increase with wind speed, which also suggests that wind shear effect is more apparent when the average wind speed is larger.

4) By comparing Condition I with III and comparing Condition II with IV, it is found that whether wind shear effect is considered or not, the displacement difference and the stress difference both present the same rule with those described in (3), indicating that the contribution of wind speed increment to displacement increases along wingspan and the contribution to stress reduces from the blade middle to each end. In addition, the displacement difference and the stress difference with wind shear are both greater than those without wind shear, showing that the increase of average wind speed in consideration of wind shear effect contributes more to displacement and stress.

Fig. 6Comparative curves under four conditions

a) Displacement response

b) Mises stress response

Fig. 7Changing curves with relative span length

a) Maximum displacement distribution

b) Maximum Mises stress distribution

5.2. The combined action of wind shear and fluctuating wind

Aiming at the rated rotating speed, response curves of displacement and Mises stress for the blade under average wind speed of 11.4 m/s with the wind shear alone, the fluctuating wind alone or the combined action and under average wind speed of 19.8 m/s with the combined action are respectively plotted in Fig. 8, and the corresponding four working conditions are renamed Condition I, II, III and IV. It is not difficult to see that:

1) Response curves of displacement and Mises stress under four conditions all show a fluctuating attenuation tendency and the decay rates are nearly the same, implying that the effects of wind shear, fluctuating wind or the combined action on vibration attenuation are almost equal.

2) From Condition I to IV, the moments when displacement and stress are up to the maximum are 1.9 s, 2.0 s, 1.9 s and 1.7 s respectively, the maximum values of displacement are 2.00 m, 1.74 m, 1.99 m and 2.52 m respectively, and the ones of stress are respectively 18.56 MPa, 16.24 MPa, 18.60 MPa and 23.55 MPa. In contrast to Fig. 6, it can be observed that the combined action makes the maximum displacement and stress advance by 0.1s compared with those under only average wind speed effect.

3) The maximum displacement and Mises stress under Condition I increase by 13.00 % and 12.50 % by comparing with those under Condition II, which indicates the effect of wind shear on vibration characteristic is larger than that of considering fluctuating wind separately under the same average wind speed. Besides, the maximum displacement and Mises stress under Condition IV increase by 26.63 % and 26.61 % compared with Condition III, which implies that the increase of average wind speed makes the maximum displacement and Mises stress under the combined action increase comparatively obviously.

4) The maximum displacement, the maximum Mises stress and the response curves during the whole time history under Condition I are very close to those under Condition III, which reveals that the fluctuating wind included in the combined condition has a quite limited effect on the maximum value of displacement and Mises stress. The maximum displacement and Mises stress under Condition III increase by 12.56 % and 12.69 % compared with Condition II, meaning the wind shear included in the combined condition has a relatively evident effect on them.

Fig. 8Comparative curves of the maximum value

a) Displacement response

b) Mises stress response

At the moments corresponding to the maximum values of the response curves under four conditions, the maximum displacement and Mises stress in blade section along the wingspan are depicted in Fig. 9, and it can be easily found that:

1) Distribution trends of the maximum displacement and Mises stress in blade section along the wingspan under Conditions III and IV do not change greatly compared with those under Conditions I and II, that is to say, the maximum displacement in the blade section along wingspan increases nonlinearly and reaches the maximum at the tip, and the maximum Mises stress occurs at the relative wingspan of 0.55 and reduces nonlinearly along both sides of the blade, indicating that the combined action has no effect on the distribution trends of displacement and stress.

2) The displacement difference between Condition I and Condition II increases along the wingspan and the corresponding stress difference decreases from the middle to each blade end, showing the effect of wind shear in the whole span is greater than that of fluctuating wind alone. Moreover, by comparing Condition IV with Condition III, it is not difficult to find that the contribution of wind speed increment to the displacement under the combined action increases along the wingspan, while the contribution to the stress decreases gradually from the blade middle to both sides.

3) The distribution curves of the maximum displacement and Mises stress in the blade section along the wingspan under Condition I are almost coincided with those under Condition III, implying that the influence of fluctuating wind in the combined action on the distribution of displacement and stress along the wingspan is very small. Furthermore, by comparing Condition II with Condition III, the displacement difference increases along the wingspan and the stress difference decreases nonlinearly from the middle to each blade end, which reveals that the effect of wind shear in the combined action on the displacement distribution and Mises stress distribution is more obvious.

Fig. 9Variation curves along the wingspan

a) Displacement distribution

b) Mises stress distribution

The maximum Mises stress locates at the relative wingspan of 0.55, so the corresponding blade section was selected to carry out the comparative analysis for Mises stress distribution under two kinds of average wind speeds. Under Conditions III and IV, distribution cloud pictures at the moment when Mises stress is up to the maximum are given in Fig. 10, and it can be obtained that:

1) In cross section, one can set the part from the leading edge to the centre as the front half section, and the rest is the back half section. Figs. 10(a) and (b) both show that Mises stress of the front half section presents a gradient distribution in blade thickness direction and reaches the maximum on the blade surface, while the one of the back half section shows a gradient distribution in the blade chord direction and achieves the maximum near the center of chord length, indicating that the front half section is mainly influenced by the wind load in the flapping direction and the back half section is mainly subjected to the wind load in the torsion direction.

2) In the maximum stress section, the corresponding distribution trends of Mises stress under the two wind speeds are nearly accordant, and the maximum value increases with wind speed.

Fig. 10Distribution cloud pictures of Mises stress in the blade cross section

a)${\stackrel{-}{V}}_{10}=$ 11.4 m/s

b)${\stackrel{-}{V}}_{10}=$ 19.8 m/s