Vibration and aeroelastic control of wind turbine blade based on B-L aerodynamic model and LQR controller

Tingrui Liu1

1College of Mechanical and Electronic Engineering, Shandong University of Science and Technology, Qingdao, China

Journal of Vibroengineering, Vol. 19, Issue 2, 2017, p. 1074-1089. https://doi.org/10.21595/jve.2016.17230
Received 2 June 2016; received in revised form 1 September 2016; accepted 19 September 2016; published 31 March 2017

Copyright © 2017 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License
Abstract.

Vibration and aeroelastic control of anisotropic composite wind turbine blade modeled as symmetric layup beam analysis have been investigated based on Beddoes-Leishman (B-L) dynamic stall aerodynamic model and linear quadratic controller. The blade is modeled as single-cell thin-walled beam structure, exhibiting flap bending-lag bending-twist coupling deformation, with constant pitch angle set. The stall flutter and aeroelastic control of composite blade are investigated based on some structural and dynamic parameters, with structural damping computed. The aeroelastic partial differential equations are reduced by Galerkin method, with the nonlinear aerodynamic forces computed by the method of fitting static elastic coefficients. Linear quadratic controller is applied to enhance the vibrational behavior in stall situation under divergent conditions and stabilize displacements that might be unstable in the absence of control.

Keywords: flap bending-lag bending-twist coupling, Beddoes-Leishman model, linear quadratic controller, Galerkin method, aeroelastic control.

1. Introduction

Flap bending-Lag bending-twist (BBT) flutter in dynamic stall state is an important reason of fatigue damage for wind turbine blades. How to effectively realize aeroelastic control has become an important research needed to be investigated in the last few years. For the active aeroelastic control methods of blades, one category is load reduction based on structural trailing edge flap, microtab and piezoelectric actuator [1-4], the other is the direct application of intelligent control theories [5]. Although in these works, the blade properties dynamically represent a real rotor blade; the analytical object is the flap-wise or torsional behavior of a helicopter blade or the vibration behavior of a wing airfoil in which stall nonlinear flutter of BBT coupling action is not mentioned.

Gebhardt et al. present an aeroelastic system of three-blade wind turbines which results from the coupling of aerodynamic model and a structural model based on a segregated formulation derived in an index-based notation [6]. However, the viscous effects of aerodynamic model in [6] can be confined to those regions close to solid surfaces, which is unable to study in depth the stall condition. Wang et al. present a nonlinear aeroelastic model for large wind turbine blades by combining blade element momentum theory and mixed-form formulation of geometrically exact beam theory [7]. However, compared with the solution of the aerodynamic analytical model (such as B-L model), the error of calculation result is still a little large. Khazar et al. present load mitigation of bending-twist coupling blade based on aeroelastic tailoring via unbalanced laminates [8]. Compared with the bending-twist coupling study of aeroelastic tailoring blade, research on BBT coupling is more general especially in some stall divergent cases. Dong et al. do predict BBT coupling blade loads based on the CFD and CSD solvers, but it is in a loosely coupled manner [9]. As for aeroelastic intelligent control, Boukhezzar et al. provide comparison between linear and nonlinear control strategies for variable speed wind turbines [10]. However, stall flutter is not involved in [10].

In present work vibration and aeroelastic control of stall nonlinear flutter of wind turbine blade with BBT coupling have been investigated for composite single-cell thin-walled structure based on B-L aerodynamic model [11-12] and linear quadratic controller (LQC). In view of the usual circumstances, there being a variable pitch motion in wind turbine, the structural equations are complicated and include nonlinear terms. The solution of stability analysis is also a complex process with nonlinear aeroelastic equations linearized. To simplify analysis, a set of BBT coupling models for composite rotating blade are considered, with constant pitch angle set. The analysis is applied to a laminated construction of symmetric layup beam structure with NACA 6315 aerofoil profile. The B-L stall aerodynamic model is applied based on the method of fitting static elastic coefficients to avoid the aerodynamic decomposition along blade span. The stall nonlinear flutter stability and LQC aeroelastic control effects based on different parameters are analyzed here by dynamic responses.

2. Equations of motions

The single-cell thin-walled composite structure is widely used in construction of wind turbine blade for theoretical study. Fig. 1 shows the blade sectional parameters and aerodynamic forces. The sectional middle-line expression of the aerofoil profile is approximately described as [13]:

(1)
y = c [ e x p ( 0.05 ( ( 1 - c o s θ φ ) + s i n θ φ ) ) + e x p ( - 0.05 ( ( 1 - c o s θ φ ) + s i n θ φ ) ) c o s θ φ ] / 4 ,
x = c [ e x p ( 0.05 ( ( 1 - c o s θ φ ) + s i n θ φ ) ) - e x p ( - 0.05 ( ( 1 - c o s θ φ ) + s i n θ φ ) ) s i n θ φ ] / 4 .

