Vibration and aeroelastic control of wind turbine blade based on BL aerodynamic model and LQR controller
Tingrui Liu^{1}
^{1}College of Mechanical and Electronic Engineering, Shandong University of Science and Technology, Qingdao, China
Journal of Vibroengineering, Vol. 19, Issue 2, 2017, p. 10741089.
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
JVE Conferences
Vibration and aeroelastic control of anisotropic composite wind turbine blade modeled as symmetric layup beam analysis have been investigated based on BeddoesLeishman (BL) dynamic stall aerodynamic model and linear quadratic controller. The blade is modeled as singlecell thinwalled beam structure, exhibiting flap bendinglag bendingtwist 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 bendinglag bendingtwist coupling, BeddoesLeishman model, linear quadratic controller, Galerkin method, aeroelastic control.
1. Introduction
Flap bendingLag bendingtwist (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 [14], 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 flapwise 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 threeblade wind turbines which results from the coupling of aerodynamic model and a structural model based on a segregated formulation derived in an indexbased 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 mixedform formulation of geometrically exact beam theory [7]. However, compared with the solution of the aerodynamic analytical model (such as BL model), the error of calculation result is still a little large. Khazar et al. present load mitigation of bendingtwist coupling blade based on aeroelastic tailoring via unbalanced laminates [8]. Compared with the bendingtwist 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 singlecell thinwalled structure based on BL aerodynamic model [1112] 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 BL 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 singlecell thinwalled 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 middleline expression of the aerofoil profile is approximately described as [13]:
$x=c\left[\mathrm{e}\mathrm{x}\mathrm{p}\right(0.05\left(\right(1\mathrm{c}\mathrm{o}\mathrm{s}{\theta}_{\phi})+\mathrm{s}\mathrm{i}\mathrm{n}{\theta}_{\phi}))\mathrm{e}\mathrm{x}\mathrm{p}(0.05\left(\right(1\mathrm{c}\mathrm{o}\mathrm{s}{\theta}_{\phi})+\mathrm{s}\mathrm{i}\mathrm{n}{\theta}_{\phi})\left)\mathrm{s}\mathrm{i}\mathrm{n}{\theta}_{\phi}\right]/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 ${\theta}_{z}$, with constant pitch angle of the whole blade $\beta $ set. The relative wind angle is denoted by $\psi $. 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 $\lambda =$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 $\mathrm{\Omega}=\lambda U/L$, normal to the plane of rotation. Nonlinear aerodynamic lift is denoted by ${L}_{C}$, drag by $D$, and moment by $M$ respectively. The composite properties of blade are: the circumferentially symmetric layup configuration used here consists of ${\left[\theta \right]}_{2n}$ in both the top side above the chord, and the bottom side, with ply number being 10, ply thickness being 1.27×10^{3}^{}m; the maximum exterior height perpendicular to chord direction is 1.506×10^{1}^{}m, blade density $\rho =$1672 kg/m^{3}, and ${G}_{12}=$3.5 GPa, ${v}_{12}=$0.34, ${E}_{11}=$ 25.8 GPa, ${E}_{22}=$8.7 GPa.
Fig. 1. 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 boxbeams. 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:
Flapwise motion:
$\frac{1}{2}m{\mathrm{\Omega}}^{2}\left({L}^{2}{z}^{2}\right)\mathrm{cos}\left(\beta \right){y}^{\text{'}\text{'}}\frac{1}{2}m{\mathrm{\Omega}}^{2}\left({L}^{2}{z}^{2}\right)\mathrm{sin}\left(\beta \right){x}^{\text{'}\text{'}}+m{\mathrm{\Omega}}^{2}z\left({L}^{2}{z}^{2}\right)\mathrm{cos}\left(\beta \right){y}^{\text{'}}$
$+m{\mathrm{\Omega}}^{2}z\left({L}^{2}{z}^{2}\right)\mathrm{sin}\left(\beta \right){x}^{\text{'}}={L}_{c}\mathrm{cos}\left(\psi \right)+D\mathrm{sin}\left(\psi \right).$
Lagwise motion:
$+\left[\frac{1}{2}m{\mathrm{\Omega}}^{2}\left({L}^{2}{z}^{2}\right)\mathrm{sin}\left(\beta \right)m{\mathrm{\Omega}}^{2}\mathrm{sin}\left(\beta \right){K}_{x}\right]{y}^{\text{'}\text{'}}$
$+\left[\frac{1}{2}m{\mathrm{\Omega}}^{2}\left({L}^{2}{z}^{2}\right)\mathrm{cos}\left(\beta \right)+m{\mathrm{\Omega}}^{2}\mathrm{cos}\left(\beta \right){K}_{x}\right]{x}^{\text{'}\text{'}}$
$m{\mathrm{\Omega}}^{2}z\mathrm{sin}\left(\beta \right){y}^{\text{'}}+m{\mathrm{\Omega}}^{2}z\mathrm{cos}\left(\beta \right){x}^{\text{'}}+m{\mathrm{\Omega}}^{2}\mathrm{sin}\left(\beta \right)y=m{\mathrm{\Omega}}^{2}\mathrm{cos}\left(\beta \right)x$
$=D\mathrm{cos}\left(\psi \right)+{L}_{c}\mathrm{sin}(\psi ).$
Twist motion:
$+\left[{K}_{44}+\frac{1}{2}m{\mathrm{\Omega}}^{2}\left({L}^{2}{z}^{2}\right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\beta \right)m{\mathrm{\Omega}}^{2}\mathrm{s}\mathrm{i}\mathrm{n}\left(\beta \right){K}_{x}\right]y\text{'}\text{'}+\left[\frac{1}{2}m{\mathrm{\Omega}}^{2}{K}_{m}^{2}\left({L}^{2}{z}^{2}\right)\right]{\theta}_{z}^{\text{'}\text{'}}$
$m{\mathrm{\Omega}}^{2}{K}_{\theta}\mathrm{s}\mathrm{i}\mathrm{n}\left(\beta \right)y\text{'}+m{\mathrm{\Omega}}^{2}{K}_{\theta}\mathrm{c}\mathrm{o}\mathrm{s}\left(\beta \right)x\text{'}m{\mathrm{\Omega}}^{2}{K}_{m}^{2}z{\theta}_{z}^{\text{'}}m{\mathrm{\Omega}}^{2}({K}_{m2}^{2}{K}_{m1}^{2}){\theta}_{z}$
$+m{\mathrm{\Omega}}^{2}\left({K}_{m2}^{2}{K}_{m1}^{2}\right)\beta =De\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\alpha}\right)+{L}_{c}e\mathrm{c}\mathrm{o}\mathrm{s}\left(\alpha \right)+M,$
where the items of aerodynamic forces are in the righthandside; ${K}_{x}$, ${K}_{\theta}$, and attack angle $\alpha $ are defined respectively as:
Herein, ${K}_{22}$ is the chordwise shear stiffness; ${K}_{25}$ is the spanwise bendingchordwise shear coupling stiffness; ${K}_{46}$ is the chordwise bendingtorsion coupling stiffness; ${K}_{55}$ is the spanwise (flap) bending stiffness; ${K}_{66}$ is the chordwise (lag) bending stiffness. The expressions of stiffness parameters ${K}_{ij}$ and the parameters of principle mass radii of gyration of blade cross section, ${K}_{m1}$, ${K}_{m2}$, ${K}_{m}$, are given as:
${K}_{25}=\iint {\stackrel{}{Q}}_{16}ydA+\iint \frac{{\stackrel{}{Q}}_{12}{y}^{2}dA\iint {\stackrel{}{Q}}_{26}ydA}{\iint {\stackrel{}{Q}}_{22}{y}^{2}dA},$
${K}_{46}=0.245\iint {\stackrel{}{Q}}_{26}xdA\iint {\stackrel{}{Q}}_{12}{x}^{2}dA/\iint {\stackrel{}{Q}}_{22}{x}^{2}dA,$
${K}_{55}=\iint {\stackrel{}{Q}}_{11}{y}^{2}dA{\left(\iint {\stackrel{}{Q}}_{12}{y}^{2}dA\right)}^{2}/\iint {\stackrel{}{Q}}_{22}{y}^{2}dA,$
${K}_{66}=\iint {\stackrel{}{Q}}_{11}{x}^{2}dA\iint {\stackrel{}{Q}}_{12}{y}^{2}dA\iint {\stackrel{}{Q}}_{12}{x}^{2}dA/\iint {\stackrel{}{Q}}_{22}{y}^{2}dA,$
${K}_{m1}^{2}=\frac{1}{m}\iint \rho {y}^{2}dxdy,{K}_{m2}^{2}=\frac{1}{m}\iint \rho {x}^{2}dxdy,{K}_{m}^{2}={K}_{m1}^{2}+{K}_{m2}^{2},$
where the elements of the ply stiffness matrix ${\stackrel{}{Q}}_{ij}$ are defined in texts discussing macromechanical behavior of composite piles [14].
3. BL dynamic stall aerodynamic model based on fitting static elastic coefficients
The unsteady aerodynamic model, one example is the BL 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 BL dynamic stall aerodynamic model of unsteady aerodynamics based on fitting static elastic coefficients is suggested here. The structural model in Eqs. (24) is coupled to the BL model through the aerodynamic forces [11, 12]:
where the unsteady lift coefficient ${C}_{L}^{dyn}$, the unsteady drag efficient ${C}_{D}^{dyn}$, and the unsteady moment coefficient ${C}_{M}^{dyn}$ are approximately expressed as:
${C}_{D}^{dyn}={C}_{D}^{st}\left(\frac{1}{2}\alpha \right)+\frac{1}{2}\alpha {C}_{L}^{dyn}+\left[{C}_{D}^{st}\left(\frac{1}{2}\alpha \right){C}_{D}^{st}\left({\alpha}_{0}\right)\right]\left[{\left(\frac{1\sqrt{{x}_{a2}}}{2}\right)}^{2}\left.{\left(\frac{1\sqrt{{f}^{st}\left(\frac{\alpha}{2}\right)}}{2}\right)}^{2}\right]\right.,$
${C}_{M}^{dyn}={C}_{M}^{st}\left(\frac{1}{2}\alpha \right)+{C}_{L}^{dyn}\left[{a}^{st}\left({x}_{a2}\right){a}^{st}\left({f}^{st}\left(\frac{1}{2}\alpha \right)\right)\right]\frac{\pi c\dot{\alpha}}{4U},$
were, ${\alpha}_{0}$ is the angle of attack at zero lift, ${C}_{L,\alpha}$ is linear lift slope coefficient, ${f}^{st}$ is the separation point function, ${C}_{L}^{fs}$ is the lift coefficient of fully separated flow, ${a}^{st}$ is the arm to equivalent pressure centre. The variable ${x}_{a2}$ satisfies the following equation:
where ${\alpha}_{f}={x}_{a1}/{C}_{L,\alpha}+{\alpha}_{0}$, and the variable ${x}_{a1}$ satisfies the equation:
and ${C}_{L}^{p}$ is the unsteady lift coefficient for attached flow.
Very few applications of BL model are taken into account spanwise effects, with BL model originally formulated in the statespace. During the conventional aerodynamic integration process of the aerodynamics given by the righthandside terms in Eq. (24) 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 realtime 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 ${C}_{L}^{st}$, ${C}_{D}^{st}$, and ${C}_{M}^{st}$, as the 8degree polynomial of the attack angle $\alpha $, are expressed in polynomial formulas as:
where the related coefficients ${p}_{i}$_{}($i=$ 1, 2,…, 9) are listed in Table 1.
The nonlinear aerodynamic computation is implemented based on the threefitted 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 threefitted 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}^{st}$

${C}_{D}^{st}$

${C}_{M}^{st}$

${P}_{1}$

0

–5.391×10^{}^{1}^{1}

6.662×10^{}^{1}^{1}

${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 threefitted static elastic coefficients and b) the other related intermediate parameters
a)
b)
4. Solution methodology
4.1. Discretization by Galerkin method
In order to discretize the equations of motions given by Eqs. (23), modal approximation method based on Galerkin method [16] is implemented. The first step consists of representation of displacement functions in the forms:
where ${q}_{y}$, ${q}_{x}$ and ${q}_{\theta}$ are $N\times 1$ (the number of reserved modes $N=$5) vectors of generalized coordinates; ${Y}^{T}\left(z\right)$, ${X}^{T}\left(z\right)$ and ${S}^{T}\left(z\right)$ are $1\times N$ vectors of suitable shape functions about ${y}_{j}$, ${x}_{j}$, ${s}_{j}$ ($j=$1, 2,…, $N$), respectively.
Considering the boundary conditions of clampedfree 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:
${S}^{T}\left(z\right)=\left[{s}_{1},{s}_{2},{s}_{3},{s}_{4},\cdots ,{s}_{N}\right],$
${y}_{j}\left(z\right)=\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\left({\beta}_{j}z\right)\mathrm{c}\mathrm{o}\mathrm{s}\left({\beta}_{j}z\right)+{\alpha}_{j}\left(\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}\left({\beta}_{j}z\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\beta}_{j}z\right)\right),$
${x}_{j}\left(z\right)=\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\left({\beta}_{j}z\right)\mathrm{c}\mathrm{o}\mathrm{s}\left({\beta}_{j}z\right)+{\alpha}_{j}\left(\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}\left({\beta}_{j}z\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\beta}_{j}z\right)\right),$
${s}_{j}\left(z\right)=\sqrt{2}\mathrm{s}\mathrm{i}\mathrm{n}\left({\gamma}_{j}z\right),{\gamma}_{j}=\pi \left(j\frac{1}{2}\right),{\beta}_{j}=\left(j\frac{1}{2}\right)\pi ,{\alpha}_{j}=\frac{\mathrm{c}\mathrm{o}\mathrm{s}{\beta}_{j}+\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}{\beta}_{j}}{\mathrm{s}\mathrm{i}\mathrm{n}{\beta}_{j}+\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}{\beta}_{j}}.$
Nonrotating mode shapes are used here because of computational ease. Since they depend only on the fixed constants ${\alpha}_{j}$, ${\beta}_{j}$ and ${\gamma}_{j}$, the modal integrals that result when Galerkin method is applied need to be calculated only once. Substituting Eq. (12) into Eqs. (24, 7) with Galerkin method applied, and assuming:
result in the equations, with 3$N$ subequation structures as follows:
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. [1720] investigates the structural damping of thinwalled composite onecell beams based on modal damping ratio. Modal damping ratio of composite crosssection 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. (24) can be expressed in terms of an equivalent viscous damping matrix with 3×3$N$ as:
${k}_{12}=\left[{\oint}_{}\left(\frac{B}{C}\right)ds/{\oint}_{}\left(\frac{1}{C}\right)ds\right]{A}_{e},$
${k}_{13}={\oint}_{}\left(A\frac{{B}^{2}}{C}\right)yds\left\{\left[{\oint}_{}\left(\frac{B}{C}\right)ds{\oint}_{}\left(\frac{B}{C}\right)yds\right]/{\oint}_{}\left(\frac{1}{C}\right)ds\right\},$
${k}_{14}={\oint}_{}\left(A\frac{{B}^{2}}{C}\right)xds\left\{\left[{\oint}_{}\left(\frac{B}{C}\right)ds{\oint}_{}\left(\frac{B}{C}\right)xds\right]/{\oint}_{}\left(\frac{1}{C}\right)ds\right\},$
${k}_{22}=\left[1/{\oint}_{}\left(\frac{1}{C}\right)ds\right]{A}_{e}^{2},{k}_{23}=\left[{\oint}_{}\left(\frac{B}{C}\right)yds/{\oint}_{}\left(\frac{1}{C}\right)ds\right]{A}_{e},$
${k}_{24}=\left[{\oint}_{}\left(\frac{B}{C}\right)xds/{\oint}_{}\left(\frac{1}{C}\right)ds\right]{A}_{e},$
${k}_{34}={\oint}_{}\left(A\frac{{B}^{2}}{C}\right)xds\left\{\left[{\oint}_{}\left(\frac{B}{C}\right)xds{\oint}_{}\left(\frac{B}{C}\right)yds\right]/{\oint}_{}\left(\frac{1}{C}\right)ds\right\},$
${k}_{44}={\oint}_{}\left(A\frac{{B}^{2}}{C}\right){x}^{2}ds+\left\{{\left[{\oint}_{}\left(\frac{B}{C}\right)xds\right]}^{2}/{\oint}_{}\left(\frac{1}{C}\right)ds\right\},$
${A}_{e}=\frac{1}{2}{\oint}_{}\left(x\frac{dy}{ds}y\frac{dx}{ds}\right)ds.$
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 ${\gamma}_{11}$_{} and ${\gamma}_{12}$_{}can be deduced as [20]:
where ${C}_{C}$ is 4×4 symmetric matrix with its components ${c}_{ij}$_{}embodied. Hence the equivalent damping matrix can be similarly written as [20]:
where the expressions of ${c}_{ij}$ are similar to those of ${k}_{ij}$, but the terms $A\left(s\right)$, $B\left(s\right)$, and $C\left(s\right)$ in ${k}_{ij}$ should be replaced by the equivalent parameter terms ${A}_{d}\left(s\right)$, ${B}_{d}\left(s\right)$, and ${C}_{d}\left(s\right)$ in ${c}_{ij}$, respectively. The expressions of ${A}_{d}\left(s\right)$, ${B}_{d}\left(s\right)$, and ${C}_{d}\left(s\right)$ can be found in [20].
Damping parameter factors (${a}_{s}$) 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°
Integrating structural damping ${C}_{s}$ 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:
with total 3$N$×3$N$ damping matrix ${C}_{SM}$:
4.3. LQC control
In order to realize aeroelastic control for aeroelastic system Eq. (8), defining the state vector $Y={\left[{X}^{T},{\dot{X}}^{T}\right]}^{T}$ and adjoining the identity equation $\dot{X}=\dot{X}$, Eq. (7) can be converted to the state space expression:
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]:
where 𝑄 and 𝑅 are the weighting matrices of system outputs and control inputs, respectively, with control signal set as:
Herein, $P$ is a real symmetric matrix, satisfies the Riccati equation:
The controlled state equation of the closed loop system [23] with state feedback is obtained as:
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 servocontroller that uses the setpoint command $r$ and measurements ${Y}_{0}$_{}to generate the control signal $u$ has integral action to ensure that the output ${Y}_{0}$ tracks the command $r$ with smaller controller amplitude. Here the control error $r{Y}_{0}$ is integrated and the regulator $K$ is designed for the plant augmented with integrated states.
The Kalman filter has the structure of an ordinary stateestimator or observer characterized by:
where ${x}_{0}$ is the state variable in LQG control with its optimal estimate being ${\stackrel{`}{\widehat{x}}}_{0}$; the optimal choice of ${K}_{f}$ minimizes $E\left\{\right[{x}_{0}{\widehat{x}}_{0}{]}^{T}[{x}_{0}{\widehat{x}}_{0}]\}$, which is independent of constant power spectral density matrices $W$ and $V$ in Gaussian stochastic noise processes and given by:
where ${Y}_{x}$ is the unique positivesemidefinite solution of the algebraic Riccati equation:
Fig. 4. LQG controller with integral action and reference input
5. Results and discussions
Aeroelastic stability analysis can be implemented by time domain responses of Eqs. (910). 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. 56 show the uncontrolled blade tip responses of flap bending $y\left(L,t\right)$, lag bending $x\left(L,t\right)$, and twist ${\theta}_{z}\left(L,t\right)$, for ply angles ($\theta $) and pitch angles ($\beta $) 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 $\theta =$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 $\beta =$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 followup study for stability analysis and control of divergent cases will take the fixed angles of $\theta =$30° and $\beta =$45°.
Fig. 5. The responses of bending ($y$), bending ($x$), and twist (${\theta}_{z}$), for ply angles in the range of 0°90°
Fig. 6. The responses of bending ($y$), bending ($x$), and twist (${\theta}_{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, $y\left(L,t\right)$, $x\left(L,t\right)$, and ${\theta}_{z}\left(L,t\right)$, 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
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 $ABK$ in Eqs. (21, 25), respectively) defined in statespace for $U=$ 10 m/s and $U=$ 15 m/s, respectively. Table 23 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 23 that the homogeneous equation system obtained by control of LQR controller is stable, with all the Realparts of eigenvalues being less than zeros. Compared with the controlled system, the Realparts 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 Realparts of uncontrolled system in Table 2 is 0.4278, with the maximum value of Realparts of controlled system being –0.0285. We can get very similar stability conclusions in Table3 with $U=$15 m/s. The maximum value of Realparts of uncontrolled system in Table3 is 2.7322, with the maximum value of Realparts 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 23 are in qualitative agreement with those from time responses in Figs. 78. 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
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 realtime 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 (6$N$) 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 realtime operation.
Fig. 9. a) The constant LQR controller amplitudes of 30 (6$N$) state variables $Y$, b) and LQG controller amplitudes as function of time, in the divergent case of $U=$15 m/s, respectively
a) LQR
b) LQG
In Fig. 10(ac), 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 23) 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 (Eig4Eig1) are very close, which demonstrates the stability rules coincide with those in Fig. 10(ac) obtained via the numerical responses of the three displacements.
Anyway, Fig. 910 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
a) Flap
b) Lag
c) Twist
d) Real parts deviations
6. Conclusions
Bendingbendingtwist 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 BL 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 BL 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 realtime 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
 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. 217226. [Publisher]
 Ng B. F., Palacios R., Kerrigan E. C., et al. Aerodynamic load control in HAWT with combined aeroelastic tailoring and trailingedge flaps. Wind Energy, Vol. 19, Issue 2, 2015, p. 243263. [CrossRef]
 Li N., Mark J. B., Yang H., et al. Numerical investigation of flapwisetorsional vibration model of a smart section blade with microtab. Shock and Vibration, Vol. 2015, 2015, p. 136026. [Publisher]
 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. 587592. [CrossRef]
 Qian W., Huang R., Hu H., et al. Active flutter suppression of a multipleactuatedwing wind tunnel model. Chinese Journal of Aeronautics, Vol. 27, Issue 6, 2014, p. 14511460. [Publisher]
 Gebhardt C. G., Roccia B. A. Nonlinear aeroelasticity: an approach to compute the response of threeblade largescale horizontalaxis wind turbines. Renewable Energy, Vol. 66, 2014, p. 495514. [Publisher]
 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. 487501. [Publisher]
 Khazar H., Sung K. H. Load mitigation of wind turbine blade by aeroelastic tailoring via unbalanced laminates composites. Composite Structures, Vol. 128, 2015, p. 122133. [Publisher]
 Dong O. Y., Kwon O. J. Predicting wind turbine blade loads and aeroelastic response using a coupled CFDCSD method. Control Engineering Practice, Renewable Energy, Vol. 70, 2014, p. 184196. [Publisher]
 Boukhezzar B., Siguerdidjane H. Comparison between linear and nonlinear control strategies for variable speed wind turbines. Control Engineering Practice, Vol. 18, 2010, p. 13571368. [Publisher]
 Hansen M. H., Gaunaa M., Madsen H. A. A BeddoesLeishman Type Dynamic Stall Model in StateSpace and Indicial Formulation. RisøR1354, Risø, Roskilde, Denmark, 2004. [CrossRef]
 Kallesøe B. S. A loworder model for analysing effects of blade fatigue load control. Wind Energy, Vol. 9, 2006, p. 421436. [Publisher]
 Wu X., Jin C., Wenzhong S. Integration study on airfoil profile for wind turbines. China Mechanical Engineering, Vol. 20, Issue 2, 2009, p. 211213. [CrossRef]
 Smith E. C., Chopra I. Formulation and evaluation of an analytical model for composite boxbeams. Journal of the American Helicopter Society, Vol. 36, Issue 3, 1991, p. 2335. [Publisher]
 Liu T., Ren Y., Yang X. Aeroelastic stability analysis based on fitted aerodynamic coefficients. ACTA Energiae Solaris Sinica, Vol. 31, Issue 4, 2010, p. 513516. [CrossRef]
 Hodges D. H., Pierce G. A. Introduction to Structural Dynamics and Aeroelasticity. Cambrige University Press, Cambrige, 2002. [Publisher]
 Ren Y., Du X., Sun S., et al. Structural damping of thinwalled composite onecell beams. Journal of Vibration and Shock, Vol. 36, Issue 3, 2012, p. 141152. [CrossRef]
 Liu T. Theoretical Modeling and Numerical Analysis for Dynamic Stall Aeroelastic Mechanics of LargeScale Wind Turbine Blade. Ph.M. dissertation, College of Mechanical and Electronic Engineering, Qingdao, China, Shandong University of Science and Technology, 2011. [CrossRef]
 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. 14291448. [Publisher]
 Ren Y., Zhang X., Liu Y., et al. Vibration and instability of rotating composite thinwalled shafts with internal damping. Shock and Vibration, Vol. 2014, 2014, p. 123271. [CrossRef]
 Bran A., Armanios A. B. Free vibration analysis of anisotropic thinwalled closedsection beams. AIAA Journal, Vol. 10, 1995, p. 19051910. [CrossRef]
 Ogata K. Modern Control Engineering. Fifth Edition, Publishing House of Electronics Industry, Beijing, China, 2012. [CrossRef]
 Xue D. Computer Aided Control Systems Design Using MATLAB Language. Second Edition, Tsinghua University Publishing Company, Beijing, 2006. [CrossRef]
 Liu T. Classical flutter and active control of wind turbine blade based on piezoelectric actuation. Shock and Vibration, Vol. 2015, 2015, p. 292368. [Publisher]
 Skogestad S., Postlethwaite I. Multivariable Feedback Control: Analysis and Design. Second Edition, Wiley, Chichester, UK, 2005. [CrossRef]
 Liu T. Stall flutter suppression for absolutely divergent motions of wind turbine blade base on Hinfinity mixedsensitivity synthesis method. The Open Mechanical Engineering Journal, Vol. 9, Issue 8, 2015, p. 752760. [Publisher]