Analysis of dynamic stability for wind turbine blade under fluid-structure interaction

2.1. Physical model

For the uni-direction fluid-structure interaction (UFSI), the force on the fluid-structure interface, which is applied to the structure to cause its deformation, is got by the calculation on fluid field. Fluid-structure interaction model of three dimensional flow fields is displayed in Fig. 1.

2.2. Control and discrete equations in fluid domain

The basic conservation laws of fluid domain include the mass conservation law, the momentum conservation law and the energy conservation law. In order to analyze these equations easily and solve them by using the same program, the general form of governing equations is established. $\mathbf{\varphi }$ is the general variable, which indicates the direction of velocity component $u^{\mathrm{*}}$, $v^{\mathrm{*}}$, $w^{\mathrm{*}}$, and it can be expressed as [15]:

where, $\rho$, $\mathrm{\Gamma }$, $S$ respectively are the air density, the generalized diffusion coefficient and the generalized source term. For the continuity equation, $\varphi =\mathrm{}$1; for the momentum equation, $\varphi =u_j^{\mathrm{*}}$, where $u_i^{\mathrm{*}}$ and $u_j^{\mathrm{*}}$ are the velocity component in each direction; for the two-equation $k$-$\epsilon$ turbulence model, $\mathbf{\varphi }=\left\{k,\epsilon \right\}$. The indicator values of $i$, $j$ are 1, 2 and 3 respectively, and $\mathrm{\Gamma }$, $S$ are ascertained according to $\mathbf{\varphi }$.

Fig. 1The meshing and interaction diagram of fluid-structure field

The control equation uses finite volume discretization according to the computational grid. Integral operation is conducted on the control volume $P$ and time frame $\mathrm{\Delta }t$, and Guass divergence theorem is introduced. The general form of the discrete equations is as follows [15]:

in which $nb$ is the adjacent node, $nb=\left\{w,e,s,n,b,t\right\}$; $a$ is the coefficient; $b$ comes from source item.

2.3. Wind load calculation

2.3.2. Turbulence determination

First of all, fluctuating wind spectrum needs to be simulated under the certain condition to consider the influence of TE on wind turbine blades. Davenport spectrum is the average value of 90 times vertical turbulent power spectrum based on different locations, different heights around the world, and the expressions are given as [20]:

4

$S_v\left({\omega }_0\right)={{\stackrel{-}{V}}_{10}}^2\frac{4k{{\stackrel{-}{x}}_0}^2}{{\omega }_0{\left(1+{x_0}^2\right)}^{4/3}},$

5

${\stackrel{-}{x}}_0=600\frac{{\omega }_0}{\pi {\stackrel{-}{V}}_{10}},$

where $S_v\left({\omega }_0\right)$ and ${\omega }_0$ are power spectrum and circle frequency of the fluctuating wind, respectively; $k$ is the ground roughness coefficient. According to the theory of random vibration and Chinese load specification, the formula of $k$ can be written as:

6

$k=\frac{1}{6\times 8.8^2}\times 35^{3.6\left(\theta -0.16\right)}.$

The frequency domain was converted into the time domain by the self-programmed program based on the harmonic synthesis method and Eq. (4), and the time domain curve of the fluctuating wind speed was obtained. When ${\stackrel{-}{V}}_{10}=\mathrm{}$20 m/s and ${\omega }_0=\mathrm{}$0.36587 Hz are selected, the variation relationship of fluctuating wind velocity $v\left(t\right)$ with time is simulated and shown in Fig. 2.

Fig. 2Time history of fluctuating wind velocity

In the case of average wind velocity (AWV), TE and WSE, Eq. (2) is solved by iterative method. Velocity field satisfying the continuity equation is obtained, and so wind pressure distribution on the blade surface is calculated.

2.4. Nonlinear dynamics equation and discrete equation in structure domain

The blade meets the geometrically nonlinear equation [21]:

where ${\sigma }_{ij}$ is stress tensor, $U_i$ is the blade displacement vector, $\stackrel{-}{\rho }$ is the blade material density and $f_i$ is the body force vector. The constitutive equation and the relationship between displacement and strain for the geometrically nonlinear elastic body are respectively written as:

in which $G$, $\lambda$ are shear modulus of elasticity and LAME coefficient respectively, and ${\delta }_{ij}$ is the unit tensor.