The origin of the rotating system (x, y, z) is located at the rigid root. The motions of blade beam have three degrees of freedom: flap bending direction, denoted by y, lag bending direction perpendicular to flap, denoted by x, and twist motion by θz, with constant pitch angle of the whole blade β set. The relative wind angle is denoted by ψ. The mass per unit length is m. The chord length is c= 1 m. The thickness of blade section is denoted by h, and the distance 1/4 chord length by e. Tip speed ratio λ= 6. The wind velocity is denoted by U. The length of the blade is L= 15 m in z direction and it has the constant rotating speed Ω=λU/L, normal to the plane of rotation. Nonlinear aerodynamic lift is denoted by LC, drag by D, and moment by M respectively. The composite properties of blade are: the circumferentially symmetric layup configuration used here consists of θ2n in both the top side above the chord, and the bottom side, with ply number being 10, ply thickness being 1.27×10-3m; the maximum exterior height perpendicular to chord direction is 1.506×10-1m, blade density ρ= 1672 kg/m3, and G12= 3.5 GPa, v12= 0.34, E11= 25.8 GPa, E22= 8.7 GPa.

Fig. 1. Blade section parameters and aerodynamic forces

 Blade section parameters and aerodynamic forces

Present analysis is formulated to estimate the characteristics of rotating blade with elastic coupling by means of a simplified deflection theory (SDT) [14]. The SDT method is developed for predicting the effective elastic stiffnesses and corresponding load deformation behavior of tailored composite box-beams. Based on SDT analysis for symmetric layup design and configuration, with constant pitch angle set, the BBT motions here with aerodynamic action are obtained respectively as follows:

Flap-wise motion:

(2)
m c o s β y ¨ + m s i n β x ¨ + K 55 cos β y ' ' ' ' + K 55 sin β x ' ' ' ' K 45 θ z ' ' ' '
            - 1 2 m Ω 2 L 2 - z 2 cos β y ' ' - 1 2 m Ω 2 L 2 - z 2 sin β x ' ' + m Ω 2 z L 2 - z 2 cos β y '
            + m Ω 2 z L 2 - z 2 sin β x ' = L c cos ψ + D sin ψ .

Lag-wise motion:

(3)
m c o s β x ¨ - m s i n β y ¨ + g 2 cos β x ' ' ' ' - K 66 sin β y ' ' ' ' + K 66 cos β x ' ' ' ' + K 46 θ z ' ' '
          + 1 2 m Ω 2 L 2 - z 2 sin β - m Ω 2 sin β K x y ' '
          + - 1 2 m Ω 2 L 2 - z 2 cos β + m Ω 2 cos β K x x ' '
          - m Ω 2 z sin β y ' + m Ω 2 z cos β x ' + m Ω 2 sin β y = m Ω 2 cos β x
          = D cos ψ + L c sin ( ψ )   .

Twist motion:

(4)
- m K m 2 θ ¨ z + [ K 45 c o s ( β ) - K 46 s i n ( β ) ] y ' ' ' + [ K 45 s i n ( β ) - K 46 c o s ( β ) ] x ' ' '
              + K 44 + 1 2 m Ω 2 L 2 - z 2 s i n β - m Ω 2 s i n β K x y ' ' + - 1 2 m Ω 2 K m 2 L 2 - z 2 θ z ' '
              - m Ω 2 K θ s i n ( β ) y ' + m Ω 2 K θ c o s ( β ) x ' - m Ω 2 K m 2 z θ z ' - m Ω 2 ( K m 2 2 - K m 1 2 ) θ z
                + m Ω 2 K m 2 2 - K m 1 2 β = D e s i n α + L c e c o s α + M ,

where the items of aerodynamic forces are in the right-hand-side; Kx, Kθ, and attack angle α are defined respectively as:

(5)
K x = - K 66 K 55 K 25 2 - K 22 K 55 ,           K θ = - K 46 K 55 K 25 2 - K 22 K 55 ,         α = - θ z - β + t a n - 1 U Ω z .

Herein, K22 is the chordwise shear stiffness; K25 is the spanwise bending-chordwise shear coupling stiffness; K46 is the chordwise bending-torsion coupling stiffness; K55 is the spanwise (flap) bending stiffness; K66 is the chordwise (lag) bending stiffness. The expressions of stiffness parameters Kij and the parameters of principle mass radii of gyration of blade cross section, Km1, Km2, Km, are given as:

(6)
K 22 = Q - 66 d A - Q - 26 d A 2 / Q - 22 d A ,  
K 25 = - Q - 16 y d A + Q - 12 y 2 d A Q - 26 y d A Q - 22 y 2 d A ,
K 46 = 0.245 Q - 26 x d A Q - 12 x 2 d A / Q - 22 x 2 d A ,
K 55 = Q - 11 y 2 d A - Q - 12 y 2 d A 2 / Q - 22 y 2 d A ,
K 66 = Q - 11 x 2 d A - Q - 12 y 2 d A Q - 12 x 2 d A / Q - 22 y 2 d A ,
K m 1 2 = 1 m ρ y 2 d x d y ,             K m 2 2 = 1 m ρ x 2 d x d y ,           K m 2 = K m 1 2 + K m 2 2 ,

