Shock and Vibration
Volume 2016 (2016), Article ID 8378161, 15 pages
http://dx.doi.org/10.1155/2016/8378161

Research Article

Flap/Lag Stall Flutter Control of Large-Scale Wind Turbine Blade Based on Robust H₂ Controller

Tingrui Liu and Wei Xu

College of Mechanical & Electronic Engineering, Shandong University of Science & Technology, Qingdao 266590, China

Received 2 April 2016; Revised 13 June 2016; Accepted 22 June 2016

Academic Editor: Mickaël Lallart

Copyright © 2016 Tingrui Liu and Wei Xu. 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.

Abstract.

Flap/lag stall nonlinear flutter and active control of anisotropic composite wind turbine blade modeled as antisymmetric beam analysis have been investigated based on robust H₂ controller. The blade is modeled as single-cell thin-walled beam structure, exhibiting flap bending moment-lag transverse shear deformation, and lag bending moment-flap transverse shear deformation, with constant pitch angle set. The stall flutter control of dynamic response characteristics of composite blade incorporating nonlinear aerodynamic model is investigated based on some structural and dynamic parameters. The aeroelastic partial differential equations are reduced by Galerkin method, with the aerodynamic forces decomposed by strip theory. Robust H₂ optimal controller is developed to enhance the vibrational behavior and dynamic response to aerodynamic excitation under extreme wind conditions and stabilize structures that might be damaged in the absence of control. The effectiveness of the control algorithm is demonstrated in both amplitudes and frequencies by description of time responses, extended phase planes, and frequency spectrum analysis, respectively.

1. Introduction

Flap/lag stall flutter in dynamic stall state is an important reason of fatigue damage for large-scale wind turbine blades. How to effectively avoid such aeroelastic instability has become an important research needed to be investigated by flutter control for wind turbine blades. For the active control of aeroelastic instability for large wind turbines blade based on load reduction, some researchers have conducted a lot of work in the last few years. Cooperman et al. examined active load control concepts including trailing edge flap, flexible trailing edge, and microtab over the past decade [1]. Zhang et al. provided a numerical investigation to study the influence of various parameters of deformable trailing edge flap on the fatigue load of a large-scale wind turbine blade [2]. Ng et al. investigated the aeroelastic response and gust load alleviation of a large wind turbine using models obtained from coupling composite beams [3]. Qiao et al. [4] analyzed the ability of an adaptive wind turbine blade for vibration suppression using piezoelectric material with numerical analysis performed by finite element method modeled as box beam structure. The advanced active twist rotor blade incorporating single crystal macrofiber composite actuators and the aeroelastic analysis are designed and performed by Park and Kim [5]. Although, in these works, the active control properties dynamically represent an actual physical structure and mechanism, the intelligent control theory and algorithms are rarely used. Classical flutter and active control of single-cell thin-walled composite wind turbine blade beam with bending-twist coupling based on actual physical piezoelectric actuation and intelligent LQG algorithm are analyzed in [6]. However, stall nonlinear flutter motivated by nonlinear aerodynamic action is not mentioned. Although author’s previous work [7] investigates both stall nonlinear flutter and intelligent control, its object of analysis is only based on two-dimensional sectional study.

In present work the stall nonlinear flutter and control of large-scale wind turbine blade under extreme conditions have been investigated for composite single-cell thin-walled structure based on robust H₂ optimal control. Compared with bending-twist coupling study of wind turbine blade, research on bending-bending coupling is more important especially in some extreme wind cases. For example, in September 2013, the typhoon event in Shanwei (China), the maximum instantaneous wind speed reached 62.5 m/s, which resulted in flap/lead-lag flutter failure of a large number of wind turbine blades [8]. A set of flap/lag coupling models for composite rotating blade are considered exhibiting flap bending moment-lag transverse shear deformation and lag bending moment-flap transverse shear deformation, with constant pitch angle set. The analysis is applied to a laminated construction of antisymmetric beam. The unsteady aerodynamic is an extracted nonlinear ONERA aerodynamic model proposed by Taehyoun [7, 9]. The aerodynamics for a large-scale wind turbine blade under extreme wind conditions may involve highly separated flows. The difficulty of deep stall nonlinear aeroelastic theory and intelligent control strategies for different composite thin-walled rotating structure has far surpassed what the load reduction approaches mentioned above. Hence stall nonlinear flutter stability and robust H₂ optimal control effects based on different parameters are analyzed here by dynamic responses.

2. Aeroelastic Equations

