# Two-Degree Vibration Analysis of a Horizontal Axis Turbine Blade by Finite Differential Methods

Authors: Youfeng Zhu; Hongyu Zang; Xiaofei Kong; Peng Ban
Journal: Shock and Vibration
 (2018)

---

## Abstract

Two-degree vibration partial differential equations of large horizontal axis turbine blades were established by Kallesøe’s model and Greenberg unsteady aerodynamic theory. By means of the finite difference discretization and cantilever beam boundary condition, the equations of blades can be simplified as a general vibration system. Then a linear stationary state space on the system was built. The blade tip vibration in autonomous and nonautonomous system can be simulated by MATLAB vibration toolboxes in time domain. The convergent, flutter, and divergent vibration curves were plotted in the directions of lead-lag and flapping.


---

## Body

## 1. Introduction

Many works [1–3] had been done in dealing with the pretwisted isotropic turbine blades; however, they did not consider the effects of gravity, pitch action, and rotor speed variations. The equations they provide, which are partial differential equations, are convenient to be numerically analyzed. In fact, blades in reality are made of composite materials which are anisotropic and result in the internal elastic coupling between different forms of blade motion, which cannot be solved by the equations discussed above. At present, the large wind turbine blade aeroelastic models contain ONERA [4–6], Beddoes-Leishman [7–10], and Theodorsen [11, 12]. This paper extends Kallesøe’s equations of blade motion to hollow composite materials by equivalent elastic moduli and equivalent densities. According to Greenberg unsteady aerodynamic theory, this paper establishes a new model which describes two-degree aeroelastic equations with asymmetric rigidity and general viscous damping through finite difference discretization.The mode superposition and direct integration are two main research methods [13, 14] of flutter on large horizontal axis wind turbine blades. In the proportional damping or modal damping system, the equations can be decoupled by using a real modal method. Then the system responses are obtained by a modal superposition method. However, there is no orthogonality in the general damping problem; it cannot be solved by a real modal superposition method. It can reduce the state equations into one-order form and solve the general eigenvalue and eigenvector. After the equations are decoupled by complex modal analysis, we can get solutions by complex modal superposition. The state space analysis method does not need coordinate transformation and equation decoupling, making it suitable for general mechanical vibration systems and convenient for the computer analysis in time and frequency domains.


## 2. Establishment of System Equations

### 2.1. Structural Differential Equations

Kallesøe [15] introduced coordinate conversion and differential equation in a model of turbine blades as shown in Figure 1. The blade can be made of solid metal materials or hollow composite materials. Figure 1(a) shows the rotating motion of the blade in the plane of the wind wheel. The coordinate system of (
X
,
Y
,
Z
) is the wind wheel coordinate system of the blade. The x-axis is the horizontal axis, the y-axis coincides with the horizontal axis of the blade of the wind turbine, and the z-axis is vertically upward. Because the tower of the wind turbine and the rotating shaft of the wind turbine are fixed, the coordinate system is a fixed coordinate system. The coordinate system (
x
1
,
y
1
,
z
1
) revolves around the horizontal axis of the blade, the y
1 axis and the y-axis overlap, and the z
1 axis is along the pitch axis of the blade. The angle of two coordinate systems indicates the angle of rotation of the blade.Figure 1
The model of a horizontal axis turbine blade.
(a)
(b)Figure1(b) represents a cross section perpendicular to the z-axis and it is described in the (
x
,
y
,
z
) coordinate system. The coordinate system rotates the β along the Z
1 axis, and the β
=
β
(
s
,
t
) is the pitch angle of the blade. The cross section after the deformation is expressed by (
η
,
ζ
). It rotates θ
+
φ relative to the plane of (
x
,
z
), θ
=
θ
(
s
) is the twist angle of the blade, and φ
=
φ
(
s
,
t
) refers to the twist angle of the blade. The coordinates of the elastic center in the (
x
,
y
,
z
) coordinate system are (
u
+
l
p
i
,
v
,
w
). u
=
u
(
s
,
t
) and v
=
v
(
s
,
t
) represent the deformation in the direction of the x-axis and the y-axis, respectively. l
p
i
=
l
p
i
(
s
) is the distance between the center of the blade and the center of the elasticity.The forces at thex- and y-axes are(1a)
m
u
¨
+
F
u
1
ϕ
˙
,
u
˙
,
v
˙
,
u
′
,
v
,
θ
,
β
+
F
u
2
ϕ
,
u
′
,
v
′
,
θ
,
β
+
F
u
3
u
′
′
,
v
′
′
,
θ
=
f
u
,
(1b)
m
v
¨
+
F
v
1
ϕ
˙
,
u
˙
,
v
˙
,
u
′
,
v
,
θ
,
β
+
F
v
2
ϕ
,
u
′
,
v
′
,
θ
,
β
+
F
v
3
u
′
′
,
v
′
′
,
θ
=
f
v
,where(2a)
F
u
1
=
-
ϕ
˙
2
m
u
c
g
cos
⁡
β
-
ϕ
˙
2
m
l
c
g
w
0
cos
⁡
θ
-
θ
sin
⁡
θ
′
-
cos
⁡
θ
l
c
g
T
2
′
-
2
ϕ
˙
m
l
c
g
cos
⁡
β
u
˙
′
cos
⁡
θ
+
v
˙
′
sin
⁡
θ
-
u
′
+
l
p
i
′
∫
s
R
ϕ
˙
2
m
w
0
+
T
2
d
ρ
′
,
(2b)
F
v
1
=
-
ϕ
˙
2
m
u
c
g
cos
⁡
β
-
ϕ
˙
2
m
l
c
g
w
0
cos
⁡
θ
-
θ
sin
⁡
θ
′
-
cos
⁡
θ
l
c
g
T
2
′
-
2
ϕ
˙
m
l
c
g
cos
⁡
β
u
˙
′
cos
⁡
θ
+
v
˙
′
sin
⁡
θ
-
u
′
+
l
p
i
′
∫
s
R
ϕ
˙
2
m
w
0
+
T
2
d
ρ
′
,where(3)
u
c
g
=
u
+
l
p
i
cos
⁡
β
-
v
sin
⁡
β
+
l
c
g
cos
⁡
θ
+
β
-
l
c
g
φ
sin
⁡
θ
+
β
,
T
2
=
2
m
ϕ
˙
u
˙
cos
⁡
β
-
v
˙
sin
⁡
β
.u
c
g is the x
1 axis component of the center of gravity of the (
x
1
,
y
1
,
z
1
) coordinate system. T
2 is a leaf in the Z-axis direction of the Coriolis force. The first term (2a) refers to the centrifugal force of the rotating blade. The second and third items (2a) are caused by the inconsistency of the center of gravity and the center of the aeroelastic center. Fourth type (2a) is the rotating blades on the Z
1 axis produced by Coriolis force. The last term of formula (2a) is the bending moment of the blade position relative to the end of the blade. The effect of blade gravity is as follows:(4a)
F
μ
2
=
m
g
sin
⁡
ϕ
cos
⁡
β
+
l
c
g
u
′
+
l
p
i
′
′
cos
⁡
θ
cos
⁡
β
+
l
c
g
v
′
′
sin
⁡
θ
cos
⁡
β
m
g
sin
⁡
ϕ
,
(4b)
F
v
2
=
-
m
g
cos
⁡
ϕ
sin
⁡
β
+
l
c
g
u
′
+
l
p
i
′
′
sin
⁡
θ
cos
⁡
β
-
l
c
g
v
′
′
sin
⁡
θ
sin
⁡
β
m
g
sin
⁡
ϕ
.The first term (4a) refers to the weight of gravity in the X-axis. The second item (4a) is the torque produced by gravity. The third item (4a) is the moment caused by the inconsistency of the center of gravity and the center of the aeroelastic. The restoring force produced by the bending stiffness is expressed as follows:(5a)
F
u
3
=
E
ξ
I
ξ
c
o
s
2
θ
+
E
η
I
η
s
i
n
2
θ
u
′
′
′
′
+
E
ξ
I
ξ
-
E
η
I
η
cos
⁡
θ
sin
⁡
θ
v
′
′
′
′
,
(5b)
F
v
3
=
E
ξ
I
ξ
s
i
n
2
θ
+
E
η
I
η
c
o
s
2
θ
u
′
′
′
′
+
E
ξ
I
ξ
-
E
η
I
η
cos
⁡
θ
sin
⁡
θ
v
′
′
′
′
.Among them,(6)
I
η
=
∬
A
ξ
2
d
η
d
ξ
,
I
ξ
=
∬
A
η
2
d
η
d
ξ
.The first item (5a) is the bending stiffness in the direction of the X-axis. The second item (5a) is the coupling bending stiffness in the direction of the Y-axis. The inertial distance can be seen in the upper form.If considering the effect of gravity, according to the first- and second-order Taylor expansions, ignoring higher-order terms and aerodynamic parameters, the motion equations of the blade inx and y directions can be obtained in the following: (7a)
m
u
¨
-
u
′
∫
0
R
Ω
2
m
s
d
ρ
′
+
E
I
ξ
c
o
s
2
θ
+
I
η
s
i
n
2
θ
u
′
′
′
′
+
E
I
ξ
-
I
η
cos
⁡
θ
sin
⁡
θ
v
′
′
′
′
-
Ω
2
m
cos
⁡
β
u
cos
⁡
β
-
v
sin
⁡
β
+
f
u
=
-
m
g
sin
⁡
Ω
t
cos
⁡
β
,
(7b)
m
v
¨
-
v
′
∫
0
R
Ω
2
m
s
d
ρ
′
+
E
I
ξ
s
i
n
2
θ
+
I
η
c
o
s
2
θ
v
′
′
′
′
+
E
I
ξ
-
I
η
cos
⁡
θ
sin
⁡
θ
u
′
′
′
′
+
Ω
2
m
sin
⁡
β
u
cos
⁡
β
-
v
sin
⁡
β
+
f
v
=
m
g
sin
⁡
Ω
t
sin
⁡
β
.


