Fundamental aeroelastic properties of a bend–twist coupled blade section

The e ﬀ ects of bend – twist coupling on the aeroelastic modal properties and stability limits of a two-dimensional blade section in attached ﬂ ow are investigated. Bend – twist coupling is introduced in the sti ﬀ ness matrix of the structural blade section model. The structural model is coupled with an unsteady aerodynamic model in a linearised state – space formulation. A numerical study is performed using structural and aerodynamic parameters representative for wind turbine blades. It is shown that damping of the edgewise mode is primarily in ﬂ uenced by the work of the lift which is close to antiphase, making the stability of the mode sensitive to changes in the sti ﬀ ness matrix. The aerodynamic forces increase the sti ﬀ ness of the ﬂ apwise mode for ﬂ ap – twist coupling to feather for downwind de ﬂ ections. The sti ﬀ ness reduces and damping increases for ﬂ ap – twist to stall. Edge – twist coupling is prone to an edgetwist ﬂ utter instability at much lower in ﬂ ow speeds than the uncoupled blade section. Flap – twist coupling results in a moderate reduction of the ﬂ utter speed for twist to feather and divergence for twist to stall.


Introduction
The aeroelastic response of wind turbine blades is influenced by the structural coupling between bending and twist of the blade.Bend-twist coupling creates a feedback between the aerodynamic forces, which induce bending moments in the blade, and the blade twist, which is directly related to the angle of attack and thus the aerodynamic forces.The coupling is a result of a curved blade geometry, either from sweeping the planform (flap-twist coupling) (Liebst, 1986;Larwood and Zuteck, 2006) or from prebend or deflections (edge-twist coupling), which introduces an additional torsional component when the blade is subjected to aerodynamic loads.Or from the anisotropic properties of the fibre reinforced polymer (FRP) that is used to build wind turbine blades.When FRP laminates are loaded nonparallel to their principal axes, normal and shear stresses and strains become coupled.Thus, an asymmetric fibre layup in the spar caps and/or skin of the blade results in coupling at cross-section and beam level.While the concept of blade coupling by changing the fibre direction of the FRP originated from helicopter blade application (Mansfield and Sobey, 1979;Hong and Chopra, 1985), Karaolis et al. (1988) are commonly cited to have introduced the method to wind turbine blades.The primary goal of bend-twist coupling in wind turbine blades is to reduce the loads on the turbine (Kooijman, 1996;Lobitz et al., 1996;Lobitz and Veers, 2003).The coupling enables the alleviation of sudden inflow changes, as in gust or turbulent conditions, without the need of external control devices like pitching systems or flaps.Aside from the intended load reduction, the coupling can have a significant effect on the aeroelastic modal properties and stability of the blade (Lobitz and Laino, 1999;Lobitz and Veers, 1998;Hansen, 2011).
The first studies on bend-twist coupling were conducted on stall regulated turbines and hence coupling was introduced to pitch the blade towards stall for flapwise downwind deflections.Lobitz et al. (1996) investigate how much gain in energy production can be achieved by increasing the rotor diameter without overpowering gearbox or generator.Their findings estimate an increase in annual energy production in the order of 10-15%.However, further studies of Lobitz and Laino (1999) show that blades in a twist to stall regime suffer a significant increase in fatigue damage and are prone to stall flutter.Lobitz and Veers (1998) study the aeroelastic behaviour of a coupled blade using a finite element beam model where coupling is introduced by coefficients in the element stiffness matrix.The aerodynamic forces are based on unsteady (Theodorsen) aerodynamics.Lobitz and Veers report an increase in damping of the torsional mode for twist to stall.Increased coupling results in divergence with a significant reduction of the critical inflow speed.Hong and Chopra (1985) investigate the aeroelastic stability of coupled helicopter composite blades in hover using a finite element approach.The coupled beam formulation is obtained by integrating the strain energies over the cross section, thus explicitly considering the fibre layup.Aerodynamic forces are determined by quasi-steady strip theory.The stability of the blades are investigated by means of eigenvalue analysis around a steady state equilibrium.Hong and Chopra report that flap-twist to stall reduces the frequency, and increases the damping ratio, of the flapwise mode.
With the development towards pitch regulated turbines, twist to feather was also investigated.Lobitz and Veers (2003) review blade coupling and compare bend-twist to feather coupling for constant-speed stall-controlled, variable-speed stall-controlled and variable-speed pitch-controlled rotors.All control strategies show significant reductions in the blade root flapwise moment fatigue damage (20-80%) in turbulent inflow.Reduction in ultimate loads are also observed, especially for the variable-speed pitchcontrolled rotor.Lobitz and Veers (1998) report reduced damping of the torsional mode for pitch to feather coupled blades and a moderate reduction of the classical flutter speed.Lobitz (2004) compares the classical flutter limits of an uncoupled and flap-twist to feather coupled 1.5 MW reference blade using a finite element beam model.The stability is analysed by an iterative eigenvalue analysis.Different aerodynamic models are compared in the study and concluded that quasi-steady assumptions lead to drastic underpredictions of the flutter speed.Lobitz further reports a moderate reduction in the classical flutter limit for the coupled blade.Hansen (2011) conducts a study on swept blades.Frequencies, damping and mode shapes are calculated for a beam model of the NREL 5 MW Wind Turbine blade using eigenvalue analysis around a steady state equilibrium.The steady state aerodynamic forces are obtained using the Blade Element Momentum method while a Beddoes-Leishman type dynamic stall model (Hansen et al., 2004) is adopted for the eigenvalue analysis around the steady states.Hansen concludes that the backward sweep, which results in flap-twist to feather coupling, mainly influences the flap mode and has little influence on edgewise vibrations.Aeroelastic frequencies of the flapwise mode are higher while the damping reduces for increased sweep.A moderate reduction of the classical flutter limit is also reported for increased blade sweep.
Few studies on edgewise coupled blades have been published.Hong and Chopra (1985) report reduced frequencies when edgetwist coupling is present.The damping increase for twist to feather for edgewise deflections toward the leading edge and reduces for twist to stall.They conclude that edge-twist coupling has an appreciable influence on stability.
Most studies focus on the analysis of full blades.While this is necessary for the application of bend-twist coupled blades, a full blade analysis makes it difficult to understand the basic mechanisms that alter the aeroelastic response.Rasmussen et al. (1999) therefore investigate the damping of a blade section in attached and separated flow.The edge-and flapwise directions of vibration are prescribed and coupled with in phase and counter-phase pitch motion.Aerodynamic damping is obtained by integrating the aerodynamic work over one cycle of oscillation.For attached flow, they show that edge-twist to feather coupling reduces the damping for edgewise vibration directions between the inflow and the rotor plane.For edge-twist to stall coupling the damping increases.Flap-twist to feather reduces while flap-twist to stall increases the damping in attached flow.
The present study aims to enhance the work of Rasmussen et al. by obtaining the mode shapes from the system matrix of an aeroelastic state-space model instead of prescribing an assumed vibration mode.Consequently, the mode shapes are a combined result of the structural properties and aerodynamic forces.The modal properties for edge-and flapwise coupled sections are investigated.The work of the aerodynamic forces are calculated to examine their influence on the damping ratio.Frequency response and stability of the blade section model are also investigated.The findings are in line with previous studies on full blades.The results further show that the damping of the edgewise mode is dominated by the lift force.As the lift is near antiphase to the displacements, the stability of the mode is sensitive to changes in the stiffness matrix.Edgewise coupling is also prone to an edge-twist flutter mode with a significant reduction of the critical inflow speed.
The next section introduces the structural model of the blade section which has three degrees of freedom that are elastically coupled through the stiffness matrix.The unsteady aerodynamic model, which assumes incompressible and attached flow, is also presented.Dynamic stall effects are not included in the model because the analysis focuses on the aeroelastic properties of a blade operating in normal operation on a pitch regulated turbine.The structural and aerodynamic models are then linearised and combined into an aeroelastic model through a state-space formulation and validated against previous works on classical flutter.The following sections investigate and discuss the aeroelastic frequencies, damping, frequency response and stability limits of a typical wind turbine blade section.The article closes with a summary of the findings.