where the elements of the ply stiffness matrix Q-ij are defined in texts discussing macro-mechanical behavior of composite piles [14].

3. B-L dynamic stall aerodynamic model based on fitting static elastic coefficients

The unsteady aerodynamic model, one example is the B-L model, is firstly used to describe the unsteady aerodynamic forces on an airfoil undergoing arbitrary motion in a flow [11, 12]. The occurrence of stall flutter of wind turbine blade leads to the need for analysis of unsteady aerodynamic forces in the stalled region. The B-L dynamic stall aerodynamic model of unsteady aerodynamics based on fitting static elastic coefficients is suggested here. The structural model in Eqs. (2-4) is coupled to the B-L model through the aerodynamic forces [11, 12]:

(7)
L c = 1 2 ρ a c V 0 2 C L d y n ,           D = 1 2 ρ a c V 0 2 C D d y n ,           M = 1 2 ρ a c 2 V 0 2 C M d y n ,

where the unsteady lift coefficient CLdyn, the unsteady drag efficient CDdyn, and the unsteady moment coefficient CMdyn are approximately expressed as:

(8)
C L d y n = C L , α 1 2 α - α 0 x a 2 + C L f s 1 2 α 1 - x a 2 + π c α ˙ 2 U ,
C D d y n = C D s t 1 2 α + 1 2 α C L d y n + C D s t 1 2 α - C D s t α 0 1 - x a 2 2 2 - 1 - f s t α 2 2 2 ,
C M d y n = C M s t 1 2 α + C L d y n a s t x a 2 - a s t f s t 1 2 α - π c α ˙ 4 U ,

were, α0 is the angle of attack at zero lift, CL,α is linear lift slope coefficient, fst is the separation point function, CLfs is the lift coefficient of fully separated flow, ast is the arm to equivalent pressure centre. The variable xa2 satisfies the following equation:

(9)
x ˙ a 2 + U 3 c x a 2 = U 3 c f s t α f ,

where αf=xa1/CL,α+α0, and the variable xa1 satisfies the equation:

(10)
x ˙ a 1 + 2 U 1.7 c x a 1 = 2 U 1.7 c C L p ,

and CLp is the unsteady lift coefficient for attached flow.

Very few applications of B-L model are taken into account spanwise effects, with B-L model originally formulated in the state-space. During the conventional aerodynamic integration process of the aerodynamics given by the right-hand-side terms in Eq. (2-4) in the application of Galerkin method, aerodynamic decomposition is often carried out in place of the integral calculus with sum operation. Hence the tremendous computing workload can embarrass the real-time display in computer environment.

In order to simplify computation, we developed an approach of fitting static elastic coefficients based on experimental data of testing modeling of NACA 6315 aerofoil to reduce order of nonlinear equation, and linearize and simplify aeroelastic stability analysis [15]. The fitted static elastic coefficients CLst, CDst, and CMst, as the 8-degree polynomial of the attack angle α, are expressed in polynomial formulas as:

(11)
P s = p 1 α 8 + p 2 α 7 + p 3 α 6 + p 4 α 5 + p 5 α 4 + p 6 α 3 + p 7 α 2 + p 8 α 1 + p 9 α 0 ,

where the related coefficients pi(i= 1, 2,…, 9) are listed in Table 1.

The nonlinear aerodynamic computation is implemented based on the three-fitted static elastic coefficients, which meanwhile can greatly reduce amount of calculation for various related parameters of the intermediate process and avoid aerodynamic strip decomposition along blade span. The three-fitted static elastic coefficients, and the other related intermediate parameters curves directly deduced by the fitted coefficients, are illustrated in Fig. 2.

Table 1. The fitted static elastic coefficients

Items
C L s t
C D s t
C M s t
P 1
0
–5.391×10-11
6.662×10-11
P 2
2.726×10-10
5.704×10-9
–7.594×10-9
P 3
–4.118×10-8
–2.245×10-7
3.309×10-7
P 4
1.923×10-6
3.864×10-6
–6.604×10-6
P 5
–2.625×10-5
–2.36×10-5
5.174×10-5
P 6
–1.989×10-4
-2.08×10-5
1.611×10-5
P 7
1.468×10-3
5.415×10-4
6.842×10-4
P 8
1.155×10-1
–1.624×10-5
–1.738×10-2
P 9
3.133×10-1
5.352×10-3
–1.294×10-1

Fig. 2. a) The curves of the three-fitted static elastic coefficients and b) the other related intermediate parameters

 a) The curves of the three-fitted static elastic coefficients and  b) the other related intermediate parameters

a)

 a) The curves of the three-fitted static elastic coefficients and  b) the other related intermediate parameters

b)

4. Solution methodology

4.1. Discretization by Galerkin method

In order to discretize the equations of motions given by Eqs. (2-3), modal approximation method based on Galerkin method [16] is implemented. The first step consists of representation of displacement functions in the forms:

(12)
y z , t = Y T z q y t ,           x z , t = X T z q x t ,           θ z z , t = S T z q θ t ,