### 2.2. Aerodynamic Equations

Aerodynamic force and torque are caused by the relative motion of the blade and the unsteady aerodynamic formula is given by Greenberg. The relative velocity of the blade to the air is deduced below and put it into the aerodynamic formula. For convenience, the fourth coordinate systems are introduced(
x
2
,
y
2
,
z
2
) as shown in Figure 2.(8)
i
2
j
2
k
2
=
T
a
i
j
k
=
1
0
-
u
′
0
1
-
v
′
u
′
v
′
1
i
j
k
.Figure 2
Coordinates of the blade.The relative velocity of the blade to the air should be the vector sum of the absolute velocity and the velocity of the air flow. That is,(9)
V
=
ω
×
r
+
r
˙
-
v
i
n
.The velocity of air is(10)
v
i
n
=
λ
Ω
r
K
=
0
λ
Ω
r
0
i
j
k
T
.Andλ
=
v
i
n
/
Ω
r.After vector calculation and coordinate transformation, the following can be obtained:(11)
V
=
v
x
2
i
2
+
v
y
2
j
2
,where(12)
v
x
2
=
u
˙
+
Ω
s
,
v
y
2
=
v
˙
+
λ
Ω
r
+
Ω
u
v
′
.The aerodynamic force of theZ direction is ignored because it has little effect on the analysis of the vibration of the aeroelasticity. Next, the mod of V is(13)
V
=
V
=
Ω
s
+
λ
2
Ω
2
s
+
u
˙
+
λ
v
˙
.In order to apply these speeds to the need of aerodynamic formulas, the relative velocity is divided into static and dynamic parts first:(14)
V
0
=
Ω
s
i
2
+
λ
Ω
s
+
Ω
v
′
u
j
2
,
V
d
≈
u
˙
i
2
+
v
˙
j
2
.Next, with the static speedV
0, normal n, and tangent direction t being coordinates and n and t being unit vectors, the dynamic velocity is decomposed as shown in Figure 3.(15)
1
V
0
≈
1
Ω
s
1
-
0.5
λ
2
,
t
=
V
V
0
=
1
-
0.5
λ
2
i
2
+
λ
j
2
.Figure 3
Velocity vectors of the blade.The unit vectorn of the normal phase is perpendicular to the tangent unit vector t:(16)
n
=
-
λ
i
2
+
1
-
0.5
λ
2
j
2
.With the tangent and normal unit vector being available, the following can be obtained:(17)
h
˙
=
λ
u
˙
-
1
-
0.5
λ
2
v
˙
,
h
¨
=
λ
u
¨
-
1
-
0.5
λ
2
v
¨
,
p
˙
=
1
-
0.5
λ
2
u
˙
+
λ
v
˙
,
p
¨
=
1
-
0.5
λ
2
u
¨
+
λ
v
¨
.The total tangential velocity can be obtained by adding the static velocityv
0 to the dynamic velocity:(18)
v
t
=
v
t
x
2
i
2
+
v
t
y
2
j
2
=
Ω
s
+
u
˙
+
λ
v
˙
i
2
+
λ
Ω
s
+
Ω
v
′
u
+
λ
u
˙
j
2
.In addition to the velocity above, the angle of attack is also given before the aerodynamic load is calculated and we also need to find the mod and the reciprocal ofv
t:(19)
v
t
=
v
t
≈
Ω
s
+
0.5
λ
2
Ω
s
+
u
˙
,
1
v
t
≈
1
Ω
s
1
-
0.5
λ
2
-
u
˙
/
Ω
s
.As shown in Figure4,(20)
v
t
x
2
v
t
=
cos
⁡
θ
i
-
α
=
cos
⁡
θ
i
cos
⁡
α
+
sin
⁡
θ
i
sin
⁡
α
,
v
t
y
2
v
t
=
sin
⁡
θ
i
-
α
=
sin
⁡
θ
i
cos
⁡
α
-
cos
⁡
θ
i
sin
⁡
α
.Figure 4
Pitch angel of the blade.The aerodynamic expression given by the Greenberg aerodynamic force theory is(21)
F
L
n
=
ρ
b
c
2
8
h
¨
+
v
t
α
˙
+
v
˙
t
α
+
c
4
α
¨
,
F
L
c
=
ρ
b
c
v
t
2
v
t
α
+
c
2
α
˙
+
h
˙
,
F
D
=
ρ
c
c
d
u
2
2
.Acyclic lift of aerodynamic forceF
L
n and circulation lift F
L
c are perpendicular to the blade string and the velocity vector u. The resistance F
D is parallel to the velocity vector u, but the direction is opposite. The torque is acting on the z
2 axis.(22)
F
x
2
=
-
F
D
v
x
2
u
-
F
L
c
v
y
2
u
-
F
L
n
sin
⁡
θ
i
,
F
y
2
=
-
F
D
v
y
2
u
+
F
L
c
v
x
2
u
+
F
L
n
cos
⁡
θ
i
.The aerodynamic lift with circulation, acyclic lift, and drag can be gotten by using Greenberg aerodynamic theory. The changes of inertia and stiffness as a result of aerodynamic forces have little effect on the stability of turbine blades, so they can be ignored. According to the coordinate conversion and high order simplification, the aerodynamic forces can be derived on small angle of attack wind turbine blade inx and y directions:(23a)
f
u
=
0.5
ρ
a
c
sin
⁡
θ
v
i
n
cos
⁡
β
u
˙
+
0.5
ρ
a
c
v
i
n
cos
⁡
β
v
i
n
Ω
s
1
-
c
o
s
θ
u
˙
+
0.125
ρ
a
c
2
sin
⁡
θ
Ω
β
u
˙
+
0.5
ρ
a
c
sin
⁡
θ
Ω
s
cos
⁡
β
v
˙
+
0.5
ρ
a
c
sin
⁡
θ
v
i
n
cos
⁡
β
v
i
n
Ω
s
v
˙
-
0.5
ρ
a
c
v
i
n
cos
⁡
β
v
˙
+
0.25
ρ
a
c
2
Ω
β
v
˙
,
(23b)
f
v
=
-
ρ
a
c
sin
⁡
θ
Ω
s
cos
⁡
β
u
˙
+
0.5
ρ
a
c
v
i
n
cos
⁡
β
u
˙
-
0.375
ρ
a
c
2
Ω
β
u
˙
+
0.5
ρ
a
c
Ω
s
cos
⁡
β
1
-
0.5
v
i
n
Ω
s
2
v
˙
-
0.5
ρ
a
c
sin
⁡
θ
Ω
s
cos
⁡
β
v
i
n
Ω
s
v
˙
,where m
=
m
(
s
), θ
=
θ
(
s
), I
ζ
=
I
ζ
(
s
), I
η
=
I
η
(
s
), C
=
C
(
s
).