Aeroelastic model
The aeroelastic model for the modal analysis is derived in this section.The model is split into a structural part and an aerodynamic part.The structural equations of motions are derived by means of Lagrange's equations.The unsteady aerodynamic model is a time domain formulation of Theodorsen's theory.The structural and aerodynamic models are combined in a linearised state-space formulation.

Equations of motion
The blade section was modelled in a Cartesian coordinate system with the origin at the elastic centre (centre of rotation) of the section (see Fig. 1).The section has three degrees of freedom: Edgewise translation parallel to the chord and positive towards the leading edge denoted x, flapwise translation normal to the chord and positive towards the suction side of the aerofoil denoted y, and twist rotation about an axis perpendicular to the section and positive nose-up denoted θ.Consistent with this coordinate system, twist to feather coupling refers to a negative (nose-down) twist, and twist to stall coupling to a positive (nose-up) twist for positive edgewise (towards leading edge) and flapwise (towards suction side) deflections.The chord length of the section is c and the aerodynamic centre (AC) is defined at quarter chord.The twist rotation is around the elastic centre, which is located on the chord at distance e ac aft of the aerodynamic centre.The centre of gravity (CG) is also on the chord e cg aft of the elastic centre.The collocation point (CP) of the aerofoil, at which the quasi-steady angle of attack is evaluated, is at the three-quarter chord (Katz and Plotkin, 2001).
The equations of motion were obtained with Lagrange's equations.Assuming small rotations, the displacement of CG can be written as x x = where k x , k y and k θ are the edge, flap and rotational stiffness.The coupling stiffness are k xθ and k yθ for edge-twist and flap-twist, respectively.For the numerical analysis of this study the stiffness terms were defined as (1) where ω ω ω , , x y θ are the natural frequencies of the uncoupled blade section and γ γ , x y are edge-and flap-wise coupling coefficients as proposed by Lobitz and Veers (1998).In theory the coupling coefficient is limited by the requirement of a non-singular stiffness matrix and can thus range between γ γ −1.0 < , < 1.0 x y .In practice, coupling coefficients in this order are not achievable due to material and manufacturing constraints.To ensure constant structural damping ratios of the individual modes irrespective of the coupling, the damping was introduced by creating a proportional damping matrix C Φ CΦ = ∼ s T s − − 1 from the modal matrix Φ and a diagonal damping matrix C ∼ s containing the modal damping ratios (see (A.2) in Appendix A).The equations of motion in matrix form are

