Dynamic stability analysis on large wind turbine blade under complicated offshore wind conditions

2.1. Random wind speed simulation

According to more than 90 strong wind records measured from different locations and heights in the world, Davenport assumed that the turbulence integral scale of the horizontal gust spectrum is unchanged along the height and average measured wind speeds of different heights above ground. Therefore he put forward horizontal pulsating wind speed spectrum, which can be expressed as [11]:

where $S_v\left({\omega }_0\right)$ is the pulsating power spectrum, ${\omega }_0$is the pulsating wind circular frequency, ${\stackrel{-}{V}}_{10}$ is the wind speed in ten meters standard height, $k$ is the ground roughness coefficient. According to the random vibration theory and the China load specification, the coefficient $k$ is written as [12]:

where ${\alpha }_0=0.12$ is the wind shear index obtained from experience.

Matlab, which is an interactive advanced scientific computing software produced by MathWorks, is employed to simulate the wind speed spectrum, and 20 m/s average wind speed curve with time is plotted in Fig. 1.

Fig. 1Wind speed time curve at 20 m/s

2.2. Blade aerodynamic loads

Based on momentum blade element theory used to solve the blade aerodynamic loads, the blade is divided into a number of infinitesimal segments called blade element. Suppose the flow between elements does not interfere with each other, and then the blade element is regarded as a two-dimensional blade profile. The blade force can be figured out after the force on blade element is integrated along wingspan direction. For each blade element, the velocity $V_0$ can be separated into $V_{0x}$ and $V_{0y}$, which are components perpendicular and parallel to rotating surface of the wind wheel, respectively. In detail, the force analysis is shown in Fig. 2.

Fig. 2Airflow speed triangle and aerodynamic component on blade element

According to momentum theory, aerodynamic force $d{F}_a$ on blade element with a $dr$ length caused by synthesis airflow velocity $V_0$, can be divided into normal force $d{F}_n$ and tangential force $d{F}_t$, and the expressions are given by [13]:

where $\rho$ is the air density, $c$ is the section chord length of blade element, $C_l$ and $C_d$ are the lift coefficient and the drag coefficient of blade element airfoil respectively, $\phi$ is the inflow angle of blade element, $\phi =\alpha +\theta$, $\alpha$ and $\theta$ are attack angle and geometric twist angle respectively. The synthesis airflow velocity $V_0$ of blade element can be expressed as:

where $a=v_a/V_1$, $b={\omega }_z/\left(2\omega \right)$, $v_a$ is the axial induced speed of wind wheel, $V_1$ is inflow wind speed, ${\omega }_z$ is the circumferential induced angular velocity at the diameter $r$ position of wind wheel, $\omega$ is rotational speed of wind wheel.

2.3. Structural motion equation

Differential motion equation of the blade is:

where $\left[M\right]$, $\left[C\right]$ and $\left[K\right]$ are mass matrix, damping matrix, stiffness matrix respectively, $\left[F\right]$ is blade load vector which can be obtained from formulas (4) and (5). Neglecting the influence of damping on the blade and introducing Newmark method into Eq. (7), one can get [14, 15]:

where $\alpha$ is the parameter determined by integral accuracy and stability.

Fig. 3Blade solid model and finite element meshing model

a) Blade solid model

b) Finite element model of entire blade

c) Enlarged model

4.1. Numerical simulation process

Numerical simulation includes wind load calculation and dynamic analysis. Matlab and Profili (an aircraft airfoil design and aerodynamic analysis software) are employed to find the wind speed and dynamic parameters in the process of wind load calculation. Finite element analysis is carried out by using ANSYS. The detail of calculation process is shown in Fig. 4.

Fig. 4Numerical simulation flowchart

In the numerical simulation on the dynamical response of the blade, as described by above flowchart, the following calculation procedure is performed.

Step 1: Establish the theoretical model of Davenport wind spectrum, and build the wind field at different working conditions through simulation of random wind speed using Matlab.

Step 2: Obtain detailed data of wind turbine blade elements, such as chord length, attack angle, airfoil, and so on; create table files including attack angle $\alpha$ in direction of blade wingspan and chord length $c$.

Step 3: Calculate the Reynolds number for blade elements along the wingspan direction at different wind speeds by using Profili, and then find the lift coefficient $C_l$ and the drag coefficient $C_d$ for blade elements under the corresponding working conditions.