### 2.3. Finite Difference Method

This paper uses the fourth-order derivative inx-direction to show the process of finite different method. The equations in y-direction are the same as those in x-direction.Before the blade model by using finite difference method being simplified, the technical parameters of wind turbine blades must be solved through differential solution or curve fitting to get values of different nodes. The boundary conditions are expressed as follows:(24)
u
0
,
t
=
u
′
0
,
t
=
v
0
,
t
=
v
′
0
,
t
=
0
,
u
′
′
R
,
t
=
v
′
′
R
,
t
=
u
′
′
′
R
,
t
=
v
′
′
′
R
,
t
=
0
.Based on Taylor expansions, the general difference form of each derivative can be obtained through two-order central difference. According to the first-order central difference formula and (24), (25) can be obtained:(25)
u
0
′
=
1
2
h
-
u
-
1
+
u
1
=
0
,
u
-
1
=
u
1
.According to the second- and third-order difference formulas and (24), (26) can be obtained:(26)
u
n
′
′
=
1
12
h
2
-
u
n
-
2
+
16
u
n
-
1
-
30
u
n
+
16
u
n
+
1
-
u
n
+
2
=
0
,
u
n
′
′
′
=
1
2
h
3
-
u
n
-
2
+
2
u
n
-
1
-
2
u
n
+
1
+
u
n
+
2
=
0
.Equation (27) can be obtained by substituting the boundary conditions into (26):(27)
u
n
+
1
=
1
7
u
n
-
2
-
9
u
n
-
1
+
15
u
n
,
u
n
+
2
=
1
7
9
u
n
-
2
-
32
u
n
-
1
+
30
u
n
.Finally, the fourth-order derivative difference equations are as follows:(28)
u
2
′
′
′
′
=
1
h
4
-
4
u
1
+
6
u
2
-
4
u
3
+
u
4
,
u
n
′
′
′
′
=
1
h
4
u
n
-
2
-
4
u
n
-
1
+
6
u
n
-
4
u
n
+
1
+
u
n
+
2
,
u
N
-
1
′
′
′
′
=
1
7
h
4
7
u
N
-
3
-
27
u
N
-
2
+
33
u
N
-
1
-
13
u
N
,
u
N
′
′
′
′
=
12
7
h
4
u
N
-
2
-
2
u
N
-
1
+
u
N
,
3
≤
n
≤
N
-
2
.Similarly, the first-order derivative difference equations are as follows:(29)
u
1
′
=
1
12
h
u
1
+
8
u
2
-
u
3
,
u
2
′
=
1
12
h
-
8
u
1
+
8
u
3
-
u
4
,
u
n
′
=
1
12
h
u
n
-
2
-
8
u
n
-
1
+
8
u
n
+
1
-
u
n
+
2
,
u
N
-
1
′
=
1
12
h
u
N
-
3
-
57
7
u
N
-
2
+
9
7
u
N
-
1
+
41
7
u
N
,
u
N
′
=
1
12
h
6
7
u
N
-
2
-
96
7
u
N
-
1
+
90
7
u
N
,
3
≤
n
≤
N
-
2
.The second-order derivative difference equations are as follows:(30)
u
1
′
′
=
1
12
h
2
-
31
u
1
+
16
u
2
-
u
3
,
u
2
′
′
=
1
12
h
2
16
u
1
-
30
u
2
+
16
u
3
-
u
4
,
u
n
′
′
=
1
12
h
2
-
u
n
-
2
+
16
u
n
-
1
-
30
u
n
+
16
u
n
+
1
-
u
n
+
2
,
u
N
-
1
′
′
=
1
84
h
2
-
7
u
N
-
3
+
111
u
N
-
2
-
201
u
N
-
1
+
97
u
N
,
u
N
′
′
=
0
,
3
≤
n
≤
N
-
2
.The third-order derivative difference equations are as follows:(31)
u
1
′
′
′
=
1
2
h
3
-
u
1
-
2
u
2
+
u
3
,
u
2
′
′
′
=
1
2
h
3
2
u
1
-
2
u
3
+
u
4
,
u
n
′
′
′
=
1
2
h
3
-
u
n
-
2
+
2
u
n
-
1
-
2
u
n
+
1
+
u
n
+
2
,
u
N
-
1
′
′
′
=
1
14
h
3
-
7
u
N
-
3
+
15
u
N
-
2
-
9
u
N
-
1
+
u
N
,
u
N
′
′
′
=
0
,
3
≤
n
≤
N
-
2
.


## 2.1. Structural Differential Equations