M q
Cq Kq 0 ¨+ ˙+ = s s s s s s (2) where x y θ q = { } s T is the vector of the structural degrees of freedom, and M s , C s and K s are the structural mass, damping and stiffness matrices provided in Appendix A.

Aerodynamic forces
The aerodynamic forces acting on the section are derived under the assumption that the unsteady response can be characterised by a two-dimensional flat-plate aerofoil in attached flow conditions.The geometry of the aerofoil (i.e.camber, thickness) is considered in the steady state response only through the aerodynamic coefficients C L , C D and C M .Before deriving the aerodynamic forces, inflow velocities and angles are defined.The squared relative inflow velocity as shown in Fig. 1 is where the inflow component along the chord line (x-axis) U 0 is constant, the inflow perpendicular to the chord (y-axis) consists of mean inflow V 0 and small variations V 1 , and x ˙and y ˙are the velocities of the blade section as defined above.The geometric angle of attack between chord and free-stream flow defines the direction of the aerodynamic forces: The quasi-steady angle of attack 0 is calculated at the collocation point and used to determine the magnitude of the forces.

Lift
The unsteady aerodynamic lift acting on the aerofoil can be split into three parts (Fung, 1969).Lift due to the circulatory airflow around the section L circ , a force due to apparent mass L iner and a force of centrifugal nature L cent .For attached flow, the circulatory contribution is obtained from where ρ is the density of the air and C α ( ) L E is the lift coefficient at the effective angle of attack α E .To account for the shed vorticity of the circulatory lift, the Theodorsen part of the Beddoes-Leishman type dynamic stall model proposed by Hansen et al. (2004)

is implemented. The effective angle of attack thus becomes α
2 where z i for i ∈ {1, 2} are state variables of the unsteady aerodynamic model obtained from the first order ordinary differential equations The variables A i and b i depend on the aerofoil, for a flat plate they were derived by Jones (1940) as: The apparent mass term L iner is caused by the inertial force of the air surrounding the section.The inertial force acts at mid-chord and equals the mass of air in a cylinder with the diameter of chord c times the vertical mid-chord acceleration: The force of centrifugal nature L cent is caused by the directional change of the apparent mass.This force is acting at the collocation point and obtained as With the three contribution the total lift becomes

Drag
The total drag is modelled as where C D is the function of the quasi-steady drag coefficient.The first term is the viscous drag.Under attached flow conditions, the viscous drag is dominated by friction drag and it changes little with the angle of attack.The second term is the induced drag which, in two dimensional flow, is caused by the effective angle of attack α E lagging behind the geometric angle of attack α.Thus, the unsteady lift, which is perpendicular to α E , has a component in the drag direction defined by the geometric angle of attack.

Moment
For symmetric aerofoils the aerodynamic centre is positioned at quarter chord.If the aerofoil is cambered, this centre moves aft which is incorporated by a moment coefficient C M .Apart from this offset, the lift forces L iner , L cent and the apparent moment of inertia ρπc 128 4 (see e.g.Fung, 1969) contribute to the moment about quarter chord which is

Linearised aerodynamic forces
As this study aims to investigate the linear dynamic response around a steady state equilibrium, the flow velocities, angles and aerodynamic forces are linearised.Ignoring higher order terms for small variations and velocities V x y , ˙, ˙⪡1 1 , the inflow (3) can be split into a mean W 0 and variable part W 1 as And similarly for the geometric angle of attack (4) using a first degree Taylor series With the lift-curve slope C α ( ) L α , 0 at the mean angle of attack α 0 , the steady state value L 0 and first order variation L 1 of the lift force (10) are A linearisation of the unsteady aerodynamic model ( 6) has been derived by Hansen et al. (2004) as where is the first order variation of the quasi-steady angle of attack and z i1 are the first order variations of the aerodynamic states.With a constant drag coefficient C D , the first order variation of the drag force (11) becomes where in the last term it was accounted for that the variation of the geometric angle of attack α 1 around the steady state angle of attack α 0 , relative to which the linearised forces are defined, induces a lift component in the steady state drag direction of magnitude L α − 0 1 .Assuming a constant moment coefficient C M , the first order variation of the moment ( 13) is

