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

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 2 controller. The blade is modeled as single-cell thin-walled beam structure, exhibiting flap bendingmoment-lag transverse shear deformation, and lag bendingmoment-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 2 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.


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 singlecell 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 twodimensional 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 2 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 largescale 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 2 optimal control effects based on different parameters are analyzed here by dynamic responses.

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 middleline 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  = 30 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  = 12.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 [] 2 in the top side above the chord and [−] 2 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 0.2643 × 10 −3 m, blade density   = 1672 kg/m 3 , and  12 = 3.5 GPa, V 12 = 0.34,  11 = 25.8GPa, and  22 = 8.7 GPa.
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.

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 1× ( 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] 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].

Shock and Vibration
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, with 2 + 3 subequation 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.

Robust H 2 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 timevarying 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  LQG being 2(2 + 3) × 2(2 + 3).
Here the feedback design problem is cast as a robust optimal problem with H 2 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]  () =  11 () where () is the designed feedback controller.Hence the standard H 2 optimal control problem is to find a stabilizing controller  which minimizes In addition, to simplify the formulas in H 2 algorithms, let which can make  11 and  22 strictly proper.Meanwhile we can computes H 2 optimal controller for the augmented plant with suitable loop-shaping weighting functions such that H 2norm 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  1 = 1; complementary weight  3 = 1 − 10 (not equal to zero to avoid singular problem in numerical simulation); control weight  2 is where the weight parameters are  0 = 1.5,  0 = 3, and As for H ∞ controller in present study, when solving Riccati equations and performing H-infinity existence tests,  11 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).

Time Responses Based on Nonzero Initial Condition.
Some cases intended to highlight the effects of stall flutter suppression under extreme conditions played by H 2 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  = 65 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  = 30 m, so actually the flap motion has been in the unstable state of divergence.From the whole responses trends in the case of  = 65 m/s, the two displacements all are divergent with tiny amplitudes after 32 s or so.
Figure 4 shows the controlled responses of flap/lag motions under  = 65 m/s.It can be seen that, compared with uncontrolled responses, the responses controlled by H 2 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 2 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  = 65 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.
In order to test such H 2 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  = 15 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  = 15 m/s can still be considered as an undesirable wind speed situation.Almost the same conclusions of flutter control by H 2 controller can be obtained from Figure 6.In addition, it seems that the adaptive effect of the controller under more complex conditions (such as  = 65 m/s) is better than that of the others.
In order to verify H 2 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 ∘ ( = 65 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 ∘ ( = 65m/s and  = 15 ∘ ).Symbol " * " represents the situation without control and "o" represents the situation under H 2 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 2 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.

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:  = 65 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  = 65m/s with zero initial condition.From the whole response trend in the case of  = 65 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  = 30 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 9 shows the controlled responses and phase plane of flap/lag motions under extreme wind condition of  = 65 m/s.As can be seen from the closed phase plane, flap movement is between −13−3 m to 2−3 m, with the stability greatly improved.In addition, compared with the amplitudes of controlled responses of flap/lag motions under nonzero initial condition of  = 65 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 10 shows amplitudes and bounded phase plane of controller  for flap/lag motions under  = 65 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.

Frequency Spectrum Analysis.
For the time domain responses of ( 10) and (11), the frequency spectrum is often influenced by H 2 controller.Still take the previous case with the zero initial condition in Section 5.2:  = 65 m/2,   = 30 ∘ , and  = 15 ∘ ; for example, Figure 11 illustrates flap()/lag() 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 secondorder 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.

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  = 0,  = 0, Ω = 0, 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 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 [] 4 in the top side above the horizontal center line and [−] 4 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.
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 Microsoft5 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.
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.

Conclusion
Flap/lag stall nonlinear flutter behavior and robust control of large-scale wind turbine blade based on H 2 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.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 2 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) where the aerodynamic parameters not indicated can be found in [7].

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

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

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

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

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

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

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

Table 1 :
Simulation conditions of solver options.
The nonlinear aerodynamic forces are decomposed by strip theory.Flutter control is realized by robust H 2 optimal controller with the three specially customized weights.The study displays the obvious control effect of H 2 controller on flutter control.(2)The effect of H 2 controller is confirmed by time responses under nonzero initial condition, responses