where qy, qx and qθ are N×1 (the number of reserved modes N= 5) vectors of generalized coordinates; YTz, XTz and STz are 1×N vectors of suitable shape functions about yj, xj, sj (j= 1, 2,…, N), respectively.

Considering the boundary conditions of clamped-free flexible blade [16], the assumed mode shapes for the three displacements are the standard nonrotating, uncoupled mode shapes for a uniform cantilever structure, which are depicted as:

(13)
Y T z = y 1 , y 2 , y 3 , y 4 , , y N ,           X T z = x 1 , x 2 , x 3 , x 4 , , x N ,
S T z = s 1 , s 2 , s 3 , s 4 , , s N ,
y j z = c o s h β j z - c o s β j z + α j s i n h β j z - s i n β j z ,
x j z = c o s h β j z - c o s β j z + α j s i n h β j z - s i n β j z ,
s j z = 2 s i n γ j z ,           γ j = π j - 1 2 ,           β j = j - 1 2 π ,           α j = - c o s β j + c o s h β j s i n β j + s i n h β j .

Nonrotating mode shapes are used here because of computational ease. Since they depend only on the fixed constants αj, βj and γj, the modal integrals that result when Galerkin method is applied need to be calculated only once. Substituting Eq. (12) into Eqs. (2-4, 7) with Galerkin method applied, and assuming:

(14)
X = q y N × 1 T ,   q x N × 1 T , q θ N × 1 T T ,

result in the equations, with 3N sub-equation structures as follows:

(15)
M M X ¨ + C M X ˙ + K M X = Q M .

4.2. Assumption of structural damping

In order to accurately estimate the reliability of the blade motion, composite structural damping should be estimated. Whereas there is still no uniform structural damping estimation method. [17-20] investigates the structural damping of thin-walled composite one-cell beams based on modal damping ratio. Modal damping ratio of composite cross-section beam is defined as ratio of dissipated energy to maximum strain energy in the cycle of vibration. Based on research of [17], we investigate a structural analysis model which can determine the influences of different order modal damping factors for different structural and material parameters, and gives structural damping matrix based on Raleigh assumption [18]. To study the flexural vibration of the slender blade, based on modal damping factors above mentioned, another method is developed to compute structural damping as depicted in [19]. Structural damping from structural parts of blade section in Eq. (2-4) can be expressed in terms of an equivalent viscous damping matrix with 3×3N as:

(16a)
k 11 =   A - B 2 C d s +   B C d s 2 /   1 C d s ,
k 12 =   B C d s /   1 C d s A e ,
k 13 = -   A - B 2 C y d s -   B C d s   B C y d s /   1 C d s ,
k 14 = -   A - B 2 C x d s -   B C d s   B C x d s /   1 C d s ,
k 22 = 1 /   1 C d s A e 2 , k 23 = -   B C y d s /   1 C d s A e ,
k 24 = -   B C x d s /   1 C d s A e ,
(16b)
k 33 =   A - B 2 C y 2 d s +   B C y d s 2 /   1 C d s ,
k 34 = -   A - B 2 C x d s -   B C x d s   B C y d s /   1 C d s ,
k 44 =   A - B 2 C x 2 d s +   B C x d s 2 /   1 C d s ,
A e = 1 2   x d y d s - y d x d s d s .

Meanwhile similar to the previous variation of the strain energy formulation, the variation of the dissipated energy of the cross section in terms of the strains γ11 and γ12can be deduced as [20]:

(17)
δ W s = γ 11 γ 12 A d ( s ) B d ( s ) B d ( s ) C d ( s ) γ 11 γ 12 d s = δ Δ T C C Δ ,

where CC is 4×4 symmetric matrix with its components cijembodied. Hence the equivalent damping matrix can be similarly written as [20]:

(18)
C - = c 11 E i j c 12 F i j c 13 G i j c 14 G i j c 12 I i j c 22 J i j c 23 K i j c 24 K i j c 13 N i j c 23 O i j c 33 P i j c 34 P i j c 14 N i j c 24 O i j c 34 P i j c 44 P i j ,

where the expressions of cij are similar to those of kij, but the terms As, Bs, and Cs in kij should be replaced by the equivalent parameter terms Ads, Bds, and Cds in cij, respectively. The expressions of Ads, Bds, and Cds can be found in [20].

Damping parameter factors (as) are illustrated with different ply angles ranging from 0°-90° in Fig. 3.

Fig. 3. a s as a function of ply angles ranging from 0°-90°

as as a function of ply angles ranging from 0°-90°

Integrating structural damping Cs into Eq. (15) and carrying out the indicated variations and the required integration in Galerkin method, result in the matrix equations governing the system motion for blade tip as follows:

(19)
M M X ¨ + C S M X ˙ + K M X = Q M ,

with total 3N×3N damping matrix CSM:

(20)
C S M = C M + C s 11 y j d z C s 12 y j d z C s 13 y j d z C s 21 x j d z C s 22 x j d z C s 23 x j d z C s 31 s j d z C s 32 s j d z C s 33 s j d z 3 N × 3 N .