State-space representation
The linearised aerodynamic model ( 18) can be written in matrix form as  19) and (20) the first order variations of the aerodynamic forces can be expressed as where T is a matrix that transforms the aerodynamic forces, which are relative to the steady state inflow, into the chord coordinate system of the blade section.The structural state, aerodynamic state and wind input vectors are q s , q a and V V u = { ˙}T 1 1 .The aerodynamic mass, damping and stiffness matrices are M a , C a and K a , and A f and W f represent the contribution of unsteady flow model and wind variation to the aerodynamic forces.The matrices T M C K A , , , , a a a f and W f are also provided in Appendix A. The equations of motion (2) with the first order variation of the aerodynamic forces (22) on the right hand side and the aerodynamic model ( 21) can be rearranged to The coupled aeroelastic (23) and unsteady flow (24) model can be rewritten in first-order form as where A is the state matrix of the aeroelastic model and B is the input matrix of the wind speed variations.A Python implementation of the state-space model is available on GitHub. 1

Aerodynamic work and damping ratio
To investigate the individual contributions of lift, drag and moment to the damping ratios of the aeroelastic modes, the work of the aerodynamic forces over one cycle of harmonic oscillation at modal frequency ω d is obtained from where q T q ˙= ∼s s −1 are the velocities in the direction of the aerodynamic forces.The aerodynamic forces and the velocities can be expressed as where ∼ are the amplitudes of the aerodynamic forces and displacements (in force direction), and φ φ φ , , are the phases of the displacements and φ φ φ Δ , Δ , Δ D L M are the phase differences between the displacements and the forces.Substituting ( 27) into (26) and integrating yields Eq. ( 28) shows the contribution of each force component to the aerodynamic work.The magnitude of the force contribution depends on the amplitude of the force and the amplitude of the mode shape in force direction.The sign of the force contribution depends on the phase between the displacement and the force.If the force has a positive phase angle relative to the displacement, the contribution is negative and energy is accumulated in the system.For negative phase angles energy is extracted from the system.The aerodynamic work, together with the structural damping, contributes to the energy dissipated over one cycle of harmonic oscillation E dis which can be related to the damping ratio by where E tot is the total (kinetic + potential) energy of the system.The approximation has been derived in Appendix B.

Validation
To validate the present model, the following non-dimensional parameters have been used to compare classical flutter speeds with results provided by Zeiler (2000): Fig. 2 shows that the flutter speeds obtained with the present model are in good agreement with the results by Zeiler.The small discrepancies are probably related to the digitising process to obtain the data from the original graph.

Numerical analysis and results
The present model is analysed using the following parameters:  Turbine (Bak et al., 2013).The section radius was chosen on the outer third of the blade where most of the aerodynamic work is performed.However, it is not necessarily the best representation of the actual mode shape of the blade.The structural frequencies ω x , ω y , ω θ and damping ratios ζ x ∼, ζ y ∼, ζ θ ∼ correspond to the first flapwise bending, first edgewise bending and first torsional modes of the DTU 10 MW RWT blade at standstill.The inflow speed is set to 45 m/s with a steady state angle of attack of α = 7°0 which corresponds approximately to the inflow at 75% blade radius at a wind speed of 8 m/s with an induction factor of a=0.3 and a tip speed ratio of 7.5.The 75% blade radius is chosen because it is considered to be the point with the highest aerodynamic forces.Inboard the inflow velocity reduces with the radial position, and outboard the shorter chord and the tip losses reduce the lift.The lift, drag and moment coefficient were assumed as C α = 7.15