Thin-walled composite beams with closed cross sections are widely used in the construction of wind turbine blade. Figure 1 shows the basic parameters and stall aerodynamic forces of thin-walled blade section. The sectional middle-line expression is an integration-line equation [6]. The origin of the rotating axis system is located at the rigid root. The motions of blade beam have two bending degrees of freedom: flap direction, denoted by , and lag direction perpendicular to flap, denoted by . The length of the blade is  m in direction and it has the constant rotating speed , normal to the plane of rotation. The mass per unit length is . The characteristic cross-sectional dimension, chord length, is . The thickness of blade section is denoted by , and the distance 1/4 chord length is denoted by . β represents constant pitch angle of the whole blade, the wind velocity is denoted by . Tip speed ratio of large-scale wind turbine . Nonlinear aerodynamic lift is denoted by and drag is denoted by , respectively. The composite properties of blade are as follows: the circumferentially asymmetric stiffness configuration used here consists of in the top side above the chord and in the bottom side, with ply number being 10 and ply thickness being 1.27 × 10−3 m; the maximum exterior width of section is 2.421 × 10−3 m, the maximum exterior height is  m, blade density  kg/m³, and  GPa, ,  GPa, and  GPa.

Figure 1: Coordinate system and aerodynamic forces of cross section.

The analysis is formulated to estimate the characteristics of rotating beam with elastic coupling by means of a simplified deflection theory [10]. Present study is confined to flap/lag vibrations; hence both bending-twist coupling and extension-twist are zero and corresponding displacement relations are not mentioned here. It is assumed that two moments , , axial force , and transverse shear , , including two inertial loadings and , are acting on the beam element [11]. For antisymmetric beam analysis and configuration, the flap/lag motions and force-displacement relations are obtained as follows [10, 11].

Flap motion is

Lag motion is

The relation of flap bending moment and lag transverse shear is

The relation of lag bending moment and flap transverse shear is

where the coefficients are defined in [11].

Inserting (2a) and (2b) into (1a) and (1b) and considering aerodynamic forces in directions and constant pitch angle β in Figure 1, the aeroelastic equations in and directions are obtained, respectively, as follows:

where the related parameters and can be found in Appendix A.

3. Solution Methodology

In order to solve the aeroelastic equations given by (3a) and (3b), the following steps will be implemented [6]. The first step consists of representation of displacement functions in the forms

where and are ( is the number of reserved modes) vectors of suitable shape functions; and are vectors of generalized coordinates.

Discretization of aerodynamics at the right-hand-side terms in (3a) and (3b) is by strip theory. The blade is divided into spanwise sections,. Considering constant pitch angle and ignoring elastic twist behavior, the simplified aerodynamic model extracted from [7] is applied in (3a) and (3b), which is suitable for pure pitching motion. is the circulatory part in lift equation of the simplified ONERA model; D is the drag part with a cubic form of quasistatic drag curve [9]. The aerodynamic model is [7]

with the decomposed nonlinear aerodynamic variables expressed as

where

Herein, the nonlinear items and represent the deviation of the linear static curve from the quasistatic curve for the lift and drag coefficient, respectively. The other nonlinear items , , , and can be found in [7].

Substituting (4) into (3a) and (3b) with Galerkin method applied [6] and carrying out the strip theory result in the following equations:

where the related coefficients and constants are listed in Appendix B.

Considering the two aeroelastic equations in (8a) and (8b) and the three decomposed nonlinear aerodynamic variables expression in (6a) and (6b) and assuming

result in the equations governing the whole aeroelastic system, withsubequation structures as follows:

The time response of the dynamic system of (10) is implemented in SIMULINK environment [12] with flap/lag responses obtained by integrating (4). The simulation conditions of Solver options are displayed in Table 1.

Table 1: Simulation conditions of solver options.

Items
Parameters

Solver
ode23tb (stiff/TR)-BDF2

Type
Variable-step

Relative tolerance
1 × 10⁻¹

Absolute tolerance
Auto

Shape preservation
Disable all

Initial conditions
0 in Section 5.2; 0.01 for others

4. Robust H₂ Optimal Controller

In order to control the flutter and divergent state in aeroelastic system of (10), using the method described in [6], and defining the state vector and adjoining the identity equation , (10) can be converted to the state space expression:

Since the flap/lag structures often experience time-varying dynamic stall loads, their safe and effective design requires accurate responses. The linear quadratic controller in feedback design is often used to analyze vibration control of rotating composite beam [6], which might provide a systematic way of computing the state feedback control gain matrix. In order to suppress the classical flutter and decrease the influence of measurement noise, the linear quadratic Gaussian (LQG) controller is used in [6]. However in present study, the linear quadratic controller failed to stabilize the plant or find an optimal feedback gain. The reason might lie in the fact that there are nonlinear links and items in the system as depicted in (6a) and (6b), not all unstable poles of matrix are controllable through B, or it is difficult to modify the weights in LQG controller to make related process positively definite. After all in some sense, the LQG controller in fact is man-made, with weighting matrix being.

Here the feedback design problem is cast as a robust optimal problem with H₂ controller implemented. It is often useful to have a standard problem formulation into which the particular problem might be manipulated. Such a general formulation is afforded by introduction of an augmented matrix model , which is described as [12]

The corresponding augmented state space based (11) can be expressed as

which can result in the closed-loop transfer function [13]

where is the designed feedback controller. Hence the standard H₂ optimal control problem is to find a stabilizing controller which minimizes

In addition, to simplify the formulas in H₂ algorithms, let

which can make and strictly proper. Meanwhile we can computes H₂ optimal controller for the augmented plant with suitable loop-shaping weighting functions such that H₂-norm of the closed-loop transfer function is minimized. In order to reduce the simulation time and to ensure a certain sensitivity and accuracy, sensitivity weight in Figure 2; complementary weight (not equal to zero to avoid singular problem in numerical simulation); control weight is

where the weight parameters are , and .

Figure 2: Double ports block diagram structure of H₂ control system.

As for H∞ controller in present study, when solving Riccati equations and performing H-infinity existence tests, is always not small enough, so matrix is singular to working precision during numerical simulation. It is difficult to find the optimal H∞ controller in the simulation process if structural damping [14, 15] is not estimated. This leads to the zero damping term in the original physical structure system (without action of the aerodynamic forces).

5. Numerical Simulation and Discussion

5.1. Time Responses Based on Nonzero Initial Condition

Some cases intended to highlight the effects of stall flutter suppression under extreme conditions played by H₂ optimal controller will be presented in the following. Here basic testing parameters are as follows: ply angle = 30°; constant pitch angle β = 15°. The state variables in (10) and in (11) are all in line with the nonzero initial condition being 0.01.

Figure 3 shows the uncontrolled responses of flap/lag motions under extreme wind condition of  m/s. It is noticed that for the flap displacement y, within 2 s time, the amplitude of the vibration quickly exceeded the length of the blade  m, so actually the flap motion has been in the unstable state of divergence. From the whole responses trends in the case of  m/s, the two displacements all are divergent with tiny amplitudes after 32 s or so.

Figure 3: Uncontrolled responses of flap/lag motions under extreme wind condition of  m/s.

Figure 4 shows the controlled responses of flap/lag motions under  m/s. It can be seen that, compared with uncontrolled responses, the responses controlled by H₂ controller are of more advantages from the viewpoint of amplitude suppression. It also shows the flutter amplitudes of both flap and lag displacements decrease rapidly with the change of time and tend to be steady within 8~9 s. It demonstrates obvious flutter suppression effect of H₂ controller on stall instability. Figure 5 displays the amplitudes of numerical controller for flap/lag motions. It is obvious that, even under extreme wind condition  m/s, the amplitudes of the controller will not be too large, in the range of controllable operation, and can be soon stabilized, which reflects the good control performance.

Figure 4: Controlled responses of flap/lag motions under  m/s.

Figure 5: Amplitudes of controller for flap/lag motions under  m/s.

In order to test such H₂ controller it can be commonly used at that given time and place. Another case of different wind speed is investigated. Figure 6 shows both uncontrolled and controlled responses of flap/lag motions under  m/s, respectively, with controller values attached. Note that typhoon is generally defined as the wind speed over 17 m/s. Hence the wind speed  m/s can still be considered as an undesirable wind speed situation. Almost the same conclusions of flutter control by H₂ controller can be obtained from Figure 6. In addition, it seems that the adaptive effect of the controller under more complex conditions (such as  m/s) is better than that of the others.

Figure 6: Both uncontrolled and controlled responses of flap/lag motions under  m/s, with controller values attached.

In order to verify H₂ controller law of universal, further validation concerning the cases of different structural and aerodynamic parameters is implemented.