4.3. LQC control

In order to realize aeroelastic control for aeroelastic system Eq. (8), defining the state vector Y=XT,X˙TT and adjoining the identity equation X˙=X˙, Eq. (7) can be converted to the state space expression:

(21)
Y ˙ = A Y + B u , Y o = C Y + D u .

The LQC controller is often used to analyze vibration control of rotating composite beam, one example is linear quadratic Regulator (LQR). An advantage of LQR controller is that it provides a systematic way of computing the state feedback control gain matrix K of of the optimal control vector, so as to minimize the performance index [22]:

(22)
J = 0 ( Y T Q Y + u T R u ) d t ,

where 𝑄 and 𝑅 are the weighting matrices of system outputs and control inputs, respectively, with control signal set as:

(23)
u * = - K Y t = - R - 1 B T P Y t .

Herein, P is a real symmetric matrix, satisfies the Riccati equation:

(24)
P A + A T P - P B R - 1 B T P + Q = 0 .

The controlled state equation of the closed loop system [23] with state feedback is obtained as:

(25)
Y ˙ = A - B K Y + B u , Y o = C - D K Y + D u .

In order to suppress the too large initial vibration amplitude in LQR process, another LQC controller, the conventional linear quadratic Gaussian (LQG) controller is used to stabilize the displacements [24]. It consists of first determining the optimal control K to a deterministic LQR problem, the second step to find an optimal estimate (of the state Y) given by a Kalman filter and independent of Q and R, and the final connection of the Kalman filter and linear quadratic optimal gain K. However, the conventional LQG algorithm based on directly optimal control does not suppress the amplitude of the controller itself.

Here in order to suppress the too large amplitude of controller K, a special LQG controller with integral action for inverse response process is designed to decrease the amplitude of controller response [25]. LQG controller with integral action and reference input is illustrated in Fig. 4. The standard LQG design does not give a controller with integral action, so the plant here is augmented with an integrator before designing the state feedback regulator. A LQG servo-controller that uses the setpoint command r and measurements Y0to generate the control signal u has integral action to ensure that the output Y0 tracks the command r with smaller controller amplitude. Here the control error r-Y0 is integrated and the regulator K is designed for the plant augmented with integrated states.

The Kalman filter has the structure of an ordinary state-estimator or observer characterized by:

(26)
x ^ ` 0 = A x ^ 0 + B u + k f Y 0 - C x ^ 0 ,

where x0 is the state variable in LQG control with its optimal estimate being x^`0; the optimal choice of Kf minimizes E{[x0-x^0]T[x0-x^0]}, which is independent of constant power spectral density matrices W and V in Gaussian stochastic noise processes and given by:

(27)
K f = Y x C T V - 1 ,

where Yx is the unique positive-semidefinite solution of the algebraic Riccati equation:

(28)
Y x A T + A Y x - Y x C T V - 1 C Y x + W = 0 .

Fig. 4. LQG controller with integral action and reference input

 LQG controller with integral action and reference input

5. Results and discussions

Aeroelastic stability analysis can be implemented by time domain responses of Eqs. (9-10). Influences of ply angles and pitch angles on aeroelastic stability are investigated. The obvious effects of LQR and LQG controllers are demonstrated by vibration suppression for divergent cases.

5.1. Influences of ply angles and pitch angles on aeroelastic stability

Generally, with the range of 30°-60° of both ply angles and pitch angles, good aeroelastic stability of wind turbine blade can be obtained [18, 26]. Some cases intended to highlight the effects of ply angles and pitch angles under condition of U= 5 m/s are presented.

Fig. 5-6 show the uncontrolled blade tip responses of flap bending yL,t, lag bending xL,t, and twist θzL,t, for ply angles (θ) and pitch angles (β) in the range of 0°-90°, respectively. It is noticed from Fig. 5 that for all the three displacements within 40 s time, the amplitude of the vibration of θ= 30° quickly decreases in the state of convergence, and stabilizes at a lower numerical point. Comprehensive evaluation being from Fig. 6, it can be seen that the performance of β= 45° is in the middle of all pitch angle instances, with all the three displacements convergent, and with moderate amplitudes after 25 s or so. Hence the follow-up study for stability analysis and control of divergent cases will take the fixed angles of θ= 30° and β= 45°.

Fig. 5. The responses of bending (y), bending (x), and twist (θz), for ply angles in the range of 0°-90°

 The responses of bending (y), bending (x), and twist (θz), for ply angles in the range of 0°-90°  The responses of bending (y), bending (x), and twist (θz), for ply angles in the range of 0°-90°  The responses of bending (y), bending (x), and twist (θz), for ply angles in the range of 0°-90°

Fig. 6. The responses of bending (y), bending (x), and twist (θz), for pitch angles in the range of 0°-90°

 The responses of bending (y), bending (x), and twist (θz), for pitch angles in the range of 0°-90°  The responses of bending (y), bending (x), and twist (θz), for pitch angles in the range of 0°-90°  The responses of bending (y), bending (x), and twist (θz), for pitch angles in the range of 0°-90°