Frequencies, damping and aerodynamic work
To investigate the aeroelastic modal and stability properties of the blade section, eigenvalue analysis of the state matrix A of Eq. ( 25) was used to calculate modal frequencies and damping ratios for different coupling coefficients.Fig. 3 shows contour plots of aeroelastic frequencies (left plots) and damping ratios (right plots) of edgewise (top), flapwise (middle) and torsional (bottom) mode for varying edge-twist (ordinate) and flap-twist (abscissa) coupling coefficients.Looking at the top left contour plot, the frequency of the edgewise mode reduces for edge-twist coupling to both feather and stall and changes little for flap-twist coupling.Edgewise damping in the top right plot becomes negative for edge-twist to stall coupling if the flap-twist coupling coefficients is γ < 0.2 Looking at the bottom plots the frequency and damping of the torsional mode are proportional over the full range of flap-twist coupling coefficients.The frequency increases while the damping decreases for flap-twist to feather coupling.The effects of edgetwist coupling on frequency and damping of the torsional mode are small.Fig. 4 shows 2D plots of selected sections of the contour plots in Fig. 3, as indicated in those plots.The sections show both frequencies and damping ratios for both structural and aeroelastic modes over the coupling parameters. .Comparing the frequencies of the structural and aeroelastic modes, it is seen that coupling of the aerodynamic forces with the motion increases the frequency of the mode for flap-twist to feather coupling and reduce it for flaptwist to stall.The damping of the aeroelastic flapwise mode is significantly higher than for the structural flapwise mode due to the well known high aerodynamic damping of flapwise vibrations in attached flow (Rasmussen et al., 1999).The aeroelastic damping reduces slightly for flap-twist to feather and increases for flap-twist to stall coupling.Fig. 4d shows a cut through the contour plots of the torsional mode along the flap-twist coupling direction for an edge-twist coupling coefficient of γ = 0.0 x .The frequency slope is associated with the structural mode while the aerodynamic forces reduce the frequency by a constant 0.25 Hz.The aeroelastic damping is above the structural and proportional to the coupling.
To understand the effects of the couplings on the aeroelastic damping of the three modes, the work over one period of oscillation is computed for each mode and for each force and then normalised with the total energy π E 4 tot of the mode so that the total normalised work approximates the damping ratio as shown in Eq. ( 29).Positive work means that the forces extract energy from the system (stabilising) while negative work means that the energy in the system is increased (destabilising).
Fig. 5 shows the normalised aerodynamic work (left), structural damping work (middle), and total (aerodynamic + structural) work (right) for the edgewise, flapwise, and torsional mode for varying edge-twist (ordinate) and flap-twist (abscissa) coupling coefficients.Except for the edgewise mode, the aerodynamic work is always positive.The structural damping work is always positive.The total normalised work is a good approximation of the damping ratio shown in Fig. 3 if the damping ratio is small.Fig. 6 shows the normalised aerodynamic work of drag (left), lift (middle), and moment (right) for the edgewise (top), flapwise (middle), and torsional (bottom) mode for varying edge-twist (ordinate) and flap-twist (abscissa) coupling coefficients.The aerodynamic work of the edgewise mode is dominated by the lift force and it is negative or close to zero except for strong edge-twist coupling if flap-twist to stall coupling is strong.Eq. ( 28) shows that the phase determines the sign of the aerodynamic work and for the edgewise mode the phase between lift force and motion is close to −180°.The lift work is therefore sensitive to changes in the structural stiffness which can cause the phase angle to become positive, resulting in negative damping and instability.The drag work of the edgewise mode increases with edge-twist to feather coupling and reduces for edge-twist to stall.The work of the moment increases with edge-twist coupling.The aerodynamic work of the flapwise mode is also dominated by the lift work.However, the phase angle is around −90°(displacement leads force) and changes in the mode shape have little influence on the damping ratio.The drag work of the flapwise mode is mostly negative.The moment work is negative for flap-twist to feather coupling and positive for flap-twist to stall coupling.The torsional mode is dominated by the positive work of the moment which can largely be attributed to the pitch rate term associated with θ ˙in Eq. ( 20).The lift work is negative for all coupling coefficients, its magnitude is about 30% of the moment work.Drag work is small and positive except for strong edge-twist to feather coupling.

Frequency response
Fig. 7 shows the frequency responses of the blade section due to mean wind speed variations V 1 for the uncoupled and flap-twist coupled sections.The corresponding H ∞ -norms, which are the peak gain of the frequency response and related to ultimate loads, and H 2 -norms, which are an average gain over all frequencies and related to fatigue loads, are given in Table 1.Due to the high aerodynamic damping there is no distinct peak in the flapwise response and the H ∞ -norm is less relevant.Flap-twist coupling has little influence on the edgewise magnitude.The flapwise magnitudes is reduced for twist to feather and increased for twist to stall compared to the uncoupled section.The torsional magnitude increases significantly around the flapwise frequency (0.61 Hz) for flap-twist to feather and stall coupling.The H 2 -norm of the edgewise response increases (24%) for twist to feather and decreases (10%) for twist to stall.For the flapwise response the H 2 -norm reduces (10%) for twist to feather and increases (31%) for twist to stall.The H 2 -norm increases for the torsional response, both for twist to feather and stall.The change of the frequency response is different for flap-twist coupling to feather and stall and therefore not linear.