In general case, the discrete motion differential equation of geometrically nonlinear structure domain is:

where $\left[\mathbf{M}\right]$, $\left[\mathbf{C}\right]$ are mass matrix and damping matrix, respectively; $\left[\mathbf{x}\right]=[u,\mathrm{}\mathrm{}v,\mathrm{}\mathrm{}w{]}^T$; $\left[\mathbf{F}\right(t){]}_{t+\mathrm{\Delta }t}$ is the wind load column vector at $t+\mathrm{\Delta }t$; $[\stackrel{-}{\mathbf{P}}{]}_{t+\mathrm{\Delta }t}$ is the nonlinear term of equation, it can be expressed as:

where $[{\mathbf{K}}_T{]}_t$ is the tangential stiffness matrix calculated by node displacement at $t$.

For the discrete motion differential Eq. (10), the displacement, velocity and acceleration at $t+∆t$ are calculated by Newmark method and Newton-Raphson iterative method [22].

2.5. Modal analysis

When the natural frequency of the blade structure is solved, the classical eigenvalue equation is expressed as [23]:

where $\left\{{\mathbf{\varphi }}_i\right\}$ and ${\omega }_i$ ($i=$ 1, 2, 3,…) are the feature vector and the natural frequency of the $i$th order modal.

3.1. Entity modeling and parameter setting

The correct establishment of entity model directly affects the accuracy and reliability of numerical results in the process of the finite element analysis by using ANSYS software. The design parameters of 1.5 MW wind turbine chosen in this paper are shown in Table 1. The blade entity model set up by the Pro/E software is listed in Fig. 3.

Nowadays, the glass-fiber reinforced plastics (GRPs) are commonly used to manufacture the commercial wind turbine blade, because of their light weight, high strength and good rigidity. Specific material parameters are shown in Table 2 [24].

Fig. 3The entity model of the 1.5 MW wind blade

3.2. Finite element mesh

Model meshing, which directly affects the accuracy of the subsequent numerical results, is the basis for differential equation discretization and numerical solution. The overall grid of fluid-structure field is shown in Fig. 1.

With the increase of the unit grid quantity, it is found that the maximum brandish displacement of the blade is convergent, as shown in Table 3. On the premise of achieving both the calculation precision and the speed, the grid quantities of fluid-structure domain are determined to be 430847 and 9656, respectively.

4.1. Model validation

The numerical program is degenerated into the cantilever blade structure corresponding to the reference [25] in which the maximum Mises stress of blade root portion is 42.70 MPa, while the one calculated in this work is 43.53 MPa as shown in Fig. 4. The relative error is only 1.94 %, which verifies the reliability of the present theoretical method, the entity modeling and finite element method.

4.2. The natural frequencies of wind turbine blade

The actual percentage contribution of each order mode calculated by modal analysis is different from the system brandish vibration, and the brandish vibration mode is mainly concentrated in the first six modes. The first six natural frequencies are listed in Table 4.

Fig. 4Maximum Mises stress contour

4.3. Blade vibration analysis under periodic wind speed

The dynamic stability of the wind turbine blade under the action of fluid structure interaction is emphatically discussed in this work. The stability region is determined by calculating the critical curves of dynamic instability, and the characteristic of instability is the vibration divergence. Therefore, the frequency $\omega$ of the periodic wind speed is set close to the first-order natural frequency ${\omega }_1$ of the blade, and the dynamic stability analysis is carried out with UFSI taken into account. AWV effect ($\omega ={\omega }_1=\mathrm{}$0.36587 Hz), TE (${\omega }_0={\omega }_1=\mathrm{}$0.36587 Hz) and WSE & TE (${\omega }_0={\omega }_1=\mathrm{}$0.36587 Hz) are respectively denoted by Condition I, Condition II and Condition III. The vibration divergence curves of maximum displacement and Mises stress are obtained under above three conditions, as shown in Fig. 5.

Fig. 5 depicts the vibration divergence curves of maximum displacement and Mises stress under three conditions with UFSI considered, and it is clear that:

1) The peak of maximum displacement gradually increases with time and the vibration divergence phenomenon appears under each condition, as shown in Fig. 5(a), (c) and (e).

2) The maximum Mises stress of the blade also increases with time and divergence phenomenon becomes visible (see Fig. 5(b), (d) and (f)).