Figure 7 demonstrates the maximum values of all the real parts of eigenvalues in the three cases: different pitch angles ranging 5°~90° at interval of 5° ( m/s and = 30°), different wind velocities ranging 5~90 m/s at interval of 5 m/s ( = 30° and β = 15°), and different ply angles ranging 5°~90° at interval of 5° ( m/s and β = 15°). Symbol “” represents the situation without control and “o” represents the situation under H₂ control at each point. In the three different cases, all the real parts of eigenvalues without control are small and fluctuate near zero. From Figure 7, it is apparent that when H₂ control is applied, almost all the real parts of eigenvalues are below zero and even some of the values have large deviation from the zero line, which indicates obvious effects on flutter control.

Figure 7: The maximum values of all the real parts of eigenvalues in different pitch angles, different wind velocities, and different ply angles.

5.2. Bounded Phase Plane Study

For the time-varying system of (10) and (11), the time domain response is often influenced by the initial condition [12]. When the initial condition is zero, the stability can be studied by the bounded phase plane [16]. Note that the phase plane here, not a description of the first-order derivative of state variable versus state variable, is an extended phase plane describes relationship between two different state variables [12]. Still take the previous case:  m/2, = 30°, and β = 15°, for example, with the initial condition of all the state variables in (10) and in (11) being zero.

Figure 8 illustrates uncontrolled responses and phase plane of flap/lag motions under extreme wind condition of  m/s with zero initial condition. From the whole response trend in the case of  m/s, the two displacements all seem to be convergent and present a state of steady convergence after 30 m/s. In fact, for the uncontrolled flap displacement , within 3 s time, the amplitude of the vibration quickly exceeds the length of the blade  m/s, so actually the flap () motion has been in a divergent state, which is in line with the previous analysis in Figure 3. In addition, as can be seen from the bounded phase plane, flap movement which is between −80 m and 80 m, far beyond the length of the blade, in fact, is an unstable state.

Figure 8: Uncontrolled responses and phase plane of flap/lag motions under extreme wind condition of  m/s with zero initial condition.

Figure 9 shows the controlled responses and phase plane of flap/lag motions under extreme wind condition of  m/s. As can be seen from the closed phase plane, flap movement is between  m to  m, with the stability greatly improved. In addition, compared with the amplitudes of controlled responses of flap/lag motions under nonzero initial condition of  m/s in Figure 4, amplitudes of controlled responses here under zero initial condition are much smaller, which indicates superior control performance under zero initial condition.

Figure 9: Controlled responses and phase plane of flap/lag motions under extreme wind condition of  m/s with zero initial condition.

Figure 10 shows amplitudes and bounded phase plane of controller for flap/lag motions under  m/s. From the view of the phase plane, the amplitudes of the controller for both flap and lag motions are reduced to a smaller range, so that it can easily be realized in theory. In addition, one obvious conclusion is that the controlled responses and phase plane in Figure 9 fluctuate almost completely centering on the similar track responses of the controller. The overall trends of the controlled object and controller are consistent.

Figure 10: Amplitudes and phase plane of controller for flap/lag motions under  m/s with zero initial condition.

5.3. Frequency Spectrum Analysis

For the time domain responses of (10) and (11), the frequency spectrum is often influenced by H₂ controller. Still take the previous case with the zero initial condition in Section 5.2:  m/2, = 30°, and β = 15°; for example, Figure 11 illustrates flap()/lag() spectrum analysis based on uncontrolled responses and controlled responses, respectively.

Figure 11: Flap/lag frequency spectrum analysis based on uncontrolled responses and controlled responses, respectively.

In view of the aforementioned time domain responses which are based on the unequal time points and cannot be directly carried out by Fourier transformation analysis, a spline interpolation of equal interval [17] based on discrete data points is implemented.

It can be seen that frequency distributions of uncontrolled flap/lag motions including the first-order and second-order frequencies are almost close to each other, except for the difference of amplitudes. Uncontrolled flap which has a larger amplitude, in turn, confirms the aforementioned instability conclusions in Figure 3. For the controlled flap and lag motions, the first-order frequencies (about 20 Hz or so) are much smaller than the uncontrolled ones (about 220 Hz); the second-order frequencies (about 220 Hz or so) are likewise much smaller than the uncontrolled ones (about 520 Hz). Reduced frequencies mean increasing cycles of tip responses; that is to say, the tip vibrations are gentle, with stability being greatly enhanced. At the same time, the greatly reduced spectral amplitudes (especially for the second-order flap frequency with very tiny amplitude) in turn confirm the aforementioned time domain response amplitudes reduction in Figure 9, which indicates the greatly improved stability.

6. Validation