Kallesøe [15] introduced coordinate conversion and differential equation in a model of turbine blades as shown in Figure 1. The blade can be made of solid metal materials or hollow composite materials. Figure 1(a) shows the rotating motion of the blade in the plane of the wind wheel. The coordinate system of (
X
,
Y
,
Z
) is the wind wheel coordinate system of the blade. The x-axis is the horizontal axis, the y-axis coincides with the horizontal axis of the blade of the wind turbine, and the z-axis is vertically upward. Because the tower of the wind turbine and the rotating shaft of the wind turbine are fixed, the coordinate system is a fixed coordinate system. The coordinate system (
x
1
,
y
1
,
z
1
) revolves around the horizontal axis of the blade, the y
1 axis and the y-axis overlap, and the z
1 axis is along the pitch axis of the blade. The angle of two coordinate systems indicates the angle of rotation of the blade.Figure 1
The model of a horizontal axis turbine blade.
(a)
(b)Figure1(b) represents a cross section perpendicular to the z-axis and it is described in the (
x
,
y
,
z
) coordinate system. The coordinate system rotates the β along the Z
1 axis, and the β
=
β
(
s
,
t
) is the pitch angle of the blade. The cross section after the deformation is expressed by (
η
,
ζ
). It rotates θ
+
φ relative to the plane of (
x
,
z
), θ
=
θ
(
s
) is the twist angle of the blade, and φ
=
φ
(
s
,
t
) refers to the twist angle of the blade. The coordinates of the elastic center in the (
x
,
y
,
z
) coordinate system are (
u
+
l
p
i
,
v
,
w
). u
=
u
(
s
,
t
) and v
=
v
(
s
,
t
) represent the deformation in the direction of the x-axis and the y-axis, respectively. l
p
i
=
l
p
i
(
s
) is the distance between the center of the blade and the center of the elasticity.The forces at thex- and y-axes are(1a)
m
u
¨
+
F
u
1
ϕ
˙
,
u
˙
,
v
˙
,
u
′
,
v
,
θ
,
β
+
F
u
2
ϕ
,
u
′
,
v
′
,
θ
,
β
+
F
u
3
u
′
′
,
v
′
′
,
θ
=
f
u
,
(1b)
m
v
¨
+
F
v
1
ϕ
˙
,
u
˙
,
v
˙
,
u
′
,
v
,
θ
,
β
+
F
v
2
ϕ
,
u
′
,
v
′
,
θ
,
β
+
F
v
3
u
′
′
,
v
′
′
,
θ
=
f
v
,where(2a)
F
u
1
=
-
ϕ
˙
2
m
u
c
g
cos
⁡
β
-
ϕ
˙
2
m
l
c
g
w
0
cos
⁡
θ
-
θ
sin
⁡
θ
′
-
cos
⁡
θ
l
c
g
T
2
′
-
2
ϕ
˙
m
l
c
g
cos
⁡
β
u
˙
′
cos
⁡
θ
+
v
˙
′
sin
⁡
θ
-
u
′
+
l
p
i
′
∫
s
R
ϕ
˙
2
m
w
0
+
T
2
d
ρ
′
,
(2b)
F
v
1
=
-
ϕ
˙
2
m
u
c
g
cos
⁡
β
-
ϕ
˙
2
m
l
c
g
w
0
cos
⁡
θ
-
θ
sin
⁡
θ
′
-
cos
⁡
θ
l
c
g
T
2
′
-
2
ϕ
˙
m
l
c
g
cos
⁡
β
u
˙
′
cos
⁡
θ
+
v
˙
′
sin
⁡
θ
-
u
′
+
l
p
i
′
∫
s
R
ϕ
˙
2
m
w
0
+
T
2
d
ρ
′
,where(3)
u
c
g
=
u
+
l
p
i
cos
⁡
β
-
v
sin
⁡
β
+
l
c
g
cos
⁡
θ
+
β
-
l
c
g
φ
sin
⁡
θ
+
β
,
T
2
=
2
m
ϕ
˙
u
˙
cos
⁡
β
-
v
˙
sin
⁡
β
.u
c
g is the x
1 axis component of the center of gravity of the (
x
1
,
y
1
,
z
1
) coordinate system. T
2 is a leaf in the Z-axis direction of the Coriolis force. The first term (2a) refers to the centrifugal force of the rotating blade. The second and third items (2a) are caused by the inconsistency of the center of gravity and the center of the aeroelastic center. Fourth type (2a) is the rotating blades on the Z
1 axis produced by Coriolis force. The last term of formula (2a) is the bending moment of the blade position relative to the end of the blade. The effect of blade gravity is as follows:(4a)
F
μ
2
=
m
g
sin
⁡
ϕ
cos
⁡
β
+
l
c
g
u
′
+
l
p
i
′
′
cos
⁡
θ
cos
⁡
β
+
l
c
g
v
′
′
sin
⁡
θ
cos
⁡
β
m
g
sin
⁡
ϕ
,
(4b)
F
v
2
=
-
m
g
cos
⁡
ϕ
sin
⁡
β
+
l
c
g
u
′
+
l
p
i
′
′
sin
⁡
θ
cos
⁡
β
-
l
c
g
v
′
′
sin
⁡
θ
sin
⁡
β
m
g
sin
⁡
ϕ
.The first term (4a) refers to the weight of gravity in the X-axis. The second item (4a) is the torque produced by gravity. The third item (4a) is the moment caused by the inconsistency of the center of gravity and the center of the aeroelastic. The restoring force produced by the bending stiffness is expressed as follows:(5a)
F
u
3
=
E
ξ
I
ξ
c
o
s
2
θ
+
E
η
I
η
s
i
n
2
θ
u
′
′
′
′
+
E
ξ
I
ξ
-
E
η
I
η
cos
⁡
θ
sin
⁡
θ
v
′
′
′
′
,
(5b)
F
v
3
=
E
ξ
I
ξ
s
i
n
2
θ
+
E
η
I
η
c
o
s
2
θ
u
′
′
′
′
+
E
ξ
I
ξ
-
E
η
I
η
cos
⁡
θ
sin
⁡
θ
v
′
′
′
′
.Among them,(6)
I
η
=
∬
A
ξ
2
d
η
d
ξ
,
I
ξ
=
∬
A
η
2
d
η
d
ξ
.The first item (5a) is the bending stiffness in the direction of the X-axis. The second item (5a) is the coupling bending stiffness in the direction of the Y-axis. The inertial distance can be seen in the upper form.If considering the effect of gravity, according to the first- and second-order Taylor expansions, ignoring higher-order terms and aerodynamic parameters, the motion equations of the blade inx and y directions can be obtained in the following: (7a)
m
u
¨
-
u
′
∫
0
R
Ω
2
m
s
d
ρ
′
+
E
I
ξ
c
o
s
2
θ
+
I
η
s
i
n
2
θ
u
′
′
′
′
+
E
I
ξ
-
I
η
cos
⁡
θ
sin
⁡
θ
v
′
′
′
′
-
Ω
2
m
cos
⁡
β
u
cos
⁡
β
-
v
sin
⁡
β
+
f
u
=
-
m
g
sin
⁡
Ω
t
cos
⁡
β
,
(7b)
m
v
¨
-
v
′
∫
0
R
Ω
2
m
s
d
ρ
′
+
E
I
ξ
s
i
n
2
θ
+
I
η
c
o
s
2
θ
v
′
′
′
′
+
E
I
ξ
-
I
η
cos
⁡
θ
sin
⁡
θ
u
′
′
′
′
+
Ω
2
m
sin
⁡
β
u
cos
⁡
β
-
v
sin
⁡
β
+
f
v
=
m
g
sin
⁡
Ω
t
sin
⁡
β
.