Stability
The stability analysis was conducted by stepwise increasing the inflow speed at zero steady state angle of attack and checking for eigenvalues with positive real parts in each step.Fig. 8 shows the critical inflow speeds for both edge-and flap-twist coupling over the coupling range of γ −0.5 < < 0.5 where the other coupling is zero.In the middle range of edge-twist coupling with coefficients of around γ −0.2 < < 0.2 x , the critical inflow speed of 185 m/s remains constant.The instability is classical flutter (predominantly flapwise and torsional mode shape with torsion leading flap towards 90°).For coupling coefficients γ | | > 0.2 x , the critical inflow speed for the edge-twist coupled section undergoes a steep drop.This drop is caused by an edge-twist flutter mode resulting from interaction of the flapwise mode and the torsional component of the coupled edge-twist mode.
The stable domain of the flap-twist coupled section in the right plot of Fig. 8 is limited by classical flutter with slightly decreasing critical speed towards pitch to feather and divergence for flap-twist to stall which becomes critical at a coupling coefficient of about γ = 0.03 y .Fig. 9 shows a contour plot of the critical inflow speed for the entire domain of edge-and flap-twist coupling.The reduced inflow speed in the right half of the plot is caused by divergence as also shown in Fig. 8.The reduced inflow speeds in the top and bottom thirds of the left half of the plot are caused by edge-twist flutter.The remaining plateau represents the classical flutter limit with a slight decrease in critical inflow speed towards flap-twist to feather.
From the computation of the aerodynamic work over one cycle of oscillation presented in Fig. 6, it was shown that the negative damping of the edgewise modes is caused by the lift force, which can be explained by looking at the mode shape of the edge-twist flutter mode.Figs. 10 and 11 show the motion of and the forces on the section over one period for the edge-twist flutter modes.Instead of torsion, as in classical flutter, the coupled edgewise motion and torsional rotation is ahead of the flapwise motion.This phase difference results in a small phase shift of the lift, which is now leading the flapwise displacements, resulting in the aerofoil extracting energy from the air stream.For aeroelastic flutter to occur the torsional component of a coupled edge-twist mode hast to be sufficiently large.For the present blade section, the structural mode shape of the coupled edge-twist mode with a coupling coefficient of γ = 0.2 x has a torsional component of about 2°at 1 m edgewise deflection.The inflow speed at which edge-twist flutter becomes critical also depends on the frequency difference between edge-and flapwise mode.A larger difference delays the onset of the edge-twist flutter mode.The main difference of the critical modes for edge-twist to feather and stall is the direction of rotation.For aeroelastic flutter to occur, the lift, and thus the angle of attack, has to lead the flapwise motion which results in an anti-clockwise circulating mode for edge-twist to feather and a clockwise circulating mode for edge-twist to stall.The high coupling coefficients in Figs. 10 and 11 where chosen to illustrate the modes.For more moderate coupling coefficients the edgewise and torsional amplitude reduce.
Fig. 12 shows frequencies and damping ratios of the three aeroelastic modes over inflow speed for zero and different edge-and flap-twist couplings.Classical flutter occurs where the torsional and flapwise modal frequencies approach each other and the mode shapes interact, over the aerodynamic and structural couplings, leading to a negatively damped combined mode.The figure shows Fig. 6.Normalised work over one cycle of harmonic oscillation at modal frequency of aerodynamic drag (left), lift (middle), and moment (right) for the edgewise (top), flapwise (middle), and torsional (bottom) mode for varying edge-twist (ordinate) and flap-twist (abscissa) coupling coefficients.Negative coupling coefficients denote twisting to feather for edgewise deflection towards the leading edge and flapwise deflection towards the suction side respectively.The inflow velocity is 45 m/s with a 7°angle of attack.that a classical flutter mode forms at an inflow speed of approximately 160 m/s with a sudden decrease in the previously increasing torsional damping for all chosen coupling coefficients.The damping continues to decrease until the mode becomes negatively damped at around 190 m/s inflow.For the edge-twist coupled cases (2nd and 3rd row), the flapwise frequency slightly increases with wind speed and when it gets in the vicinity of the edgewise frequency, at an inflow speed of about 90 m/s, the torsional component of the coupled edge-twist mode leads to flutter with slightly negative damping of the edge-twist mode.The flap-twist to feather case (4th row) shows that the reduced classical flutter speed can be attributed to an increased slope of the torsional damping after the flapwise and torsional modes have approached each other.This increased slope leads to a negative damping ratio at lower inflow speeds.The point at which flap and torsional mode join, however, changes little with coupling and remains around an inflow speed of 160 m/s.For flap-twist to stall (5th row) the damping slope flattens and the critical inflow speed is increasing.Classical flutter for this case, however, is well beyond the inflow speed at which the flapwise mode becomes divergent.