The effectiveness of the stability analysis has been confirmed by many of the above methods. However, the reliability of the structural modeling and the effectiveness of the Galerkin method should be fully demonstrated. An experimental platform of cantilever beam which intended to highlight the effects of a simplified structural model from (3a) and (3b) will be presented in the following based on the first-order frequency analysis.

Consider structural equations in (3a) and (3b), let, and (3a) and (3b) can be reduced to

which implies an equivalent model for the rectangular laminated cantilever beam structure without structural damping (both bending-twist coupling and extension-twist are zero and corresponding displacement relations are not mentioned). It also means a free vibration model with no aerodynamic effect.

In order to be consistent with the experimental results as far as possible, structural damping in system (18a) and (18b) should be correctly estimated. To study the flexural vibration of the slender blade, structural damping can be expressed in terms of an equivalent viscous damping matrix, which depends on the damping coefficient , a ratio of the first-order natural frequency of the reduced structure to the corresponding first-order modal damping ratio [14]. We calculated the values of versus ply angles [14], which are listed in Table 2. Integrating structural damping into (18a) and (18b), with Galerkin method applied, (18a) and (18b) is deduced as

where

Table 2: The values of and 1st frequencies (Hz) versus ply angles (°).

Ply angles

1st frequency

0
0.0070
12.815

15
0.0135
27.106

30
0.0544
78.005

45
0.0839
98.366

60
0.0727
66.442

75
0.0523
32.636

90
0.0433
20.531

Since (19) represents a homogeneous equation system, each order frequency of this system can be directly derived from the eigenvalue analysis. The first-order frequencies versus ply angles are listed in Table 2. The specific structural parameters are the same as the experimental parameters in next paragraph.

In order to test the first-order natural frequency on a real-time basis, an experimental platform of cantilever beam with rectangular laminated cross section is built by project team members [18] in Figure 12. The composite properties of blade are as follows: the circumferentially asymmetric stiffness configuration used here consists of in the top side above the horizontal center line and in the bottom side, with ply angle being zero; the exterior width of section is 3.04 × 10−1 m and the exterior height is 1.18 × 10−1 m. The selection of width and height is characterized as far as possible to reflect the following mechanical behavior: both bending-twist coupling and extension-twist are negligible, with corresponding displacement relations not involved.

Figure 12: The experimental platform of cantilever beam.

The experimental system is devoted to the determination of the frequency response function [18]. The experimental process is briefly described as follows: the first step is to hammer the cantilever beam, and then the acceleration signal is detected by the piezoelectric sensor; the second step is to amplify the signal, remove the low frequency shaking of the clamping device through the filter, and then send it to the A/D analog module; finally after CPU treatment by PLC controller, transfer it to the Matlab environment through the Open Process Control (OPC). OPC is a series of seven specifications defined by the OPC Foundation for supporting open connectivity in industrial automation and uses Microsoft® DCOM technology to provide a communication link between OPC servers and OPC clients, detailed implementation of which can be found in [18].

The acceleration value of the first experimental process is illustrated in Figure 13. The curve of vertical flap displacement is obtained after two integral operations in Matlab environment. Finally the frequency response function can be obtained by FFT analysis. The first-order natural frequency can be affirmed by the maximum value of frequency response function, which is indicated about 12 Hz or so. Compared with 1st frequency in first line of Table 2, the measured value of 1st frequency in Figure 13 reveals an excellent agreement.

Figure 13: The acceleration value and the frequency response function of the first experimental process.

In order to test the universality, a further validation of the modeling and solution is confirmed by another experimental process in Figure 14. Almost the same experimental results can be obtained. The results displayed in Figure 14 coincide with those in Figure 13, which verify both the modeling and numerical procedures.

Figure 14: The acceleration value and the frequency response function of the second experimental process.

7. Conclusion

Flap/lag stall nonlinear flutter behavior and robust control of large-scale wind turbine blade based on H₂ controller are investigated. An analytical study devoted to modeling and simulation is presented. The related remarks can be outlined from the results.

(1)Flutter behavior under extreme conditions is investigated and discussed based on antisymmetric beam configuration with bending-bending coupling. The nonlinear aerodynamic forces are decomposed by strip theory. Flutter control is realized by robust H₂ optimal controller with the three specially customized weights. The study displays the obvious control effect of H₂ controller on flutter control.(2)The effect of H₂ controller is confirmed by time responses under nonzero initial condition, responses under zero initial condition, and closed phase plane and spectrum analysis under zero initial condition, respectively. The effectiveness of the control algorithm is demonstrated in both amplitudes and frequencies.(3)It should be stated that, due to the fact that H₂ controller is based on the full state feedback including the aerodynamic variables, the realization of this full state feedback controller is difficult in practice. However for this simplified stall aerodynamic model, the error between full state feedback and physically achievable structural feedback is relatively small [7]. So present simulation is not only the theoretical research, but also an object of reference in order to make comparison with the physically achievable feedback controllers for structural variables.