5.2. Influences of LQR control on aeroelastic stability

5.2.1. Influences of LQR control based on divergent cases

Fig. 7 shows the controlled responses of the three motions by LQR controller under U= 10 m/s based on basic angles setting mentioned above. It shows that the three uncontrolled displacements, yL,t, xL,t, and θzL,t, all show divergent states, with the amplitude of the vibration quickly exceeding the length of the blade L= 15 m within time 15 s. Also in contrast with uncontrolled cases, the displacements controlled by LQR controller are of more advantages from the viewpoint of amplitude suppression. It can also be seen that the controlled flutter amplitudes of all the three motions decrease rapidly with the change of time, and tend to be steady within 15 s with tiny displacement deflection. It demonstrates apparent aeroelastic control effect of LQR controller on divergent aeroelastc instability.

In order to test such a LQR controller can be commonly used. Another divergent case of U= 15 m/s is investigated in Fig. 8. Compared with the uncontrolled case of U= 10 m/s, uncontrolled displacements especially for flap and lag motions under U= 15 m/s, have a larger divergence, with the vibration amplitude more than the length of the blade within 3 s time. Also, Fig. 8 displays the amplitudes of LQR control for the three motions. It is obviously demonstrated the same conclusion as Fig. 7, namely that, the influences of LQR control on divergent cases are apparent.

Fig. 7. Controlled responses of the three motions under U= 10 m/s

 Controlled responses of the three motions under U= 10 m/s  Controlled responses of the three motions under U= 10 m/s  Controlled responses of the three motions under U= 10 m/s

5.2.2. Analytical proof of the stability of LQR controller

The effectiveness of the stability analysis by LQR controller can be confirmed by eigenvalue analysis for homogeneous equation systems (only concerning the characteristic matrices of A and A-BK in Eqs. (21, 25), respectively) defined in state-space for U= 10 m/s and U= 15 m/s, respectively. Table 2-3 show comparisons of real parts of eigenvalues under uncontrolled cases and those controlled by LQR for U= 10 m/s and U= 15 m/s, respectively. If real parts of all the eigenvalues of the aeroelastic system are negative, the system is significantly stable. Furthermore, if the real parts of eigenvalues are arranged in order from large to small and named in order as Eig1, Eig2, Eig3, etc. (Note that the real parts are characterized by a total of thirty data, with only the first eight data listed in Table), the maximum aeroelastic instability might be decided. So, stability rules influenced by LQR controller might be clarified by comparison of corresponding eigenvalue pairs.

It can be seen from Tables 2-3 that the homogeneous equation system obtained by control of LQR controller is stable, with all the Real-parts of eigenvalues being less than zeros. Compared with the controlled system, the Real-parts of some characteristic values of the uncontrolled system are greater than zeros, which means that the uncontrolled system is not stable. The maximum value of Real-parts of uncontrolled system in Table 2 is 0.4278, with the maximum value of Real-parts of controlled system being –0.0285. We can get very similar stability conclusions in Table3 with U= 15 m/s. The maximum value of Real-parts of uncontrolled system in Table3 is 2.7322, with the maximum value of Real-parts of controlled system being –0.1403. Moreover, for the comparison of every corresponding eigenvalue pair, the value under controlled situation is always less than that under uncontrolled situation, which means the excellent stability of the proposed LQR controller.

The stability conclusions from eigenvalue analysis in Tables 2-3 are in qualitative agreement with those from time responses in Figs. 7-8. It should be stated that the systems of Eqs. (21, 25) all are forced systems; compared with the stability conclusions from responses analysis of the forced systems, the stability conclusions from eigenvalues analysis of the homogeneous systems are more conservative.

Fig. 8. Controlled responses of the three motions under U= 15 m/s

 Controlled responses of the three motions under U= 15 m/s  Controlled responses of the three motions under U= 15 m/s  Controlled responses of the three motions under U= 15 m/s

Table 2. Comparisons of real parts under uncontrolled cases and those controlled by LQR for U= 10 m/s

U =   10m/s
Eig1
Eig2
Eig3
Eig4
Eig5
Eig6
Eig7
Eig8
Uncontrolled
0.4278
–0.0272
–0.1397
–0.1402
–0.1405
–0.1407
–0.1415
–0.1454
Controlled
–0.0285
–0.1397
–0.1402
–0.1405
–0.1407
–0.1416
–0.1451
–0.1484

Table 3. Comparisons of real parts under uncontrolled cases and those controlled by LQR for U= 15 m/s

U =   15 m/s
Eig1
Eig2
Eig3
Eig4
Eig5
Eig6
Eig7
Eig8
Uncontrolled
2.7322
0.4279
–0.1403
–0.1414
–0.1421
–0.1425
–0.1445
–0.1532
Controlled
–0.1403
–0.1414
–0.1421
–0.1426
–0.1447
–0.1525
–0.1604
–0.1853

5.3. Influences of LQG control on divergent cases based on integral action for inverse response process