Step 4: Calculate the normal force and the tangential force of the blade by substituting the data obtained from above steps into Eqs. (4) and (5). Then, the solution of the blade aerodynamic load can be achieved.

Step 5: Establish the solid model of the blade using UG. Mesh the solid model and define constraints by using ANSYS.

Step 6: Load the aerodynamic forces obtained from Step 4 onto the stress surface of the blade, which has been meshed by employing table loading function in ANSYS. Set a reasonable load step and solve the deflection and the stress. Carry out the post-processing and analysis.

4.2. Reliability verification of simulation results

In order to verify the reliability of numerical results, the testing data in reference [17], where the wind turbine model is FD2-300 with a total power of 300 W, an airfoil of NACA632 series and a blade length of 0.85 meters, were compared with the data obtained from simulation method in this work. Aerodynamic tangential component, centrifugal force of rotating blade and blade weight were ignored in the test, and the aerodynamic force along blade wingspan direction was simplified into 7 concentrated forces distributed uniformly, as shown in Fig. 5. Combining geometric parameters with boundary conditions, the boot stress distribution of blade was found and is shown in Fig. 6. The maximum stress of the blade root under simulation is 41.4 MPa, which is very close to the test data value of 42.7 MPa with a relative error about 3.04 %. The result confirms the reliability of calculation program, the theoretical and numerical models used in this work.

Fig. 5Loading diagram of concentrated force on blade

Fig. 6The maximum stress cloud picture of blade root

5.1. Vibration and stress analysis under periodic wind speeds

Fig. 7 shows the time spectrum of waving vibration of blade tip and middle part when the average speed of the complicated offshore wind is up to 20 m/s, which is close to the first natural frequency. It is evident from Fig. 7 that the waving vibrations diverge and the amplitudes increases gradually from the blade boot to the blade tip.

It was found from Fig. 7(b) and Fig. 8 that the maximum Mises stress increases with the increase of the waving amplitude, also showing divergence phenomena.

Fig. 7Waving vibration characteristics of different blade parts at periodic wind speed

a) The tip part of blade

b) The middle part of blade

Fig. 8The divergence characteristics of the maximum Mises stress at periodic wind speed

5.2. Determination of critical curve

In applied wind farms, there are many kinds of wind speeds which may cause blade vibration to diverge. Therefore, it is very important to study variation laws of the amplitude and the stress at different wind speeds when the blade is unstable. The outcomes from the study can provide guidance for strength optimization design of the wind turbine blade.

Fig. 9Characteristic curves of maximum Mises stress and maximum waving amplitude vs. average wind speed with different elastic modulus

a) Maximum waving amplitude vs. average wind speed

b) Maximum Mises stress vs. average wind speed

Fig. 9 gives the relationships among maximum Mises stress, maximum waving amplitude and average wind speed with different elastic modulus. Both the maximum Mises stress and the maximum waving amplitude of the blade increase nonlinearly with the average wind speed increasing. In addition, the maximum waving amplitude decreases and the maximum Mises stress increases, respectively, when elastic modulus is increased.

In view of the practical situation where the blade waving amplitude can not exceed the horizontal distance between the blade and the tower, the blade maximum waving amplitude is set at 1.5 meters according to a blade length of 40 meters, and it is then reduced to 1.0 meter for a safer design. This implies that the wind turbines cannot operate safely when the waving amplitude is more than above two set amplitudes. Under this condition, the quantitative relationship between maximum Mises stress and elastic modulus was studied.

Fig. 10Characteristic curves of the critical value of maximum Mises stress vs. elastic modulus with different horizontal distances

Fig. 10 plots the characteristic curves of the critical value of maximum Mises stress as a function of elastic modulus for different horizontal distances between the blade and the tower when the vibration diverges. Fig. 10 exhibits that the maximum Mises stress increases almost linearly with elastic modulus, and this linear trend does not vary with the horizontal distance. The dynamic stability region of blade safe operation is displayed in Fig. 10, which indicates that the elastic dynamic instability occurs when the maximum Mises stress at corresponding elastic modulus lies above the critical curve. Otherwise, the dynamic response of the blade is stable.

Although some progress has been made in dynamic stability analysis on wind turbine blade in this work, there exist some limitations to the present simulation as follows: 1) aerodynamic loads were calculated by using empirical formula, which cannot represent the real wind field completely; 2) the influences of fluid-structure coupling on the blade vibration were not taken into account.