## 2.2. Aerodynamic Equations

Aerodynamic force and torque are caused by the relative motion of the blade and the unsteady aerodynamic formula is given by Greenberg. The relative velocity of the blade to the air is deduced below and put it into the aerodynamic formula. For convenience, the fourth coordinate systems are introduced(
x
2
,
y
2
,
z
2
) as shown in Figure 2.(8)
i
2
j
2
k
2
=
T
a
i
j
k
=
1
0
-
u
′
0
1
-
v
′
u
′
v
′
1
i
j
k
.Figure 2
Coordinates of the blade.The relative velocity of the blade to the air should be the vector sum of the absolute velocity and the velocity of the air flow. That is,(9)
V
=
ω
×
r
+
r
˙
-
v
i
n
.The velocity of air is(10)
v
i
n
=
λ
Ω
r
K
=
0
λ
Ω
r
0
i
j
k
T
.Andλ
=
v
i
n
/
Ω
r.After vector calculation and coordinate transformation, the following can be obtained:(11)
V
=
v
x
2
i
2
+
v
y
2
j
2
,where(12)
v
x
2
=
u
˙
+
Ω
s
,
v
y
2
=
v
˙
+
λ
Ω
r
+
Ω
u
v
′
.The aerodynamic force of theZ direction is ignored because it has little effect on the analysis of the vibration of the aeroelasticity. Next, the mod of V is(13)
V
=
V
=
Ω
s
+
λ
2
Ω
2
s
+
u
˙
+
λ
v
˙
.In order to apply these speeds to the need of aerodynamic formulas, the relative velocity is divided into static and dynamic parts first:(14)
V
0
=
Ω
s
i
2
+
λ
Ω
s
+
Ω
v
′
u
j
2
,
V
d
≈
u
˙
i
2
+
v
˙
j
2
.Next, with the static speedV
0, normal n, and tangent direction t being coordinates and n and t being unit vectors, the dynamic velocity is decomposed as shown in Figure 3.(15)
1
V
0
≈
1
Ω
s
1
-
0.5
λ
2
,
t
=
V
V
0
=
1
-
0.5
λ
2
i
2
+
λ
j
2
.Figure 3
Velocity vectors of the blade.The unit vectorn of the normal phase is perpendicular to the tangent unit vector t:(16)
n
=
-
λ
i
2
+
1
-
0.5
λ
2
j
2
.With the tangent and normal unit vector being available, the following can be obtained:(17)
h
˙
=
λ
u
˙
-
1
-
0.5
λ
2
v
˙
,
h
¨
=
λ
u
¨
-
1
-
0.5
λ
2
v
¨
,
p
˙
=
1
-
0.5
λ
2
u
˙
+
λ
v
˙
,
p
¨
=
1
-
0.5
λ
2
u
¨
+
λ
v
¨
.The total tangential velocity can be obtained by adding the static velocityv
0 to the dynamic velocity:(18)
v
t
=
v
t
x
2
i
2
+
v
t
y
2
j
2
=
Ω
s
+
u
˙
+
λ
v
˙
i
2
+
λ
Ω
s
+
Ω
v
′
u
+
λ
u
˙
j
2
.In addition to the velocity above, the angle of attack is also given before the aerodynamic load is calculated and we also need to find the mod and the reciprocal ofv
t:(19)
v
t
=
v
t
≈
Ω
s
+
0.5
λ
2
Ω
s
+
u
˙
,
1
v
t
≈
1
Ω
s
1
-
0.5
λ
2
-
u
˙
/
Ω
s
.As shown in Figure4,(20)
v
t
x
2
v
t
=
cos
⁡
θ
i
-
α
=
cos
⁡
θ
i
cos
⁡
α
+
sin
⁡
θ
i
sin
⁡
α
,
v
t
y
2
v
t
=
sin
⁡
θ
i
-
α
=
sin
⁡
θ
i
cos
⁡
α
-
cos
⁡
θ
i
sin
⁡
α
.Figure 4
Pitch angel of the blade.The aerodynamic expression given by the Greenberg aerodynamic force theory is(21)
F
L
n
=
ρ
b
c
2
8
h
¨
+
v
t
α
˙
+
v
˙
t
α
+
c
4
α
¨
,
F
L
c
=
ρ
b
c
v
t
2
v
t
α
+
c
2
α
˙
+
h
˙
,
F
D
=
ρ
c
c
d
u
2
2
.Acyclic lift of aerodynamic forceF
L
n and circulation lift F
L
c are perpendicular to the blade string and the velocity vector u. The resistance F
D is parallel to the velocity vector u, but the direction is opposite. The torque is acting on the z
2 axis.(22)
F
x
2
=
-
F
D
v
x
2
u
-
F
L
c
v
y
2
u
-
F
L
n
sin
⁡
θ
i
,
F
y
2
=
-
F
D
v
y
2
u
+
F
L
c
v
x
2
u
+
F
L
n
cos
⁡
θ
i
.The aerodynamic lift with circulation, acyclic lift, and drag can be gotten by using Greenberg aerodynamic theory. The changes of inertia and stiffness as a result of aerodynamic forces have little effect on the stability of turbine blades, so they can be ignored. According to the coordinate conversion and high order simplification, the aerodynamic forces can be derived on small angle of attack wind turbine blade inx and y directions:(23a)
f
u
=
0.5
ρ
a
c
sin
⁡
θ
v
i
n
cos
⁡
β
u
˙
+
0.5
ρ
a
c
v
i
n
cos
⁡
β
v
i
n
Ω
s
1
-
c
o
s
θ
u
˙
+
0.125
ρ
a
c
2
sin
⁡
θ
Ω
β
u
˙
+
0.5
ρ
a
c
sin
⁡
θ
Ω
s
cos
⁡
β
v
˙
+
0.5
ρ
a
c
sin
⁡
θ
v
i
n
cos
⁡
β
v
i
n
Ω
s
v
˙
-
0.5
ρ
a
c
v
i
n
cos
⁡
β
v
˙
+
0.25
ρ
a
c
2
Ω
β
v
˙
,
(23b)
f
v
=
-
ρ
a
c
sin
⁡
θ
Ω
s
cos
⁡
β
u
˙
+
0.5
ρ
a
c
v
i
n
cos
⁡
β
u
˙
-
0.375
ρ
a
c
2
Ω
β
u
˙
+
0.5
ρ
a
c
Ω
s
cos
⁡
β
1
-
0.5
v
i
n
Ω
s
2
v
˙
-
0.5
ρ
a
c
sin
⁡
θ
Ω
s
cos
⁡
β
v
i
n
Ω
s
v
˙
,where m
=
m
(
s
), θ
=
θ
(
s
), I
ζ
=
I
ζ
(
s
), I
η
=
I
η
(
s
), C
=
C
(
s
).