In view of the real-time implementation of LQR control, controller amplitudes of state Y should be kept within a certain range to avoid hardware damage problems. Fig. 9(a) shows the constant LQR controller amplitudes of 30 (6N) state variables Y in the divergent case of U= 15 m/s in section 5.2.1. Notice that the horizontal coordinate values are represented by 30 state variables in turn. The order of magnitude of controller amplitudes reaches an incredible degree, so it is difficult to achieve in controller hardware design. Especially for seventh to tenth state variables, which imply the displacements of lag motion, the huge amplitude values also indicate that it is difficult to realize the flutter suppression for lag displacements.

Fig. 9(b) depicts LQG controller amplitudes in the divergent case of U= 15 m/s as a function of time. For the three different motions, the amplitude values of the three controllers all are small, each stable at a very small value, which indicates the feasibility of real-time operation.

Fig. 9. a) The constant LQR controller amplitudes of 30 (6N) state variables Y, b) and LQG controller amplitudes as function of time, in the divergent case of U= 15 m/s, respectively

 a) The constant LQR controller amplitudes of 30 (6N) state variables Y, b) and LQG controller amplitudes as function of time, in the divergent case of U= 15 m/s, respectively

a) LQR

 a) The constant LQR controller amplitudes of 30 (6N) state variables Y, b) and LQG controller amplitudes as function of time, in the divergent case of U= 15 m/s, respectively

b) LQG

In Fig. 10(a-c), the control effects of both LQR and LQG in the divergent case of U= 15 m/s are involved, where the comparison of the three controlled displacements is indicated. In contradistinction with the trends occurring in LQR control, the results of LQG are stable in larger numerical values with longer steady state time, which means that the reduction of the amplitude of the controller is accompanied by the attenuation of the control performance. In Fig. 10(d), the real parts deviations of the first eight largest eigenvalues (similar to Tables 2-3) between LQR and LQG are illustrated. The very small deviation data show relatively close control performance between LQR and LQG, especially the first four decisive values (Eig4-Eig1) are very close, which demonstrates the stability rules coincide with those in Fig. 10(a-c) obtained via the numerical responses of the three displacements.

Anyway, Fig. 9-10 show an expected conclusion, namely that for stability analysis of stall flutter, LQR and LQG, each can play an important role in enhancing the vibration behavior of the structure.

Fig. 10. Effects of both LQR and LQG in the divergent case of U= 15 m/s

 Effects of both LQR and LQG in the divergent case of U= 15 m/s

a) Flap

 Effects of both LQR and LQG in the divergent case of U= 15 m/s

b) Lag

 Effects of both LQR and LQG in the divergent case of U= 15 m/s

c) Twist

 Effects of both LQR and LQG in the divergent case of U= 15 m/s

d) Real parts deviations

6. Conclusions

Bending-bending-twist stall flutter behavior and linear quadratic control of wind turbine blade are investigated. The analytical study and numerical illustrations reveal the results:

1) Vibration behavior is investigated and discussed based on symmetric layup configuration. The nonlinear aerodynamic forces are characterized by B-L model based on fitted static aeroelastic coefficients. Stall flutter and aeroelastic stability are analyzed based on ply angles and pitch angles. Aeroelastic control is carried out by LQR and LQG controllers.

2) The method of fitting aeroelastic coefficients in B-L model can greatly reduce the complexity of processing, avoid the nonlinear aerodynamic strip decomposition and reduce the number of equations, and realize Galerkin method directly through the integral action along the blade span.

3) The obvious effect of LQR controller is demonstrated by vibration suppression for divergent cases, with validation implemented by eigenvalue analysis. In order to avoid excessive controller amplitude in LQR, LQG control in applied. However, the reduction of the controller amplitude is accompanied by the attenuation of the control performance in LQG, with larger time response amplitude and longer steady state time obtained.

4) Control performance of LQG based on integral action for inverse response process is illustrated, accompanied by excellent real-time controller amplitude. In contrast to responses in the divergent cases, the effectiveness of LQR and LQG is reflected in both decreasing amplitudes and reducing steady state time.

Acknowledgements

This work is supported by the National Natural Science Foundation of China (Grant No. 51675315), and the Natural Science Foundation of Shandong Province of China (Grant No. ZR2013AM016).