Discussion
In this section, the aeroelastic properties of the bend-twist coupled blade section are discussed and compared to previous studies.For edge-twist coupling, the section model shows that the edgewise frequency is reduced for both, twist to feather and twist to stall.The frequency change is caused by the structural coupling which reduces the stiffness of the edgewise mode (cf.Fig. 4a and  b).Similar observations are made by Hong and Chopra (1985) in their numerical study on material coupled rotor blades who attribute the decreased edgewise frequency to the reduced bending stiffness for non-zero ply angles.For coefficients until about γ | | = 0.3 x edgewise damping of the blade section increases for edge-twist to feather and reduces for edge-twist to stall coupling (cf.Fig. 4a and b).The same observation is made by Hong and Chopra who subsequently conclude that edge-twist coupling has an appreciable influence on stability.Rasmussen et al. (1999) make a different observation with their blade section model where edgetwist to feather coupling reduces the damping (and edge-twist to stall increases damping).The difference is probably related to the prescribed mode shapes in the study of Rasmussen et al. for which the edgewise, flapwise and torsional components are either in-Fig.7. Aeroelastic frequency response of edgewise (left), flapwise (middle) and torsional (right) amplitude to wind speed variations V 1 for the uncoupled and flaptwist coupled section.The mean inflow velocity is 45 m/s with a 7°angle of attack.phase or counter-phase to each other (i.e.intermediate phase shifts have not been investigated).It is interesting to observe that edge-twist coupling has little influence on the flapwise and torsional modes of the section and it is not transferred further by the inertia coupling terms in the structural mass matrix, or the coupling resulting from aerodynamic forces (cf.Fig. 3).Flap-twist to feather increases the frequency and reduces damping due to an increase in aerodynamic stiffness and a slight reduction in aerodynamic work (cf.Fig. 4c).Rasmussen et al. (1999) also observe a reduction in damping for flap-twist to feather.For backward swept blades (i.e.flap-twist to feather) Hansen (2011) observes increased frequencies and reduced damping for increased sweep.
For flap-twist to stall, the frequency of the section reduces and damping increases (cf.Fig. 4c).Hong and Chopra (1985) also show that the frequency increases for flap-twist to feather and decreases for flap-twist to stall.The torsional frequency of the blade section behaves inverse proportional to flap-twist coupling while the damping is proportional to the coupling (cf.Fig. 4d).The latter is also observed by Lobitz and Veers (1998).
The edgewise, flapwise and torsional frequency response to wind speed fluctuation is an indicator of the blade section response to turbulence or gusts and can therefore serve as a qualitative estimate for the structural loads.The reduction of the H ∞ and H 2 -norms of the flapwise frequency response by 8% and 10% for flap-twist to feather (γ = −0.3y ) can be interpreted as a reduction in ultimate and fatigue loads (cf.Table 1).A reduction of ultimate and fatigue loads is also observed by Lobitz and Veers (2003) who investigate a flap-twist to feather coupled rotor with a coupling coefficient of γ = −0.6 y .The reduction in flutter speed for flap-twist to feather (cf.Fig. 8) has previously been observed by Lobitz and Veers (1998) and Lobitz (2004).The latter reports a moderate reduction in flutter speed of about 15% for the 1.5 MW baseline WindPACT blade with a coupling coefficient of γ = −0.4y .The blade section in this study becomes prone to divergence for flap-twist to stall at relatively low coupling coefficients (cf.Fig. 8).Divergence for flap to stall coupling has previously been reported by Lobitz and Veers (1998).
While the findings obtained with the blade section model are generally in good agreement with previous studies on full blades, the results obtained with such a section model cannot necessarily be extended to full blades.It should further be noted that the flow velocities at instability are beyond the validity of the incompressible flow assumptions that underlie the aerodynamic model.While the assumption of incompressibility is valid for wind turbine blade analysis in the operational range with tip speeds typically below 100 m/s, the inflow velocities in, for example, a runaway situation can be close to transonic flow.For those high inflow velocity scenarios, results obtained with an incompressible aerodynamic model should be treated with caution.

Conclusion
The aeroelastic response of a bend-twist coupled wind turbine blade section is investigated.The blade section has three degrees of freedom: edgewise, and flapwise translation and twist rotation.The structural stiffness of the section is assumed linear and bendtwist coupling is introduced by means of a coupling coefficient in the stiffness matrix.An unsteady (Theodorsen) aerodynamic model is implemented to account for the effects of shed vorticity.Induced unsteady drag is accounted for by the phase lag of the lift force vector.The aerodynamic forces are linearised around a steady state equilibrium and the modal and stability characteristics are analysed by means of eigenvalue analysis of the coupled structural and aerodynamic models.The numerical analysis is carried out with blade section properties at approximately 75% span of the DTU 10 MW Reference Wind Turbine.
It is shown that under normal operation, structural coupling mainly affects frequency and damping of the two components that are being coupled (i.e.edge and twist mode for edge-twist coupling and flap and twist mode for flap-twist coupling), and is not transferred further by the inertia coupling terms in the structural mass matrix, or the coupling resulting from aerodynamic forces.Edge-twist to feather coupling of the blade section can increase edgewise damping.Edgewise damping is reduced for twist to stall.The damping ratio of the edgewise mode is primarily influenced by the work of the lift which is close to antiphase.The stability of the mode is therefore sensitive to changes in the stiffness matrix which can cause the lift force to lead the motion resulting in an increase in total energy of the system.
Flap-twist to feather increases the frequency and reduces damping due to an increase in aerodynamic stiffness and a slight reduction in aerodynamic work.Flap-twist to stall on the other hand is shown to reduces frequency and increase damping mainly due to reduced aerodynamic stiffness and an increase in the aerodynamic work.The flapwise damping is dominated by the work of the lift but its stability is not sensitive to the coupling coefficients as the phase difference is around −90°(displacement leading force).The torsional mode is mainly influenced by flap-twist coupling.Frequency and damping are proportional to the coupling with increasing frequency and slightly reduced damping for flap-twist to feather.The torsional damping relies on the pitch rate damping term to remain positively damped.
The flapwise frequency response of flap-twist coupled blades to wind speed variations results in reduced H ∞ (peak gain) and H 2 (an average gain) norms for twist to feather while both norms increase for twist to stall.
Edge-twist coupling has a significant effect on the stability for coupling coefficients γ | | > 0.2 x .The critical inflow speed is reduced significantly when the torsional component of the edge-twist mode becomes large enough to enable the formation of an edge-twist flutter mode.Flap-twist to feather results in a moderate reduction of the classical flutter limit.Flap-twist to stall leads to divergence with a steep decrease in critical inflow speed.
The findings show that aeroelastic tailoring of wind turbine blades is limited by stability considerations.The investigated blade section showed unstable behaviour, for flap-twist to stall and edge-twist to stall and feather, even for moderate coupling terms.Safety margins will further reduce the coupling range available for aeroelastic tailoring, limiting its potential applications.While the simple blade section model gives a physical understanding of the modal properties and stability limits of bend-twist coupled blades, the findings cannot necessarily be extended to full blades.where Φ is the modal matrix constructed from the mass normalised edgewise, flapwise and torsional mode shapes of the undamped system and ζ q ∼ and ω q ∼ for q x y θ ∈ { , , } ∼ ∼ ∼ ∼ are the damping rations and natural frequencies of the coupled blade section.The matrices of the aerodynamic model are (A.15)