## 2.3. Finite Difference Method

This paper uses the fourth-order derivative inx-direction to show the process of finite different method. The equations in y-direction are the same as those in x-direction.Before the blade model by using finite difference method being simplified, the technical parameters of wind turbine blades must be solved through differential solution or curve fitting to get values of different nodes. The boundary conditions are expressed as follows:(24)
u
0
,
t
=
u
′
0
,
t
=
v
0
,
t
=
v
′
0
,
t
=
0
,
u
′
′
R
,
t
=
v
′
′
R
,
t
=
u
′
′
′
R
,
t
=
v
′
′
′
R
,
t
=
0
.Based on Taylor expansions, the general difference form of each derivative can be obtained through two-order central difference. According to the first-order central difference formula and (24), (25) can be obtained:(25)
u
0
′
=
1
2
h
-
u
-
1
+
u
1
=
0
,
u
-
1
=
u
1
.According to the second- and third-order difference formulas and (24), (26) can be obtained:(26)
u
n
′
′
=
1
12
h
2
-
u
n
-
2
+
16
u
n
-
1
-
30
u
n
+
16
u
n
+
1
-
u
n
+
2
=
0
,
u
n
′
′
′
=
1
2
h
3
-
u
n
-
2
+
2
u
n
-
1
-
2
u
n
+
1
+
u
n
+
2
=
0
.Equation (27) can be obtained by substituting the boundary conditions into (26):(27)
u
n
+
1
=
1
7
u
n
-
2
-
9
u
n
-
1
+
15
u
n
,
u
n
+
2
=
1
7
9
u
n
-
2
-
32
u
n
-
1
+
30
u
n
.Finally, the fourth-order derivative difference equations are as follows:(28)
u
2
′
′
′
′
=
1
h
4
-
4
u
1
+
6
u
2
-
4
u
3
+
u
4
,
u
n
′
′
′
′
=
1
h
4
u
n
-
2
-
4
u
n
-
1
+
6
u
n
-
4
u
n
+
1
+
u
n
+
2
,
u
N
-
1
′
′
′
′
=
1
7
h
4
7
u
N
-
3
-
27
u
N
-
2
+
33
u
N
-
1
-
13
u
N
,
u
N
′
′
′
′
=
12
7
h
4
u
N
-
2
-
2
u
N
-
1
+
u
N
,
3
≤
n
≤
N
-
2
.Similarly, the first-order derivative difference equations are as follows:(29)
u
1
′
=
1
12
h
u
1
+
8
u
2
-
u
3
,
u
2
′
=
1
12
h
-
8
u
1
+
8
u
3
-
u
4
,
u
n
′
=
1
12
h
u
n
-
2
-
8
u
n
-
1
+
8
u
n
+
1
-
u
n
+
2
,
u
N
-
1
′
=
1
12
h
u
N
-
3
-
57
7
u
N
-
2
+
9
7
u
N
-
1
+
41
7
u
N
,
u
N
′
=
1
12
h
6
7
u
N
-
2
-
96
7
u
N
-
1
+
90
7
u
N
,
3
≤
n
≤
N
-
2
.The second-order derivative difference equations are as follows:(30)
u
1
′
′
=
1
12
h
2
-
31
u
1
+
16
u
2
-
u
3
,
u
2
′
′
=
1
12
h
2
16
u
1
-
30
u
2
+
16
u
3
-
u
4
,
u
n
′
′
=
1
12
h
2
-
u
n
-
2
+
16
u
n
-
1
-
30
u
n
+
16
u
n
+
1
-
u
n
+
2
,
u
N
-
1
′
′
=
1
84
h
2
-
7
u
N
-
3
+
111
u
N
-
2
-
201
u
N
-
1
+
97
u
N
,
u
N
′
′
=
0
,
3
≤
n
≤
N
-
2
.The third-order derivative difference equations are as follows:(31)
u
1
′
′
′
=
1
2
h
3
-
u
1
-
2
u
2
+
u
3
,
u
2
′
′
′
=
1
2
h
3
2
u
1
-
2
u
3
+
u
4
,
u
n
′
′
′
=
1
2
h
3
-
u
n
-
2
+
2
u
n
-
1
-
2
u
n
+
1
+
u
n
+
2
,
u
N
-
1
′
′
′
=
1
14
h
3
-
7
u
N
-
3
+
15
u
N
-
2
-
9
u
N
-
1
+
u
N
,
u
N
′
′
′
=
0
,
3
≤
n
≤
N
-
2
.


## 3. Modal Analysis

In order to easily calculate and analyze the blade model, this paper chose the parameters provided in [16] for the following calculations. As the taper cross-sectional blades are general, this paper chose the blade which is made of isotropic materials for the following simulation. The blade parameters are as follows:(32)
E
ζ
=
E
η
=
10
×
10
9
P
a
,
G
=
10
×
10
7
P
a
,
β
=
0.035
r
a
d
,
ρ
=
1.237
K
g
/
m
3
,
R
=
60
m
,
α
=
0.3
,
θ
=
2
π
s
-
60
2
120
360
,
c
=
4
-
3
s
×
2.8
180
,
b
=
0.25
c
,
m
=
ρ
0
b
c
.The differential equations are substituted in the differential items of (7a) and (7b). According to the chosen node number, the relevant parameters are synthesized. The expressions can be obtained finally in (33). M, K are the mass and stiffness matrices. And q is [
μ
1
,
μ
2
,
μ
3
,
…
,
μ
N
,
ν
1
,
ν
2
,
ν
3
,
…
,
ν
N
].Figure5 shows the first-order modes of vibration in the directions of lead-lad and flapping. They can be solved by vtb4-3 in the MATLAB vibration toolboxes. The first-, third-, and fifth-order modal analyses are made by the finite element method in Figure 6. The first six order natural frequencies with variable finite difference nodes are obtained by FDM in Table 1. Table 1 also shows that the natural frequencies become more accurate with increasing finite difference nodes. In order to reduce the computing time, the number of nodes N
=
16 is chosen in the following calculations. The first six order natural frequencies of blades are derived by FEM and FDM as shown in Table 2. A comparison of two results shows the validity of the model and algorithm.(33)
M
q
¨
+
K
q
=
0
.Table 1
Natural frequencies with variable finite difference nodes.
Number
6
10
16
30
50
ω
1
0.358907
0.378943
0.385872
0.385882
0.385887
ω
2
1.423122
1.234436
1.229124
1.229131
1.229135
ω
3
2.564390
2.379256
2.353636
2.353641
2.353639
ω
4
4.563212
4.136932
4.236893
4.236887
4.236885
ω
5
8.536790
7.315466
7.414138
7.414123
7.414121
ω
6
12.34890
10.06764
10.98524
10.98518
10.98519Table 2
Natural frequencies of blades.
Rotating speed (rad/s)
0
0.8
10
30
ω
1 (FEM)
0.3875
0.4256
1.5765
3.3946
ω
1 (FDM)
0.3859
0.4215
1.5552
3.2976
ω
2 (FEM)
1.2306
1.2393
2.2422
3.4022
ω
2 (FDM)
1.2291
1.2356
2.1287
3.3486
ω
3 (FEM)
2.4218
2.4488
4.7605
9.0082
ω
3 (FDM)
2.3536
2.3214
4.4561
8.6540
ω
4 (FEM)
4.3346
4.3503
5.6395
12.039
ω
4 (FDM)
4.2369
4.5123
6.1420
11.204
ω
5 (FEM)
7.3653
7.3922
10.498
18.218
ω
5 (FDM)
7.4141
7.4345
10.785
17.214
ω
6 (FEM)
10.549
10.555
11.365
20.293
ω
6 (FDM)
10.985
11.241
11.754
22.314Figure 5
The mode of vibration.
(a)
(b)Figure 6
Modal analysis by FEM.
(a)
(b)
(c)