3) The peak value of the displacement and the Mises stress under Condition III is the largest, which indicates that WSE & TE has the significant impact on the wind turbine blade. Therefore, WSE & TE should not be ignored in the process of numerical simulation, otherwise simulation result will deviate from the real one.

When the vibration divergence of the blade occurs, the relationship among the maximum displacement, the maximum Mises stress and wind speed (${\stackrel{-}{V}}_{}={\stackrel{-}{V}}_{10}$) in the specified computation period under three conditions is shown in Fig. 6 and Fig. 7. It can be seen that:

1) Both the maximum displacement and the Mises stress increase with the increase of wind speed under three conditions and the curve under each condition shows nonlinear upward trend. The value difference between two neighbouring curves also increases, which means that wind speed has significant influences on the blade stability.

2) The displacement and the stress under Condition III show a maximum, and the variation gradient also keep the maximum changing with wind speed, which implies that the TE & WSE has an obvious effect on the stability of the blade.

3) At the same wind speed, the difference in either displacement or stress between Condition II and Condition III represents the contribution from WSE, and the one under Condition I and Condition II from TE. It is obvious that the influence of WSE on displacement and stress is greater than that of TE, and this trend becomes more and more significant with the increase of wind speed.

4.4. Search of blade instability region

In order for the wind turbine to operate safely under actual working condition, it is necessary to ensure that the maximum displacement does not exceed the horizontal distance between the blade and tower. In this study, two horizontal distances of 1.5 m and 1.0 m are chosen as the maximum allowable displacement of the blade to meet the safe design requirements. When the dynamic instability phenomenon takes place, the quantitative relationship between Mises stress and elasticity modulus of blade under three different conditions is determined according to the safe horizontal distance criterion. This relationship provides the safety reference for the materials selection of wind turbine blades.

Fig. 5The divergent response curves of the blade under three different conditions

a) The maximum displacement under Condition I ($\stackrel{-}{V}=\mathrm{}$20m/s, $\omega ={\omega }_1$)

b) The maximum Mises stress under Condition I ($\stackrel{-}{V}=\mathrm{}$20 m/s, $\omega ={\omega }_1$)

c) The maximum displacement under Condition II (${\stackrel{-}{V}}_{10}=\mathrm{}$20 m/s, ${\omega }_0={\omega }_1$)

d) The maximum Mises stress under Condition II (${\stackrel{-}{V}}_{10}=\mathrm{}$20 m/s, ${\omega }_0={\omega }_1$)

e) The maximum displacement under Condition III (${\stackrel{-}{V}}_{10}=\mathrm{}$20 m/s, ${\omega }_0={\omega }_1$)

f) The maximum Mises stress under Condition III (${\stackrel{-}{V}}_{10}=\mathrm{}$20 m/s, ${\omega }_0={\omega }_1$)

Fig. 6The maximum displacement versus wind speed under three conditions

Fig. 7The maximum Mises stress versus wind speed under three conditions

Fig. 8The maximum Mises stress versus elasticity modulus under AWV (V-= 20 m/s, ω=ω1)

Fig. 9The maximum Mises stress versus elasticity modulus under TE (V-10= 20 m/s, ω0=ω1)

Fig. 10The maximum Mises stress versus elasticity modulus under TE & WSE (V-10= 20 m/s, ω0=ω1)

When the vibration diverges, Figs. 8, 9 and 10 represent the critical curves of the maximum Mises stress of the blade versus elastic modulus at two safe horizontal distances 1.5 m and 1.0 m. It can be drawn that:

1) Both two curves under three working conditions are the critical dynamic stability curves which determine whether the wind turbine blade can operate safely. When the maximum Mises stress calculated by material elastic modulus at each condition stays above the critical curve, the elastic dynamic instability happens and the blade cannot operate safely. Otherwise, the dynamic response of the blade is stable.

2) The critical curves under three working conditions have the nonlinear trend with the increase of the elasticity modulus at 1.0 m safe horizontal distance, and the maximum Mises stress is between 3.16 MPa to 13.40 MPa.

3) At 1.5 m safe horizontal distance, the critical curves under three working conditions basically present a linear upward trend with the increase of the elasticity modulus, and the gradient of the critical curve of the maximum Mises stress decreases slightly from Condition I to Condition III.