Appendix B. Damping ratio approximation
The approximated damping ratio in Eq. ( 29) is derived as follows.The total energies U (kinetic + potential) of a multiple degree of freedom system can be written as where t x( ) is the vector of generalized coordinates, t x ˙( ) its time derivative, and M and K are the mass and stiffness matrices.Assuming oscillations in a single eigenmode the generalized coordinates can be expressed as Let U 0 be the total energies at t=0.The total energies after one cycle of oscillation U T at t = the rotational inertia as J mr = 2 , with m being the mass of the section and r the radius of gyration about CG, the kinetic energy of the system is T m

Fig. 2 .
Fig. 2. Comparison of classical flutter speeds for the present model (lines) with digitised data from a plot by Zeiler (2000) (dots).

Fig. 3 .
Fig. 3. Contour plots of aeroelastic frequencies (left plots) and damping ratios (right plots) of edgewise (top), flapwise (middle) and torsional (bottom) mode for varying edge-twist (ordinate) and flap-twist (abscissa) coupling coefficients.Negative coupling coefficients denote twisting to feather for edgewise deflection towards the leading edge and flapwise deflection towards the suction side.Positive coupling coefficients twist towards stall.The inflow velocity is 45 m/s with a 7°angle of attack.
-twist to stall coupling the damping increases.Damping also increases for edge-twist to feather coupling if γ > − 0.2 y .For strong edge-twist to feather coupling and γ < − 0.2 y damping becomes negative.Looking at the middle left plot the frequency of the flapwise mode increases slightly for flap-twist to feather and reduces for flaptwist to stall.Flapwise damping in the middle right plot decreases slightly for flap-twist to feather and increases for flap-twist to stall.Edge-twist coupling has little effect on both frequency and damping of the flapwise mode.
Figs. 4a and b show the 2D cuts through the contour plots of the edgewise mode along the edge-twist coupling direction for flap-twist coupling coefficients of γ = −0.2y and γ = 0.2 y .It can be seen in both plots that the frequency change of the aeroelastic edgewise mode follows the frequency change of the structural mode which shows that these effects of the edge-twist coupling are not governed by the coupling with the aerodynamic forces.Aeroelastic damping of the edgewise mode for flap-twist coupling of γ = −0.2y increases in Fig. 4a for edgetwist to feather coupling until a coefficient of γ = −0.3x , and decreases for edge-twist to stall.In Fig. 4b the damping of the edgewise mode for flap-twist coupling of γ = 0.2 y increases for edge-twist to feather and reduces for edge-twist to stall until a coefficient of γ = 0.3 x .Fig. 4c shows a cut through the contour plots of the flapwise mode along the flap-twist coupling direction for an edge-twist coupling coefficient of γ = 0.0 x

Fig. 4 .
Fig. 4. Selected sections of the contour plots in Fig. 3.The sections show aeroelastic frequencies and damping ratios and structural frequencies and damping ratios of different mode shapes over coupling parameters γ = 0.0 x

Fig. 5 .
Fig.5.Normalised work over one cycle of harmonic oscillation at modal frequency of aerodynamic forces (left), structural damping (middle), and total (aerodynamic + structure) (right) for the edgewise (top), flapwise (middle), and torsional (bottom) mode for varying edge-twist (ordinate) and flap-twist (abscissa) coupling coefficients.Negative coupling coefficients denote twisting to feather for edgewise deflection towards the leading edge and flapwise deflection towards the suction side respectively.The inflow velocity is 45 m/s with a 7°angle of attack.

Fig. 8 .Fig. 9 .
Fig. 8. Critical inflow speed (m/s) for edge-(left) and flap-twist (right) coupling.Edge-twist coupling leads to a sudden drop of the critical inflow speed due to edge-twist flutter at a coupling coefficient of about γ | | = 0.2 x
transform the aerodynamic forces, which are relative to the steady state flow direction, is is the complex eigenvector and λ β iω = − + d the complex eigenvalue with attenuation rate β and modal frequency ω d .With Eq. (B.2) the total energies of a single mode can be expressed as systems the modal frequency ω d is close to the natural frequency ω n , hence ζ = ≈ degree Taylor series of the logarithm x ln around x=1 yields the approximation x x ln(1 + ) ≈ .With both approximations and by introducing the dissipated energy over one cycle of oscillation U (B.6) one obtains Eq. (29).

Table 1 H
∞ and H 2 -norm of frequency response to wind speed variations V 1 .