## 4. The Autonomous System

### 4.1. Stability of the Blade

After the simplification by finite difference method, if ignoring the effect of gravity, the turbine blades discussed above can be regarded as an autonomous system. The differential equations are substituted in the differential items of (1a), (1b), (2a), (2b), (4a), (4b), (5a), and (5b). According to the chosen node number, the relevant parameters are synthesized. The expressions in the autonomous system can be obtained finally in the following:(34)
x
¨
j
=
X
j
x
1
,
x
2
,
…
,
x
2
N
,
x
˙
1
,
x
˙
2
,
…
,
x
˙
2
N
,
j
=
1,2
,
…
,
2
N
,
x
=
u
1
,
…
,
u
N
;
v
1
,
…
,
v
N
.Supposingy
j
=
x
˙
j, (34) is transformed into the following:(35)
y
˙
j
=
Y
j
y
1
,
y
2
,
…
,
y
4
N
,
j
=
1,2
,
…
,
4
N
,
y
=
u
1
,
…
,
u
N
;
v
1
,
…
,
v
N
;
u
˙
1
,
…
,
u
˙
N
;
v
˙
1
,
…
,
v
˙
N
.After wind turbine blades are simplified, they become a linear time-invariant system. So the equations can easily be solved by state space expressions in (36). The blade parameters can be optimized by the linear state space. But any real system is always working under all kinds of occasional and continuous interference. After the system is subjected to interference, it can return to the stable work. This is a nonlinear stability problem and is also the focus of future research.(36)
x
˙
=
A
x
+
B
u
,
y
=
C
x
+
D
u
.If all poles in the state space fall on the left plane, the system is stable or unstable. In order to study the relationship between the variable parameters and stability, the maximal real part of pole is introduced. If the maximal real part of the pole is less than zero, the system is stable or unstable. The stable regions of the variable parameters such as the wind velocity, rotational speed, elastic modulus in lead-lag, elastic modulus in flapping, density of the blade, and the pitch angle are shown in Figure7. Taking Figure 7(a) as an example, we select the design and operating parameters of the blade according to experience and keep wind speed as the only variable and then choose different wind speeds to calculate the eigenvalues of the system. According to the criterion theorem of eigenvalues, if all the real part of eigenvalue is less than 0, the system is stable, or the system is unstable. In order to make it convenient for programming, we estimate the stability of system by using the real maximum value of the eigenvalues. The stability area of wind speed of the wind turbine can be determined in this way to provide basis for the design and usage of the wind turbine. The principle of (b) to (f) in Figure 7 is the same as that in Figure 7(a).Figure 7
Stability of variable turbine blade parameters.
(a)
(b)
(c)
(d)
(e)
(f)The flutter speeds in different rotational speeds are also derived in Figure8. Flutter is a critical state of convergence and divergence, that is, a demarcation point of stability or instability. The relationship between wind speed and rotation speed can be obtained in Figure 8. For example, when the wind speed is 20 m/s, the rotation speed of wind turbine should be less than 0.4 rad/s. Once it is larger than this value, the wind turbine will be unstable. This provides a basis for the operation of wind turbine.Figure 8
The flutter speeds in different rotational speeds.


### 4.2. Numerical Simulation

Assisted by MATLAB vibration toolboxes, the vibration on the turbine blade tip can be simulated in the nonautonomous system through time disperse. Select the number of nodesN
=
16, spin speed Ω
=
0.8 rad/s, and the initial displacements of blade tip to be μ
1
6
=
-
0.1885 and ν
16
=
0.6986. As the displacement of a continuous bar is continuous, the initial displacement of other nodes should be selected according to its linear vibration mode. Figures 9, 10, and 11 plot the displacements and phase tracks of the blade tip in lead-lag and flapping directions when ν
i
n
=
12 m/s, 18 m/s, and 22 m/s, respectively. From the three figures, what can be concluded is that when ν
i
n
=
12 m/s, 18 m/s, and 22 m/s, the blade vibration is convergent, flutter, and divergent, respectively.Figure 9
The convergent displacements and phase tracks of blade tip in autonomous system.
(a)
(b)
(c)
(d)Figure 10
The flutter displacements and phase tracks of blade tip in autonomous system.
(a)
(b)
(c)
(d)Figure 11
The divergent displacements and phase tracks of blade tip in autonomous system.
(a)
(b)
(c)
(d)The blade vibration will show the convergence, flutter, and divergence characteristics with increasing wind speeds. When the operating environment is harsh, wind turbines must stop working. In order to improve the maximum operating wind speed and optimize the design of blade, dynamic adjustment of operation parameters should be employed. As the air force also has effect on the convergence, flutter, and divergence, the balance position of turbine blade vibration is not zero. Owing to the limitation of space, only the vibration simulation curves on the tip of the blade are attached in this paper (Figures9, 10, and 11). The blade deformation will decrease with the decreasing length of blades. In addition, the changing initial condition is helpful to understand the coupling vibration.


## 4.1. Stability of the Blade