References

  1. Zhang M., Tan B., Xu J. Parameter study of sizing and placement of deformable trailing edge flap on blade fatigue load reduction. Renewable Energy, Vol. 77, 2015, p. 217-226. [Publisher]
  2. Ng B. F., Palacios R., Kerrigan E. C., et al. Aerodynamic load control in HAWT with combined aeroelastic tailoring and trailing-edge flaps. Wind Energy, Vol. 19, Issue 2, 2015, p. 243-263. [CrossRef]
  3. Li N., Mark J. B., Yang H., et al. Numerical investigation of flapwise-torsional vibration model of a smart section blade with microtab. Shock and Vibration, Vol. 2015, 2015, p. 136026. [Publisher]
  4. Qiao Y., Han J., and Zhang C., et al. Active vibration control of wind turbine blades by piezoelectric materials. Chinese Journal of Applied Mechanics, Vol. 30, Issue 4, 2013, p. 587-592. [CrossRef]
  5. Qian W., Huang R., Hu H., et al. Active flutter suppression of a multiple-actuated-wing wind tunnel model. Chinese Journal of Aeronautics, Vol. 27, Issue 6, 2014, p. 1451-1460. [Publisher]
  6. Gebhardt C. G., Roccia B. A. Non-linear aeroelasticity: an approach to compute the response of three-blade large-scale horizontal-axis wind turbines. Renewable Energy, Vol. 66, 2014, p. 495-514. [Publisher]
  7. Wang L., Liu X., Nathalie R., et al. Nonlinear aeroelastic modelling for wind turbine blades based on blade element momentum theory and geometrically exact beam theory. Energy, Vol. 76, 2014, p. 487-501. [Publisher]
  8. Khazar H., Sung K. H. Load mitigation of wind turbine blade by aeroelastic tailoring via unbalanced laminates composites. Composite Structures, Vol. 128, 2015, p. 122-133. [Publisher]
  9. Dong O. Y., Kwon O. J. Predicting wind turbine blade loads and aeroelastic response using a coupled CFD-CSD method. Control Engineering Practice, Renewable Energy, Vol. 70, 2014, p. 184-196. [Publisher]
  10. Boukhezzar B., Siguerdidjane H. Comparison between linear and nonlinear control strategies for variable speed wind turbines. Control Engineering Practice, Vol. 18, 2010, p. 1357-1368. [Publisher]
  11. Hansen M. H., Gaunaa M., Madsen H. A. A Beddoes-Leishman Type Dynamic Stall Model in State-Space and Indicial Formulation. Risø-R-1354, Risø, Roskilde, Denmark, 2004. [CrossRef]
  12. Kallesøe B. S. A low-order model for analysing effects of blade fatigue load control. Wind Energy, Vol. 9, 2006, p. 421-436. [Publisher]
  13. Wu X., Jin C., Wenzhong S. Integration study on airfoil profile for wind turbines. China Mechanical Engineering, Vol. 20, Issue 2, 2009, p. 211-213. [CrossRef]
  14. Smith E. C., Chopra I. Formulation and evaluation of an analytical model for composite box-beams. Journal of the American Helicopter Society, Vol. 36, Issue 3, 1991, p. 23-35. [Publisher]
  15. Liu T., Ren Y., Yang X. Aeroelastic stability analysis based on fitted aerodynamic coefficients. ACTA Energiae Solaris Sinica, Vol. 31, Issue 4, 2010, p. 513-516. [CrossRef]
  16. Hodges D. H., Pierce G. A. Introduction to Structural Dynamics and Aeroelasticity. Cambrige University Press, Cambrige, 2002. [Publisher]
  17. Ren Y., Du X., Sun S., et al. Structural damping of thin-walled composite one-cell beams. Journal of Vibration and Shock, Vol. 36, Issue 3, 2012, p. 141-152. [CrossRef]
  18. Liu T. Theoretical Modeling and Numerical Analysis for Dynamic Stall Aeroelastic Mechanics of Large-Scale Wind Turbine Blade. Ph.M. dissertation, College of Mechanical and Electronic Engineering, Qingdao, China, Shandong University of Science and Technology, 2011. [CrossRef]
  19. Kumar L. R., Datta P. K., Prabhakara D. L. Dynamic instability characteristics of laminated composite plates subjected to partial follower edge load with damping. International Journal of Mechanical Sciences, Vol. 45, 2003, p. 1429-1448. [Publisher]
  20. Ren Y., Zhang X., Liu Y., et al. Vibration and instability of rotating composite thin-walled shafts with internal damping. Shock and Vibration, Vol. 2014, 2014, p. 123271. [CrossRef]
  21. Bran A., Armanios A. B. Free vibration analysis of anisotropic thin-walled closed-section beams. AIAA Journal, Vol. 10, 1995, p. 1905-1910. [CrossRef]
  22. Ogata K. Modern Control Engineering. Fifth Edition, Publishing House of Electronics Industry, Beijing, China, 2012. [CrossRef]
  23. Xue D. Computer Aided Control Systems Design Using MATLAB Language. Second Edition, Tsinghua University Publishing Company, Beijing, 2006. [CrossRef]
  24. Liu T. Classical flutter and active control of wind turbine blade based on piezoelectric actuation. Shock and Vibration, Vol. 2015, 2015, p. 292368. [Publisher]
  25. Skogestad S., Postlethwaite I. Multivariable Feedback Control: Analysis and Design. Second Edition, Wiley, Chichester, UK, 2005. [CrossRef]
  26. Liu T. Stall flutter suppression for absolutely divergent motions of wind turbine blade base on H-infinity mixed-sensitivity synthesis method. The Open Mechanical Engineering Journal, Vol. 9, Issue 8, 2015, p. 752-760. [Publisher]