Appendix

A. The Related Parameters in (3a) and (3b)

One has

B. The Related Coefficients and Constants in (8a) and (8b)

One has

where the aerodynamic parameters not indicated can be found in [7].

Competing Interests

The authors declare that there are no competing interests regarding the publication of this paper.

Acknowledgments

This work is supported by the Natural Science Foundation of Shandong Province of China (Grant ZR2013AM016).

References

A. Cooperman and M. Martinez, “Load monitoring for active control of wind turbines,” Renewable and Sustainable Energy Reviews, vol. 41, pp. 189–201, 2015.

M. Zhang, B. Tan, and J. Xu, “Parameter study of sizing and placement of deformable trailing edge flap on blade fatigue load reduction,” Renewable Energy, vol. 77, no. 1, pp. 217–226, 2015.

B. F. Ng, R. Palacios, E. C. Kerrigan, J. M. R. Graham, and H. Hesse, “Aerodynamic load control in horizontal axis wind turbines with combined aeroelastic tailoring and trailing-edge flaps,” Wind Energy, vol. 19, no. 2, pp. 243–263, 2016.

Y. Qiao, J. Han, C. Zhang, J. Chen, and R. Gou, “Active vibration control of wind turbine blades by piezoelectric materials,” Chinese Journal of Applied Mechanics, vol. 30, no. 4, pp. 587–592, 2013.

J.-S. Park and J.-H. Kim, “Design and aeroelastic analysis of active twist rotor blades incorporating single crystal macro fiber composite actuators,” Composites Part B: Engineering, vol. 39, no. 6, pp. 1011–1025, 2008.

T. Liu, “Classical flutter and active control of wind turbine blade based on piezoelectric actuation,” Shock and Vibration, vol. 2015, Article ID 292368, 13 pages, 2015.

T. Liu, “Aeroservoelastic pitch control of stall-induced flap/lag flutter of wind turbine blade section,” Shock and Vibration, vol. 2015, Article ID 692567, 20 pages, 2015.

China Network TV Station, Guangdong Shanwei: the instantaneous wind speed over 60 meters per second, wind power generators were cut in half, 2013, http://cen.ce.cn/more/201309/24/t201309241543628.shtml

K. Taehyoun, Nonlinear large amplitude structural and aeroelastic behavior of composite rotor blades at large static deflection [Ph.D. dissertation], Massachusetts Institute of Technology, Cambridge, Mass, USA, 1992.

R. Chandra and I. Chopra, “Experimental-theoretical investigation of the vibration characteristics of rotating composite box beams,” Journal of Aircraft, vol. 29, no. 4, pp. 657–664, 1992.

E. C. Smith and I. Chopra, “Formulation and evaluation of an analytical model for composite box-beams,” Journal of the American Helicopter Society, vol. 36, no. 3, pp. 23–35, 1991.

D. Xue, Computer Aided Control Systems Design Using MATLAB Language, Tsinghua University, Beijing, China, 2nd edition, 2006.

S. Skogestad and I. Postlethwaite, Multivariable Feedback Control: Analysis and Design, John Wiley & Sons, Chichester, UK, 2nd edition, 2005.

T. Liu, “Flutter and response of composite wind turbine blade with bending-torsional coupling,” Journal of Chemical and Pharmaceutical Research, vol. 6, no. 5, pp. 1544–1550, 2014.

T.-H. Cheng, M. Ren, Z.-Z. Li, and Y.-D. Shen, “Vibration and damping analysis of composite fiber reinforced wind blade with viscoelastic damping control,” Advances in Materials Science and Engineering, vol. 2015, Article ID 146949, 6 pages, 2015.

O. Katsuhiko, Modern Control Engineering, Publishing House of Electronics Industry, Beijing, China, 5th edition, 2012.

D. Xue and Y. Chen, Advanced Applied Mathematical Problem Solutions with MATLAB, Tsinghua University, Beijing, China, 3rd edition, 2013.

Z. Yu, Numerical simulation of free vibration of piezoelectric composite cantilever blade [Ph.D. dissertation], College of Mechanical & Electronic Engineering, Shandong University of Science & Technology, Qingdao, China, 2016.