After the simplification by finite difference method, if ignoring the effect of gravity, the turbine blades discussed above can be regarded as an autonomous system. The differential equations are substituted in the differential items of (1a), (1b), (2a), (2b), (4a), (4b), (5a), and (5b). According to the chosen node number, the relevant parameters are synthesized. The expressions in the autonomous system can be obtained finally in the following:(34)
x
¨
j
=
X
j
x
1
,
x
2
,
…
,
x
2
N
,
x
˙
1
,
x
˙
2
,
…
,
x
˙
2
N
,
j
=
1,2
,
…
,
2
N
,
x
=
u
1
,
…
,
u
N
;
v
1
,
…
,
v
N
.Supposingy
j
=
x
˙
j, (34) is transformed into the following:(35)
y
˙
j
=
Y
j
y
1
,
y
2
,
…
,
y
4
N
,
j
=
1,2
,
…
,
4
N
,
y
=
u
1
,
…
,
u
N
;
v
1
,
…
,
v
N
;
u
˙
1
,
…
,
u
˙
N
;
v
˙
1
,
…
,
v
˙
N
.After wind turbine blades are simplified, they become a linear time-invariant system. So the equations can easily be solved by state space expressions in (36). The blade parameters can be optimized by the linear state space. But any real system is always working under all kinds of occasional and continuous interference. After the system is subjected to interference, it can return to the stable work. This is a nonlinear stability problem and is also the focus of future research.(36)
x
˙
=
A
x
+
B
u
,
y
=
C
x
+
D
u
.If all poles in the state space fall on the left plane, the system is stable or unstable. In order to study the relationship between the variable parameters and stability, the maximal real part of pole is introduced. If the maximal real part of the pole is less than zero, the system is stable or unstable. The stable regions of the variable parameters such as the wind velocity, rotational speed, elastic modulus in lead-lag, elastic modulus in flapping, density of the blade, and the pitch angle are shown in Figure7. Taking Figure 7(a) as an example, we select the design and operating parameters of the blade according to experience and keep wind speed as the only variable and then choose different wind speeds to calculate the eigenvalues of the system. According to the criterion theorem of eigenvalues, if all the real part of eigenvalue is less than 0, the system is stable, or the system is unstable. In order to make it convenient for programming, we estimate the stability of system by using the real maximum value of the eigenvalues. The stability area of wind speed of the wind turbine can be determined in this way to provide basis for the design and usage of the wind turbine. The principle of (b) to (f) in Figure 7 is the same as that in Figure 7(a).Figure 7
Stability of variable turbine blade parameters.
(a)
(b)
(c)
(d)
(e)
(f)The flutter speeds in different rotational speeds are also derived in Figure8. Flutter is a critical state of convergence and divergence, that is, a demarcation point of stability or instability. The relationship between wind speed and rotation speed can be obtained in Figure 8. For example, when the wind speed is 20 m/s, the rotation speed of wind turbine should be less than 0.4 rad/s. Once it is larger than this value, the wind turbine will be unstable. This provides a basis for the operation of wind turbine.Figure 8
The flutter speeds in different rotational speeds.


## 4.2. Numerical Simulation

Assisted by MATLAB vibration toolboxes, the vibration on the turbine blade tip can be simulated in the nonautonomous system through time disperse. Select the number of nodesN
=
16, spin speed Ω
=
0.8 rad/s, and the initial displacements of blade tip to be μ
1
6
=
-
0.1885 and ν
16
=
0.6986. As the displacement of a continuous bar is continuous, the initial displacement of other nodes should be selected according to its linear vibration mode. Figures 9, 10, and 11 plot the displacements and phase tracks of the blade tip in lead-lag and flapping directions when ν
i
n
=
12 m/s, 18 m/s, and 22 m/s, respectively. From the three figures, what can be concluded is that when ν
i
n
=
12 m/s, 18 m/s, and 22 m/s, the blade vibration is convergent, flutter, and divergent, respectively.Figure 9
The convergent displacements and phase tracks of blade tip in autonomous system.
(a)
(b)
(c)
(d)Figure 10
The flutter displacements and phase tracks of blade tip in autonomous system.
(a)
(b)
(c)
(d)Figure 11
The divergent displacements and phase tracks of blade tip in autonomous system.
(a)
(b)
(c)
(d)The blade vibration will show the convergence, flutter, and divergence characteristics with increasing wind speeds. When the operating environment is harsh, wind turbines must stop working. In order to improve the maximum operating wind speed and optimize the design of blade, dynamic adjustment of operation parameters should be employed. As the air force also has effect on the convergence, flutter, and divergence, the balance position of turbine blade vibration is not zero. Owing to the limitation of space, only the vibration simulation curves on the tip of the blade are attached in this paper (Figures9, 10, and 11). The blade deformation will decrease with the decreasing length of blades. In addition, the changing initial condition is helpful to understand the coupling vibration.


## 5. The Nonautonomous System

When considering the effect of gravity, the equations of turbine blade discussed above can be regarded as a nonautonomous system. The differential equations are substituted in the differential items of (7a), (7b), (23a), and (23b). According to the chosen node number, the relevant parameters are synthesized. So the expressions in nonautonomous system can be expressed as follows.(37)
x
¨
j
=
X
j
x
1
,
x
2
,
…
,
x
2
N
,
x
˙
1
,
x
˙
2
,
…
,
x
˙
2
N
,
t
,
j
=
1,2
,
…
,
2
N
,
x
=
u
1
,
…
,
u
N
;
v
1
,
…
,
v
N
.Supposingy
j
=
x
˙
j, (37) is transformed into the following:(38)
y
˙
j
=
Y
j
y
1
,
y
2
,
…
,
y
4
N
,
t
,
j
=
1,2
,
…
,
4
N
,
y
=
u
1
,
…
,
u
N
;
v
1
,
…
,
v
N
;
u
˙
1
,
…
,
u
˙
N
;
v
˙
1
,
…
,
v
˙
N
.Assisted by MATLAB vibration toolboxes, the vibration on the turbine blade tip can be simulated in the nonautonomous system through time discretization. Select the number of nodesN
=
16, spin speed Ω
=
0.8 rad/s, and the initial conditions to be the same as those with autonomous system. As the displacement of a continuous bar is continuous, the initial displacement of other nodes should be selected according to the linear vibration mode. Figures 12, 13, and 14 plot the displacements and phase tracks of the blade tip in the lead-lag and flapping directions when ν
i
n
=
12 m/s, 18 m/s, and 22 m/s, respectively. It can be concluded that when ν
i
n
=
12 m/s and 22 m/s, the blade vibration is convergent and divergent and when ν
i
n
=
18 m/s, and the blade may introduce the low-speed flutter.Figure 12
The convergent displacements and phase tracks of blade tip in nonautonomous system.
(a)
(b)
(c)
(d)Figure 13
The flutter displacements and phase tracks of blade tip in nonautonomous system.
(a)
(b)
(c)
(d)Figure 14
The divergent displacements and phase tracks of blade tip in nonautonomous system.
(a)
(b)
(c)
(d)Wind speed has an important influence on the stability of the vibration of blades. No matter how the design and operation of the blade are optimized, there is always a wind speed at which the wind blade vibration is divergent. What we can do is to improve the critical wind speed of stable operation. As long as the blade design and operation are reasonable, the blade vibration can be controlled within the stability region. Compared with the autonomous system, the convergent vibration amplitude under nonautonomous system is larger. The lead-lag vibration is a cosine curve. The flapping vibration is basically a sine curve. The vibration cycle is directly related to spin speed because of the gravity effect. Due to the limited space, only the vibration curves on the tip of the blade were attached in this paper (Figures12, 13, and 14). The blade deformation will decrease with decreasing length of blades. In addition, the changing initial condition can also be helpful in understanding the coupling vibration of turbine blades.


## 6. Conclusions

According to Hamilton principle and Greenberg unsteady aerodynamic theory, partial differential equations of large horizontal axis turbine blades had been derived. This paper had used a finite difference method and state space method to deal with the general vibration systems which cannot be decoupled. In contrast with FEM, modal analysis had been studied. By the principle of stability analysis, the stable regions of the variable parameters such as the wind velocity, rotational speed, elastic modulus in lead-lag, elastic modulus in flapping, density of the blade, and the pitch angle had been obtained. The blade tip vibration in autonomous system and nonautonomous system can be simulated by MATLAB vibration toolboxes in time domain. The convergent, flutter, and divergent vibration curves had been plotted in the lead-lag and flapping direction. This paper can also be a reference to the design, manufacture, and operation control of a wind turbine.



---

*Source: 6087295-2018-04-24.xml*