ANALYSIS OF WIND TURBINE BLADE VIBRATION AND DRIVETRAIN LOADS
By
Venkatanarayanan Ramakrishnan

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Mechanical Engineering – Doctor of Philosophy
2017

ABSTRACT
ANALYSIS OF WIND TURBINE BLADE VIBRATION AND DRIVETRAIN LOADS
By
Venkatanarayanan Ramakrishnan
The reliability of wind turbines is a major issue for the industry. Drivetrain and blade failures are
common, costly and not fully understood. Designers must thus examine and understand the key
parameters that influence reliability. As wind turbines increase in size, the blades are designed to
be more lightweight and flexible, increasing the potential for large-displacement oscillations during
operation. This necessitates the incorporation of nonlinearity in the formulation of the blade model
to better understand the dynamics and stability characteristics. Also, oscillations in the blade impart
dynamic loading onto the gearbox. Understanding these dynamic loads is essential for the design
of reliable gears and bearings, and hence economically viable wind turbines. Traditional studies
of wind turbines have focused on the aerodynamic performance of the blades, the reliability of
gearbox and its components, grid connectivity and improvements in power distribution. The aspect
of blade vibration from a dynamics point of view has garnered interest but not been fully developed
and understood.
In this work, the partial and ordinary differential equations that govern the in-plane and outof-plane motion of a wind turbine blade subject to gravitational and aerodynamic loading are
developed using Hamilton’s principle and Lagrange formulations respectively. These differential
equations include nonlinear terms due to nonlinear curvature and nonlinear foreshortening, as well
as parametric and direct excitation at the frequency of rotation. The equations are reduced using
an assumed uniform cantilevered beam mode to produce single second-order ordinary differential
equations (ODE) to approximate the blade model for the case of constant rotation rate. Embedded
in the ODE’s are terms of a forced Mathieu equation with cubic nonlinearity. Different variations
of the forced Mathieu equation are analyzed for resonances by using the method of multiple
scales. The forced Mathieu equation has instabilities and resonances at multiple superharmonic

and subharmonic frequencies. Second-order expansions are used to unfold the expressions that
govern the amplitude of response at these critical resonances. The equations of motion (EOM’s) also
have regions of instability and we employ perturbation analysis to identify the stability transition
curves of the system. These calculations compare well with numerical simulations for simple
systems under study.
The formulation is then extended to wind turbine blades. Aerodynamic forces on the wind
turbine blades are calculated using the Blade Element Momentum (BEM) theory and its extensions.
From parametric studies of the blade EOM we can understand the parameter values at which
superharmonic resonances become dominant. It is shown that as wind turbine blades become
larger they are prone to these resonances, whose existence may not be within the scope of current
design strategies. The amplitude of response at all resonances tend to become amplified for much
larger blades. Both in-plane and out-of-plane responses will increase the loading at the rotor hub
and consequently, increase the loads and moments on the wind turbine drivetrain.
To capture the effect of increased loading on the wind turbine drivetrain, we use a commercial
software to model a 750 kW gearbox used as a part of the Gearbox Reliability Collaborative
(GRC) headed by National Renewable Energy Labs (NREL). For this, we partnered with Romax
Technology Ltd. to analyze the sensitivities of the load on the elements of the gearbox to variations
in the input loads. Using the Romax gearbox model, we suggest methods to optimize the gear
geometry to improve reliability of the drivetrain by minimizing influence of manufacturing and
assembly tolerances and misalignments. We also designed novel approaches to predict gearbox
vibration using the models and suggested changes that are required to improve the overall design
of the gearbox (these have been implemented on newer generations of the NREL GRC gearbox).
The fundamental work formulated in this thesis can be extended to more complex models to
understand other system level dynamics of interest (multi-mode interaction, multi-blade resonance,
etc.). More detailed formulation of aerodynamic loads (for example by using ONERA semiempirical approach) would also improve model fidelity for predicting the influence of aerodynamic
loads on blade vibration.

Copyright by
VENKATANARAYANAN RAMAKRISHNAN
2017

To my parents (Ramakrishnan and Chitra) and loving sister (Lalitha).

v

ACKNOWLEDGEMENTS

I take this opportunity to acknowledge the many individuals who have made an indelible impression
on my personal and professional life during my years in graduate school and have in their own ways
made this work possible.
First and foremost, my regards goes to my mentor and advisor, Professor Brian Feeny. Brian is
one of the best human beings I have known and it has been an absolute pleasure to work under his
guidance. His patience with my pace of work, insight into understanding dynamical systems, and
command over technical writing have been instrumental in completing this thesis. I have always
looked up to Brian and sought his advice on a host of academic and personal situations, and he has
guided me along every step of the way.
I would like to extend my gratitude to Professor Steve Shaw and Professor Keith Promislow. My
interest in dynamic systems, vibrations and applied mathematics stems from taking graduate level
courses taught by them. Their enthusiasm during those lectures and my subsequent interactions
with them at the lab, and other social gatherings helped me formulate a lot of ideas and develop a
good understanding of mathematical modeling.
I would like to thank Professor Farhad Jaberi and Mr. Brian Wilson from Romax Technology,
for serving on my dissertation committee. Brian was actively involved in defining the scope of my
thesis work and provided valuable insights and discussions.
I have spent three semesters as a graduate research intern with Romax Technology which
has given me a lot of exposure to real world engineering problems and has been instrumental in
developing my understanding of drivetrains and wind turbine systems. I have worked under Dr.
Ashley Crowther, Dr. Chris Halse, Dr. John Coultate and Mr. Zach Wright during various projects
and these gentleman have set very fine examples of an engineering professional I aspire to be.
Working with Professor Ranjan Mukherjee for my M.S. inculcated the necessary discipline
and work ethic to pursue another advanced academic degree. I wish to thank him as well for the
necessary guidance and help, to explore my passion.
vi

I would like to thank Mr. Tom Corden with whom I have worked on developing VIBEX, a
vibration damping gel, which has been very rewarding experience and often a good distraction.
Thanks are also due to Craig Gunn for the numerous opportunities to work with outreach
activities and to Aida Montalvo, Mary Pease and the rest of the ME staff for supporting me in
various capacities at numerous times.
I have had the pleasure to work with some of the most intelligent and fun graduate students and
the many discussions in the lab. Nick, Scott, Rickey, Brendan, Ryan, Thomas, Carl, Xing, Abhisek,
Pavel, Smruti - thank you for the support and company.
My gratitude goes to my friends, who are in every sense, my extended family at MSU. Abinand,
Mohit, Sumit, Avinash, Aritro, Satish, Arun, Chris, Marco, Sridharan, Nanda, Ren, Kaustav,
Shreelina, Disha, Anchita, Anuradha, Shaheen, Sreekala, Nikita, Shalini, Priyamvada, to name a
few with whom I have spent many eventful days and evenings. I have also met many interesting and
diverse people on campus who have had much positive influence on me and shaped me to become
a better person.
I would like to thank my mentor at Chrysler, Dr. William Resh, for the necessary support and
encouragement to finish this thesis.
Finally, I would like to recognize the love, support and encouragement I receive from the love
of my life, Sneha, whom I met at MSU during the course of my graduate study. She has a been a
constant source of support for me.
The work associated with this dissertation was supported by funding made possible by National
Science Foundation Grant (CBET-0933292) and the MSU College of Engineering. Software
support was provided by Romax Technology. Additional support came from the Graduate Research
Enhancement Fellowship.
I came to the US as a young 21 year old to pursue my dreams of a higher education in
Engineering. East Lansing, MSU and its people have a very positive and significant impact and
have helped me grow professionally and personally. I will always remember this place and the
lovely people I have known through it very dearly.

vii

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

x

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xi

CHAPTER 1 INTRODUCTION . . . . . . . . . .
1.1 Objective . . . . . . . . . . . . . . . . . .
1.2 Background and Motivation . . . . . . . . .
1.3 Wind Turbine Blade Modeling . . . . . . .
1.3.1 Aerodynamic Modeling . . . . . . .
1.3.2 Aeroelastic Modeling . . . . . . . .
1.3.3 Overview of Design Codes . . . . .
1.4 Beam Vibrations and the Mathieu Equation
1.5 Wind Turbine Drivetrain . . . . . . . . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

1
1
2
8
9
12
14
17
20

CHAPTER 2 EQUATIONS OF MOTION FOR BLADE VIBRATION .
2.1 Formulation of In-Plane EOM . . . . . . . . . . . . . . . . . .
2.1.1 Modal Reduction of the EOM . . . . . . . . . . . . . .
2.2 Formulation of Out-of-Plane EOM . . . . . . . . . . . . . . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

22
22
27
30

CHAPTER 3 FORCED NONLINEAR MATHIEU EQUATION
3.1 Perturbation Analysis: First Order . . . . . . . . . . . .
3.1.1 Non-Resonant Case . . . . . . . . . . . . . . . .
3.1.2 Superharmonic Resonances . . . . . . . . . . . .
3.1.2.1 2Ω ≈ ω . . . . . . . . . . . . . . . . .
3.1.2.2 3Ω ≈ ω . . . . . . . . . . . . . . . . .
3.1.3 Subharmonic Resonances . . . . . . . . . . . . .
3.1.3.1 Ω ≈ 2ω . . . . . . . . . . . . . . . . .
3.1.3.2 Ω ≈ 3ω . . . . . . . . . . . . . . . . .
3.2 Parameter Dependence Discussion . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

35
35
37
37
37
39
41
41
42
43

CHAPTER 4 SECOND-ORDER PERTURBATION ANALYSIS . . . . . .
4.1 Second-Order Perturbation Analysis – Hard Forcing . . . . . . . . .
4.1.1 Superharmonic Resonances at ω/3 . . . . . . . . . . . . . .
4.1.2 Simulations and Discussions . . . . . . . . . . . . . . . . .
4.2 Weak Excitation—Primary Resonance Amplification . . . . . . . .
4.2.1 Weakly Forced Linear Mathieu Equation: Resonance Cases .
4.2.2 Case 1: No Resonance at O() . . . . . . . . . . . . . . . .
4.2.3 Case 2: Primary Resonance and Parametric Amplification .
4.2.3.1 Parametric Amplification . . . . . . . . . . . . .
4.2.3.2 Stability of Primary Resonance . . . . . . . . . .
4.2.3.3 Summary of Primary Parametric Resonance . . .
4.2.3.4 Cautionary Note on Second-Order Expansions in 

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

47
48
50
53
56
57
57
59
61
65
67
67

viii

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

4.3

4.4
4.5

4.2.4 Case 3: Subharmonic Resonance of Order Two . . . .
Effect of a Direct Forcing Term F0 . . . . . . . . . . . . . . .
4.3.1 Primary Resonance . . . . . . . . . . . . . . . . . . .
4.3.2 Subharmonic Resonance . . . . . . . . . . . . . . . .
4.3.3 Superharmonic Resonance . . . . . . . . . . . . . . .
Effect of Nonlinearity on Mathieu Resonances and Instabilities
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 5 APPLICATION TO WIND TURBINE BLADES
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . .
5.2 Numerical Simulations of ODE’s . . . . . . . . . . . .
5.3 Aerodynamic Forces on the Blade . . . . . . . . . . .
5.3.1 Sandia 100 m Blade . . . . . . . . . . . . . . .
5.4 Response of Wind Turbine Blades . . . . . . . . . . .

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

68
69
73
75
75
76
84

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

85
85
85
86
90
93

CHAPTER 6 DRIVETRAIN DYNAMICS . . . . . . . . . . . . . . . . . . .
6.1 Modeling of the Gearbox Reliability Collaborative (GRC) Gearbox . .
6.2 Case Study - Influence of Carrier Bearing Clearances on Gear and
Alignment and Stresses . . . . . . . . . . . . . . . . . . . . . . . . .
6.2.1 Time-varying Misalignment and Stress . . . . . . . . . . . . .
6.2.2 Optimizing the Design to Improve Contact, Stress and Life . .
6.2.3 Redesigning to Reduce Sensitivity . . . . . . . . . . . . . . .
6.3 Prediction of Wind Turbine Gearbox Vibration . . . . . . . . . . . . .
6.4 Redesigning the Wind Turbine Gearbox: For Improved Reliability . .
6.4.1 Tapered Roller Bearing for Carrier and Planet . . . . . . . . .
6.4.2 Interference Fit Planet Pins . . . . . . . . . . . . . . . . . . .
6.4.3 Decreased Spline Crowing . . . . . . . . . . . . . . . . . . .
6.4.4 Other Changes and Recommendations . . . . . . . . . . . . .
6.5 Torsional Model of Wind Turbine Drivetrain . . . . . . . . . . . . . .
6.6 Extension to Loading from Blade Vibration . . . . . . . . . . . . . .

. . . . .
. . . . .
Bearing
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .

. . 100
. . 100
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

104
105
105
107
109
111
112
114
115
116
117
117

CHAPTER 7 CONCLUSION AND FUTURE WORK . . . . . . . . . . . . . . . . . . . 119
7.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
7.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

ix

LIST OF TABLES

Table 4.1: Stability Wedge and Resonance Chart. R1 : Resonance identified at 1st order
of MMS expansion. R2 : Resonance identified at 2nd order of MMS expansion.
W2 : Instability wedge (Arnold tongue) expression can be found at second order
of MMS expansion. −: Known resonance case/ Instability not uncovered up
to two orders of expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
Table 5.1: Blade fraction align the station, chord length,twist and airfoil description for
the Sandia 100m blade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
Table 5.2: 100m blade structural properties . . . . . . . . . . . . . . . . . . . . . . . . . . 97
Table 6.1: Validation of gearbox housing natural frequencies . . . . . . . . . . . . . . . . . 102
Table 6.2: Metrics on the Quality of Gear Contact . . . . . . . . . . . . . . . . . . . . . . 105
Table 6.3: Metrics for the quality of gear contact for cases shown in figure 6.9 . . . . . . . . 108

x

LIST OF FIGURES

Figure 1.1: External loads on an off shore wind turbine (graphic from Jason Jonkman, NREL)

3

Figure 1.2: Global cumulative installed wind capacity from 1996 - 2013 [1] . . . . . . . . .

4

Figure 1.3: Electricity generation mix during 2013 in United States [2] . . . . . . . . . . .

4

Figure 1.4: Failure rate of wind turbine components and down time per failure [2] . . . . .

5

Figure 1.5: Large static and dynamic deflections of wind turbine blade (Images taken at
NREL) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

Figure 2.1: The arbitrary configuration of the position of a point P on the rotated and
deflected beam. The point P corresponds to location x on the undeflected beam. 23
Figure 2.2: The arbitrary configuration of the position of a point P on a beam rotating
in-plane and deflecting in the out-of-plane direction. The point P corresponds
to location x on the undeflected beam. . . . . . . . . . . . . . . . . . . . . . . 31
Figure 2.3: The nonlinear foreshortening due to out-of-plane blade motion at a point
along the blade at a distance x . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
Figure 3.1: Simulated response of equation (3.1) showing superharmonic resonances at
orders 1/2 and 1/3; µ̂ = 0.5,  = 0.1, α = 0.5, F = 0.1, γ = 3. The frequency
ratio sweeps up. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
Figure 3.2: Simulated response of the linear case of equation (3.1) showing superharmonic resonances at orders 1/2 and 1/3; µ̂ = 0.5,  = 0.1, α = 0, F = 0.5, γ = 3 . 40
Figure 3.3: Amplitudes of simulated response of equation (3.1) showing the effect of
the parametric forcing amplitude, with  = 0.1, µ̂ = 0.5, α = 0, F = 0.5.
Different curves depict Îł = 0.5, 1 and 3. . . . . . . . . . . . . . . . . . . . . . 44
Figure 3.4: Amplitude of simulated response of equation (3.1) showing the effect of the
direct forcing amplitude, with  = 0.1, µ̂ = 0.5, α = 0, γ = 3. Different
curves depict F = 0.5, 1 and 2. . . . . . . . . . . . . . . . . . . . . . . . . . . 45
Figure 3.5: Amplitude of simulated sweep-up response of equation (3.1) showing the
effect of the nonlinear term, with  = 0.1, µ̂ = 0.5, γ = 3, F = 0.5. Different
curves depict Îą = 0.1, 0.5 and 1. . . . . . . . . . . . . . . . . . . . . . . . . . 45

xi

Figure 4.1: Numerical evaluation of equation (4.14) shows the amplitude versus detuning
curve. Graph shown is generated using Âľ = 0.25,  = 0.1, F = 0.5, Îł = 3,
and θ = 0. amax = 0.2278 occurs at σ = −0.07. . . . . . . . . . . . . . . . . . 54
Figure 4.2: Simulated response of the linear case of equation (4.1) showing multiple
superharmonic resonances ¾ = 0.05,  = 0.1, F = 0.5, γ = 7, and θ = 0. For
the isolated superharmonic resonance at Ω/ω ≈ 1/3, amax = 2.477 by (4.17). . 55
Figure 4.3: Simulated response of (4.1) showing subharmonic resonance at 2; Âľ =
0.05,  = 0.1, F = 0.5, θ = 0, and γ = 0.5. (Low parametric excitation
is chosen so that the system doesn’t show unstable response at primary resonance. For the given set of values, system is unstable at subharmonic resonance.) 55
Figure 4.4: Variation of maximum amplitude of vibration at primary resonance as a
function of parametric excitation γ for various values of σ. (Left) Excitation
frequency is detuned to be slightly less than natural frequency. (Right)
Excitation frequency is detuned to be slightly more than natural frequency.
Parameters are F = 5, Âľ = 0.1 and  = 0.1. . . . . . . . . . . . . . . . . . . . . 63
Figure 4.5: (Left) Numerical simulations show that the predicted parametric resonance
can have errors when the amplification becomes large as Îł reaches critical
values. (Right) Graphed is zoomed in to show that numerical and analytical
results match for moderate levels of parametric excitation. . . . . . . . . . . . . 63
Figure 4.6: (Left) The numerical evaluation of equation (4.34) showing the amplitude
versus the detuning. Graph shown was generated using Âľ = 0.25,  = 0.1, F =
5, Îł ranges from 0 to 4 in 0.5 increments; as Îł increases a increases. (Right)
The numerical evaluation of equation (4.34) for Îł = 3 showing the amplitude
versus the detuning at various orders of  expansion. Note that the curves
with O() (dotted line) and O( 2 ) (dashed line) expansions are visibly erroneous. 64
Figure 4.7: Stability wedge for the linear forced Mathieu equation with damping at primary resonance ((ω +  σ)2 = δ = N 2 /4, N = 2 for traditional Mathieu
equation). Figure shows comparison of stability of numerical system response/amplitude and the analytically calculated stability boundary. . . . . . . . 66
Figure 4.8: Amplitudes of simulated responses of equation (4.45) showing unstable response at subharmonic resonance due to increase of the parametric forcing
amplitude;  = 0.1, Âľ = 0.25, F0 = 2. Different curves depict Îł = 0.5 and 1. . . 72
Figure 4.9: Amplitudes of simulated responses of equation (4.45) showing the occurrence of superharmonic resonances at high parametric forcing amplitude;
 = 0.1, Âľ = 0.25, F0 = 2. Different curves depict Îł = 2 and 4. . . . . . . . . . . 72

xii

Figure 4.10: Trend of amplitudes of simulated responses of equation (4.45) at different
critical frequencies. Order 2 subharmonic goes unstable at significantly lower
values of Îł and hence is not shown. . . . . . . . . . . . . . . . . . . . . . . . . 73
Figure 5.1: Response amplitude of system described by equation (2.8) that described inplane motion of the blade. Values used in simulation b = 0.1, d = 0.1, e1 =
1, e2 = 0.3, f1 = 0.1, f2 = 0.1, F = 1. To account for damping we add a 2µqÛ
term with Âľ = 0.25 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
Figure 5.2: The local forced on an aerosol structure of the blade, reproduced from [3] . . . 88
Figure 5.3: Velocities at the rotor plane, reproduced from [3] . . . . . . . . . . . . . . . . . 89
Figure 5.4: Screenshot of NuMAD blade geometry, reproduced from [4] . . . . . . . . . . 90
Figure 5.5: Variation of coefficient of lift and drag for the NACA 64 airfoil and DU30 airfoil 92
Figure 5.6: Amplitude response of the NREL GRC blade, blade properties from [5] . . . . 94
Figure 5.7: Response of a linear system with cyclic loading (left) and constant loading
(right), based on parameter from a 60 m blade. . . . . . . . . . . . . . . . . . . 95
Figure 5.8: Ratio of parametric excitation to the square of natural frequency as a function
of length of the blade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
Figure 5.9: Flutter speed as a function of blade length. Data obtained from [6] . . . . . . . 99
Figure 6.1: RomaxWind model of GRC gearbox . . . . . . . . . . . . . . . . . . . . . . . 101
Figure 6.2: Schematic showing GRC gearbox layout and component nomenclatures . . . . . 101
Figure 6.3: Comparing the 1st natural frequency and modal response of the gearbox
housing. Figure shows the test set-up and test result (319 Hz) on the left and
RomaxWind model and simulation result (338 Hz) on the right . . . . . . . . . 103
Figure 6.4: RomaxWind model of GRC gearbox installed in wind turbine with a load case
example: Input Torque (Mx) = 325 kNm, Off Axis Moment = -100 kNm,
Rotor Weight = -1500 kN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
Figure 6.5: Planet load share in kNm vs. carrier rotation in degrees . . . . . . . . . . . . . 106
Figure 6.6: Sun/Planet gear misalignment in microns vs. carrier rotation in degrees . . . . . 106

xiii

Figure 6.7: Results of full factorial study for upwind and downwind planet carrier bearing’s radial internal clearance (shown as PLCA RIC and PLCB RIC in figure):
(left) Sun/Planet misalignment; (right) Sun/Planet load distribution. . . . . . . . 107
Figure 6.8: Sun/Planet load distribution for (top) original case; (bottom) case with tapered
roller bearing with +100 Âľm preload and after micro-geometry optimization . . 109
Figure 6.9: Downwind planet bearing maximum contact stress with carrier rotation for
(top) original case; (bottom) case with tapered roller bearing with +100 Âľm
preload and after micro-geometry optimization . . . . . . . . . . . . . . . . . . 110
Figure 6.10: Normalized amplitude of vibration at different sensor location for High speed
shaft (HSS) and Intermediate speed shaft (IMS) gear excitations . . . . . . . . . 112
Figure 6.11: RomaxWind model with carrier and planet gear tapered roller bearing in “X"
(marked at location of bearing by rectangles) and “O" (marked at location of
bearings by circle) configurations . . . . . . . . . . . . . . . . . . . . . . . . . 113
Figure 6.12: Mesh misalignment vs. preload for initial tapered roller bearing selection . . . . 114
Figure 6.13: Mesh misalignment for increasing ring and carrier stiffness with interference
fit pins for planetary gear contacts . . . . . . . . . . . . . . . . . . . . . . . . . 115
Figure 6.14: Deflection distribution on tooth on horizontal axis for decreasing crown of
(top left) 173Âľm, (top right) 86.5Âľm, (bottom left) 43.25Âľm, (bottom right)
21.625Âľm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
Figure 6.15: 1-D torsional model of a simple wind turbine gearbox, courtesy of Romax
Technology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118

xiv

CHAPTER 1
INTRODUCTION

1.1

Objective

The aim of this work is to improve the understanding of wind turbine blade oscillations and
how they translate to dynamic loading on the gearbox. Blade oscillations can be generated in a
variety of ways. Oscillations in the blades imply dynamic loadings transferred to the gearbox.
Understanding these dynamic loadings is essential to the design of reliable gears and bearings, and
hence economically viable wind turbines.
In particular, this work focuses on horizontal axis wind turbines with aerodynamic and gravitational loading, and aims to
1. Obtain the nonlinear dynamical equations of motion of rotating blades, including nonlinearity
due to large deflections (nonlinear curvature and nonlinear foreshortening), aerodynamic and
aeroelastic loading, and gravitational loading.
2. Reduce the partial differential equations with assumed modes. Analyze the reduced model
using perturbation methods.
3. Understand the resonances and instabilities for the developed model. Isolate the roles of
parameters and each type of loading on the dynamics.
4. Understand the dynamics for single blade models with in-plane deflection and out-of-plane
deflection.
5. Interpret the resonant and unstable dynamics as input loads to the gearbox. Simulate gearbox
loading with commercial gearbox software (Romax) to understand the sensitivities of bearing
loads and gear loads for different operating conditions.

1

The approach in this thesis is to apply reduced-order modeling, and obtain analytical expressions
of response characteristics as functions of parameters, so that the roles of parameters may be
understood and used to guide more detailed numerical simulations. The numerous sources of
external loads for an example offshore wind turbine blade are shown in figure 1.1 to showcase the
complexity of these systems.

1.2

Background and Motivation

The wind industry has grown consistently at high rates every year in terms of global cumulative
installed capacity (figure 1.2). This is enabled by political and public pressure behind the wind
industry to boost this renewable energy sector with ambitious plans, such as 20% wind energy by
2030 in the USA [7] and similar continuously increasing targets in other global markets. As of
2013, wind energy accounts for 4% of the total energy generation mix in the United States (see
figure 1.3).
Accompanying the growth in demand and installed capacity is the growth in the size of the
machinery. In 2001, the global average size of installed turbines was 0.93 MW and by 2008 this
metric had reached 1.53 MW [8]. Current production wind turbines are typically in the 2 - 3 MW
range and active research and development is being done for offshore installation with a capacity of
5 - 13.5 MW per wind turbines. Alongside this rapidly expanding market and growing machinery
is a significant issue with the reliability and maintenance of rotating machinery, and there is a focus
on improving gearbox and bearing life. With the larger wind turbines, gearbox failures are too
frequent.
Reliability statistics for European databases of wind turbine failures and the analysis of the
Landwirtschaftskammer Schleswig-Holstein (LWK) [9] which survey up to 643 wind turbines is
reproduced here from the American Wind Energy Association’s annual report in figure 1.4. The
statistics show the blades and drivetrain (gearbox and main shaft) have significant failure rates and
also a gearbox failure yields the highest downtime. It is also one of the most expensive components
to replace; with any failure comes additional costs for removal and replacement, and loss of energy

2

Figure 1.1: External loads on an off shore wind turbine (graphic from Jason Jonkman, NREL)

3

Figure 1.2: Global cumulative installed wind capacity from 1996 - 2013 [1]

Figure 1.3: Electricity generation mix during 2013 in United States [2]

4

Failure/turbine/year and downtime from two large surveys of land-based European wind turbines over 13 years
Electrical System

LWK Failure Rate, approx 5800 Turbine Years

Electrical Control

WMEP Failure Rate, approx 15400 Turbine Years

Other

LWK Downtime, approx 5800 Turbine Years

Hydraulic System

WMEP Downtime, approx 15400 Turbine Years

Yaw System
Rotor Hub
Mechanical Brake
Rotor Blades
Gearbox
Generator
Drive Train
1

0.75

0.5

0.25

0

2

4

6

8

10

12

14

Downtime per failure (days)

Failure/turbine/year

Figure 1.4: Failure rate of wind turbine components and down time per failure [2]
capture. These findings have also been substantiated with data collected by other organizations
and
12
researchers (for example, [10]).
The growing market and need for improved reliability and robust designs have brought an
industry-wide efforts such as the development of a new IEC standard 61400-4, and the National
Renewable Energy Laboratory Gearbox Reliability Collaborative (NREL GRC) [5, 11], a large
multi-year project which brings together a range of participants across the industry for analysis
and testing of wind turbine gearboxes. Germanischer Lloyd, a certification body enlists a set of
design specifications [12] which designers should adhere to, in order to get the turbines certified.
In spite of these efforts, frequent component failures occurring in wind turbine systems well before
expected lifetime of operations indicate that loadings and failure mechanisms of operating wind
turbines are not yet sufficiently understood. The lack of understanding leads to potentially operating
the system under extreme loading conditions that have not been accounted for. It is likely that there
exist unknown dynamical loading cases.
There are a variety of dynamic loading cases. The turbine during steady operation has cyclic
aeroelastic and gravitational loads. The aeroelastic loads have the capacity to cause self-excitation,
much like the aeroelastic flutter observed in aircraft [13, 14]. There are sudden changes of wind
5

speed and direction, which lead to time varying aeroelastic blade loadings. The ratio of flutter
speed to wind turbine operating speed also reduces as the blades grow in size. In [4], it is noted
that as blade length increases from 5 m to 34 m the ratio drops from 6 to about 2.5. For much larger
blades, there is little to no margin on flutter speed based on some preliminary calculations [6].
When winds become high, the brakes may be applied, introducing the possibility of frictional
stick-slip, which has a myriad of complications in and of itself [15,16]. A yaw controller can induce
yaw motions, which in turn induce inertial blade and shaft excitations. It has been noted that parked
wind turbine system also exhibit unstable dynamic behavior that cause extreme loading on the blades
and contribute to failure [17]. When the turbine is on idle, there can be self-induced oscillations
in the standby position, which can then involve stick-slip, and clearance and impact events in the
gearbox. The dynamic blades are coupled to other substructures. Coupling with the gearbox allows
the gearbox modes and nonlinearities to affect the rotor and blades. Coupling between two or three
blades, each of which has nearly the same modal characteristics, through the blade hub, allows for
the possibility of internal resonances [18] and symmetry-breaking phenomena [19], which can lead
to vibration localization [20, 21], meaning extra-large vibrations locally on the coupled structure.
Finally, the total structure is also coupled with the tower, which may have nearly identical in-plane
and out-of-plane modal characteristics. Any resonances or instability can amplify the loading that
is transferred to the gearbox. It would be no wonder, if there are such dynamical events that are
currently not catalogued, that drivetrain loads go beyond their design specifications.
The large wind-turbine industry will benefit greatly from improved understanding of the dynamical loads on the gearboxes. To better understand these dynamical loads, a careful description
of the aeroelastic and gravitational loads is required, and then a careful nonlinear dynamic analysis
for instabilities and resonances is needed.
In this thesis, we will focus on the role of cyclic gravitational and aerodynamic loading on the
in-plane and out-of-plane vibration of an isolated blade. We also include nonlinearity associated
with large deflections (nonlinear curvature and nonlinear foreshortening). Very large blades exhibit
large deflections under static loads and dynamic loads (see figure 1.5). Dynamic deflections have

6

the capacity to be even larger. The blades are more flexible in the out-of-plane direction with
deflections up to 10 m during normal operation.
Accounting for nonlinearity will enable us to accurately capture the motion of the blade during
operation under specified loads and will also expose the presence of resonances other than one-toone resonances of the structural modes. Combining nonlinearity and parametric excitation terms
leads to subharmonic and superharmonic resonances. That is, such conditions can easily lead to
resonances when the excitation frequency is twice or half (or some multiple or fraction) of a modal
frequency. We first analyze the loadings on a single rotating blade.
Aeroelastic blade loadings are known to induce self-excited blade oscillations [22–26], which
tend to involve the coupling of particular linear modes, and get particularly more serious with
offshore turbines [27]. In horizontal axis wind turbines, aeroelastic loadings can be cyclic, thereby
producing forced oscillations of the blades [28]. Sources of cyclic variation are the wind speed
gradient, or wind shear, such that the blade loading varies with elevation, as well as tower pass-by.
Furthermore, the blades endure cyclic gravitational loadings which also lead to forced oscillations.
These cyclic and aeroelastic loadings can induce vibrations and instabilities. We aim to understand
these dynamic responses.
There is a large body of work on modeling of rotating blades under aerodynamic and gravitational loading. Wendell [29] developed the partial differential equations of motion for a rotating
wind turbine blade. Additional efforts have been made by to include gravity and pitch action by
KallosĂłe [30]. Caruntu [31] has developed nonlinear equations that model the flexural potential
in non-uniform beam that can be extended to large flexible blades of wind turbines. Chopra et
al. [32] have modeled the blades as a hinged structure and constructed equations of motion to
study its dynamic behavior. Jonkman [33, 34] used the modal information of blade vibration to
calculate the tip displacement in blade, the torques and moments experienced at the hub amongst
other things. These calculations have been used to develop the code FAST that is used in research
laboratories, industry and universities around the world as a basis for performance evaluations. Numerous commercial and research codes are also available that enable us to understand the behavior

7

Figure 1.5: Large static and dynamic deflections of wind turbine blade (Images taken at NREL)
of wind turbine blades and the influence of external loads on their responses. Recent academic
thesis [3, 35, 36] have addressed improving aerodynamic and aeroelastic modeling of wind turbine
blade in order to improve formulations of the loading of wind turbine systems.
Extensive literature also exists on analysis related to the problem of vibration of helicopter
blades [37–40]. Hodges and Dowell [41] also developed the nonlinear partial differential equations
for a twisted helicopter rotor blade. These rotate at much higher rotational velocities than wind
turbine blades, and in a horizontal plane, unlike the wind turbine blades. Hence the gravitational
influence on the vibration characteristics of the wind turbine blade would be much different.
However, many features of helicopter blade dynamics, such as modeling of centrifugal effects and
deflected geometry, carry over to the wind turbine blade models.

1.3

Wind Turbine Blade Modeling

Wind turbine blade modeling can be broadly classified into two categories: aerodynamic
modeling and aeroelastic modeling.

8

1.3.1

Aerodynamic Modeling

Theodorsen’s work [42] on the general theory of aerodynamic instability and mechanism of flutter
laid the foundations for significant follow up work in the area. The aerodynamic forced on an
oscillating airfoil with three degrees of freedom were determined. Numerous simplification and
alterations of these models have been used in commercial codes to predict the loads on an airfoil
structure. There are numerous other simpler aerodynamic models and a few that increase accuracy
by more detailed formulations. Some of the various types of aerodynamic modeling are listed
below:
1. Actuator disc model (momentum theory): The simplest model for wind turbine aerodynamics
that provides useful results and insights. In this, we consider a rotating annular stream tube
and apply conservation of momentum and Bernoulli’s equation to get axial and tangential
forces on annular element of fluid.
2. Glauert annulus momentum vortex theory: The next level of complexity is the annulus
momentum vortex theory. Each radial section of the blade is analyzed independently using
a 2-D airfoil data and equations are based on continuity and conservation of momentum.
Corrections based on the work of Prandtl and Goldstein have been satisfactory for extending
this theory to finite number of blades (average effects). This theory is most useful for basic
rotor design.
3. Prescribed wake vortex theory: The basic concept is that the vortices that form a spiral trail
behind the rotor-blade tips define the flow field around the rotor according to the Biot-Savart
theorem. Each section of the rotor blade generates a lift force that is proportional to its local
bound vorticity. By integrating the effects of trailing vortices over the blade, the induced
flow and hence the rotor forces and moments can be determined.
4. Free wake vortex theory: The assumption of rigid vortex theory can be eliminated and one
can attempt to find the exact path of the trailing vortices iteratively.

9

The blade element momentum (BEM) theory is a combination of the actuator disc model
(momentum theory) and the annulus momentum vortex theory (blade element theory). It is a
computationally cheaper method based on the 2D airfoil characteristics (lift and drag). BEM
theory is most commonly used in wind turbine dynamic simulations due to its computational
efficiency and overall accuracy. Limited success has been made with computational flow solvers
based on Reynolds Averaged Navier Stokes (RANS) and other similar three-dimensional models.
This is primarily due to the complexity of modeling wind turbines. Wind turbine aerodynamics
are dependent on far field conditions, several rotor diameters up and down stream, while at the
same time being dependent on small scale flow conditions at the blade. Coupled with body motion,
the need to have fine resolution and a large domain makes these models highly computationally
intensive.
BEM theory equates two methods of examining how a wind turbine operates.
1. The first method is to use a momentum balance on a rotating annular stream tube passing
through a turbine.
2. The second method is to examine forces generated by the aerofoil lift and drag co-efficient
at various sections along the blade (in an undeformed state) and these forces are integrated
over the total length of the blade to give the overall aerodynamic loads.
BEM theory is used to develop a blade design tool by considering power requirements [43]. The
rotor torque applied by the lift and drag forces (calculated using BEM) are determined and then this
torque is equated to the torque applied by rotating annulus (using Bernoulli’s equation) produced
by wind field. No aerodynamic interactions between different blade elements are considered and
the blades are assumed to be rigid.
In [44–46], BEM theory is used to optimize the blade shape and pitch control strategy. A cost
function is formulated based on the aerodynamic performance of the blade (power obtained from
lift and drag forces), weight of the blade and its rigidity. Aerodynamic properties are calculated
using BEM and blade rigidity is determined using a simple beam theory. A simplified pitch control
10

system has been formulated to relate pitch angle with power. This cost function is then used to
optimize the blade parameters like chord length, section area and twist.
In [47], BEM theory is used in integration with a Finite Element (FE) approach. The system
model is constructed by combining an aerodynamic model (obtained using BEM) and a FE blade
model. Aerodynamic model consists of a BĂŠzier surface to calculate the lift, drag and moments
caused by aerodynamics for different wind velocities and attack angles. In the model, the wind
velocity is assumed to be normal to the plane of rotation and steady state flow is considered.
Variation of the gust excitation is modeled by varying the angle of attack. Change in aerodynamic
forces due to dynamic stall (dynamic stall is the unsteady aerodynamic effect caused by the strong
vortex formation that occurs when airfoil rapidly changes the angle of attack) is neglected. The
FE model of the blade consists of several sections and each section is made up of 3D elastic beam
elements. Connection between the blade and hub is rigid and the twist angle variation of the blade
is neglected. The blade section is assumed to be piece-wise variable along its length.
BEM theories generally give acceptable approximations to aerodynamic loads under the condition where wind is normal to the plane of rotor (i.e. the turbine is un-yawed with respect to the
wind). But under yawed conditions (when wind is at certain angle to the plane of rotation) the
BEM theory is less able to model the physics due to the 3D nature of the flow.
Since the basis of the BEM theory is strictly a 2D analysis, 3D effects such as physical roll-off in
the lift as the blade tip is approached (tip loss) cannot be modeled analytically but can be treated by
using some corrections such as Prandtl’s tip loss function. Several advanced aerodynamic models
are described in [48] which take into account the 3D wind flow properties.
The model described in [48] takes into account the vertical wake system produced behind the
rotor. This wake is comprised of vortices that trail from each of the blade tip. Though various
modifications to the BEM approaches have been used to represent these wake effects, they have
not been satisfactory in determining the blade loads. An improved version models these wakes
by considering the aerodynamic lag of the inflow development over the rotor disk in response to
changes in blade pitch or rotor thrust. The mathematics is described in terms of ordinary differential

11

equations and a time constant representing the dynamic lag in the build-up of the inflow. Since the
aerodynamic equations are described in terms of ordinary differential equations, it can be easily
coupled with the aeroelastic models and the combined system can be solved simultaneously.
Another model described in [48] takes into account the unsteady aerodynamics caused by
the blade section (stall effect) and tower shadow effect (the distortion in the wind flow when the
blade passes through the tower). The model used is semi-empirical but is represented in terms of
state-space formulation and hence is easy to implement. Other advanced aerodynamic models are
based on solution of full field Navier-Stokes equation. But these methods are computationally very
expensive and are generally implemented in CFD codes.
Our work focuses on modeling the blade vibration and related dynamic instabilities. We seek
to use a simpliefied model to account for the effects of wind on the blade loading. BEM theory, in
spite of its limitations provides a valid model for blade loading due to wind flow. We incorporate
wind shear in our calculation of blade load as it induces a cyclic variation for every revolution.

1.3.2

Aeroelastic Modeling

Simple beam models (Euler-Bernoulli or Timoshenko) are used to model the blade sections in
[43, 46, 47]. These beam models take into account the blade flexibility. The degrees of freedom
involved are less, in order to make these models computationally cheaper. Typically, beam models
are all based on symmetric cross-sections and hence they neglect the effect of cross-sectional
properties like shear center, tension center, etc. Another draw back of these beam elements is the
assumption that an initially plain cross-section remains plain under deformation. This assumption
does not hold for wind turbine blades but nonetheless provide a basic foundation for studying the
dynamics of wind turbine blades.
In [30], as noted before, an advanced beam model is developed by using the principle of virtual
work (taking into account the potential and kinetic energy of the element). The equations of motion
are derived in terms of partial differential equations (using a Lagrangian formulation), relatively
small terms are neglected in a consistent manner by introducing an ordering scheme. The model
12

takes into account the effect of sectional properties like shear center and gravity center on the blade
motion. The effects of gravitational forces and Coriolis forces are also considered. The model has
full coupling between the bending and torsion modes of the blade. A pitch model is also developed
to determine the pitching torque for the given value of pitch angle and consequently the effect of
pitching action on the blade dynamics is also studied. An eigenvalue problem is solved using the
finite difference discretization of the equations of motion to get the natural frequencies and mode
shapes. However, the influence of yaw motion or tower flexibility on the blade dynamics is not
considered, and the blade is assumed to be isotropic and inextensible.
Other advanced rod models are formulated [49], that takes into account the warping effect, blade
extension and deformation of the blade section due to tilt. The equations of motion are derived in
terms of partial differential equations using the principle of virtual work. These equations are also
coupled with a warping function that takes into account the out-of-plane deformation of the cross
section. A strain-displacement relationship is also formulated that takes into account the distortion
of the blade section due to in-plane motion. Instead of using one element for a cross section; the
cross section is defined using several quadrilateral elements. Due to the higher number of degrees
of freedom, reduction methods are used to reduce the size of the model. The model is then validated
with the experimental results. In this work, the blade model is developed only to validate/justify
the improved rod element for WT applications, so the blade is assumed to be stationary (like a test
setup); therefore the effect of centrifugal forces and gyroscopic terms is neglected.
Full finite element blade models are used to investigate the modal response in [3]. But these
models are computationally very expensive. Modal analysis can be performed using these models
but it is very difficult to use them for full time domain dynamic analysis.
Other efforts have been made to do experimental modal analysis to determine the natural
frequencies, damping and mode shapes for blades [50]. The test setup is capable of testing full
sized wind turbine blades (up to 20 m in length). The blade is excited by a hammer and the
frequencies and mode shapes associated with first and second flap mode, first edge mode and first
torsional mode are calculated. Similar efforts have been made at other national labs across the

13

globe, however, with the increasing dimensions of the wind turbine blades, these methods are
restricted for use on small blades.
However, it should be noted that as wind turbine increase in size, their design will be stabilitydriven in contrast to current load driven design and analysis. Large rotating structures can have
harmonic resonances and instabilities in their operating range. Operating the wind turbine under
those parameters cannot be avoided, unless the design process accounts for existence of such
resonances and instabilities. In fact, parked (idle) wind turbine systems also have instabilities
associated with the edgewise motion of the blade and the side to side motion of the tower [17].
Hence a good aeroelastic code is essential to capture the inherent structural response properties of
wind turbine blades. We develop the structural model in this thesis, by defining the motion of a
single blade around the rotor at a constant speed. The variations in geometry along the blade are
accounted for. Large deflections induce nonlinear phenomena, however, the associated details on
the aerodynamic loads is less comprehensive, but, sufficient to capture the necessary dynamics.

1.3.3

Overview of Design Codes

Wind turbines systems are designed and analyzed using simulation tools that capable of predicting
the coupled dynamic response of a system under various operating conditions. Onshore wind
turbine analysis relies on the use of aero-servo-elastic codes, which incorporate wind inflow,
aerodynamic, control system (servo) and structural-dynamic (elastic) models in time domain in a
coupled simulation environment [33]. In recent years, a number of these codes have been expanded
to include the additional dynamics pertinent to offshore support structures. Some standard design
codes/softwares that are currently in use in industries that have been developed by NREL to calculate
loads and forces in wind turbine structures are listed below:
• AirfoilPrep
– Developed by Windward Engineering and NREL
– Generates airfoil data files from 2D data (experimental)
14

– Adjusts 2D data for rotational augmentation (3D effects) and computes dynamic stall
parameters
• WT Perf
– Developed by NREL (initially developed at OSU)
– Calculates performance of rotor geometry and airfoil data
– Steady-state calculations (no dynamics), computes power, torque, thrust and blade-root
bending moment
– Uses BEM theory
• PreComp
– Developed by NREL
– Computes coupled section properties of composite blades for beam-type models
– Inputs are blade external shape and internal lay up of composite laminas
– Section properties like direct and cross stiffnesses, inertias, locations of shear, tension
and mass centers are found directly
– Uses modified combined laminate theory with shear flow approach
– Underlying assumptions: straight blade, thin walled closed sections, no transverse
shearing, no in-plane distortion, no warping
• AeroDyn
– Developed by Windward Engineering and NREL
– Computes aerodynamic loads as part of aero-elastic solution (fully coupled)
– Uses BEM and Dynamic Wake Algorithms; hub loss, tip loss and skewed wake corrections in BEM
– Capable of simulating turbulent wind inputs
15

Apart from the methods and codes listed above, there are numerous other codes that have been
developed by national labs and universities for aerodynamic and aeroelastic modeling. FLEX5
by Technical University of Denmark, HAWC2 by Risø National Laboratory, DUWECS by Delft
University of Technology, Alcyone and GAST by National Technical University of Athens, PHATAS
by Netherlands Energy Research Foundation, to name a few have been developed in an effort to
understand various aspects of wind turbine loads. Commercial code such as ADAMS/WT developed
by MSC Software and GH Bladed by Garrad Hassan have also been extensively used in the industry
for predicting loads experienced by the blades during operation.
One limitation of these linear frequency and time domain analyses is that they cannot capture the
non-linear dynamic characteristics and transient events that are important considerations in wind
turbine analysis. Each code focuses on different aspects of the model and adopts different techniques
to assess loading. These can be broadly classified into three main categories: multi body simulation
using rigid elements (for e.g. ADAMS/WT), finite element modeling of the blade structure (for e.g.
HAWC2) and an assumed cantilever beam mode approach (for e.g. FAST). However, in our survey
of these codes, most of the codes assume that the deflection of the wind turbine blades are small
and consequently, linear. FAST applies foreshortening to the blade deflection, however, linearity is
still assumed. This assumption becomes increasingly invalid for large deflections of wind turbine
blades. Also, aerodynamic loads are applied to the static blade shape. This is a pretty valid
assumption for most cases, however, literature exists that account for aerodynamic loading with
deformed blade elements. The loads computed by the BEM theory suggested in this thesis should
be extended for local deflections of the blade element. Aerodynamic modeling using the ONERA
semi-empirical approach for dynamic stall provides a better framework for computing lift and drag
forces that are dependent on the instantaneous velocity of the blade accounting for motion due to
blade vibration. In [51], Bierbooms compares two semi-empirical approaches and mentions that
the unsteady aerodynamic effects captured by the ONERA approach may have important influence
on the fatigue life of the wind turbine and its stability.
Another aspect of code development for wind turbine loads analysis has been that blade modeling

16

and flow modeling theories have been developed typically in isolation. No model implements the
advanced aerodynamic models in conjunction with the advanced blade models. Some investigation
is also required on the effectiveness of the existing theories to model the complex non-linear
phenomenon like warping, tower shadow and dynamic stall accurately. There is very little work
done on the influence of tower flexibility and yawing motion on the blade dynamics, although some
very simple models do exists that attempt the capture these interactions.
The need to improve the modeling mechanisms of blade and wind turbine systems has generated
a lot of interest in the academic community and consequently various researchers have developed
mathematical models to best represent the wind turbine blades and drivetrains. Sant [36] explored
improving the BEM-based aerodynamic models by including the free wake vortex formulations.
Weinzierl [35] also developed a BEM-based simulation tool that incorporates active flow control
elements (for e.g., pitch control, trailing edge flaps, etc.) in the model. Larsen [52] proposed a
dynamic stall model for capturing the forces on the wind turbine air foil. Hansen and Butterfield [53]
in their review paper make note of the improved understanding of dynamics loads for HAWT’s
whose structural dynamics and airfoil characteristics are well understood. However, the prediction
of loads under off-design (say superharmonic resonant loading) conditions still remains a formidable
challenge. Blade motion under resonance has also been observed in simulations and experiments
of out-of-plane motion in work done by Ishida et al. [54, 55].
In this thesis, we model the blade structure with large nonlinear deflections through assumed
mode approximations of a beam and apply aerodynamic loads using the BEM theory. This
allows for assessing dynamic response phenomenon that have not been captured by simple linear
representations of blade motion and accounts for aerodynamic loading.

1.4

Beam Vibrations and the Mathieu Equation

Modeling and analysis of the dynamics of beam is one of the cornerstones towards improved
understanding of the behavior of various systems ranging from the stability of bridges and buildings
to nematode locomotion. Mathematical expression that detail the static and dynamic response of

17

the beam under various loading conditions enables us to explain the fundamental phenomenon that
govern the response of the more complex structures they represent under similar loading conditions.
In the presence of nonlinearities, the beam response exhibits numerous phenomenon not observed for the linear system. Harmonic resonances, hysteresis, chaos and numerous other phenomena can be characterized by the nonlinear equations of motion that are used to represent the motion
of a beam under excitation. Nonlinearities occur due to numerous factors and as such in reality
all physical system exhibit nonlinear behavior. Typical nonlinearities considered for modeling the
vibration of beams include nonlinearity arising due to inertia, geometry, large deflection and damping [56]. The manifestation of nonlinearity in the more complex systems that the beam represent
can be due to other factors, for e.g., in gear tooth modeling, backlash and deformations in the root
induce nonlinearity which can be captured via a beam model. In the formulations of nonlinear
analysis of beam motion, it is typical to consider weak nonlinearity, and use book keeping parameters to quantify relative magnitude of the terms. This help us arrive at an approximate analytical
solution to the differential equations, thereby, enabling us to understand the influence of terms that
govern the dynamics of the system.
There have been numerous studies on different formulations for modeling blade vibration.
Yao [?], Sabuncu [57], Venkatesan [58], Yoo [59] and Al-Nassar [60] have analyzed the dynamics
of rotating blades with some additional complexity viz., nonlinear flapping motion, pre-twisted
blades, varying rotational speed, stability due to torsional excitation, etc. These studies have
been motivated by the fundamental dynamics that govern the system under different conditions.
The blade models with nonlinearity and excitation are represented by the Duffing oscillator or
the nonlinear Mathieu equation with parametric excitation. Both these are fundamentally well
understood system, but, ongoing research in this area still actively explores these systems and their
characteristic responses.
In this thesis, we use the beam to model the motion of the wind turbine blade. Under steady
operation and wind conditions, in which the blade rotates at a constant rate Ω, a blade takes on cyclic
transverse loading due to both gravity and aerodynamic forces with effects of wind shear and tower

18

passing. The blade also endures gravitational axial compression (softening) when pointing up, and
axial tensioning (stiffening) when pointing down. Considering these effects on a single-degreeof-freedom description of transverse deflection motivates us to study a forced Mathieu equation
with cubic nonlinearity. Also, the terms intrinsic to a forced nonlinear Mathieu is embedded in the
in-plane and out-of-plane EOM developed for the wind turbine blade.
There have been extensive studies on systems with parametric excitation that fit in into a
minor variation of the Mathieu equation. Rhoads and Shaw [61] have studied MEMS structures
with parametric amplification, in which direct excitation occurs at half the frequency of excitation
of the parametric excitation. Parametric amplification was also demonstrated in experiments by
Miller [62]. Other work has examined nonlinear variations of the Mathieu equation, including
Van der Pol, Rayleigh oscillators and the Duffing nonlinear terms. Rand et al. [63–66] analyzed
the dynamics and bifurcations of a forced Mathieu equation and properties of superharmonic
2:1 and 4:1 resonances. Belhaq [67] studied quasi-periodicity in systems with parametric and
external excitation. Veerman [68] analyzed the dynamic response of the Van der Pol - Mathieu
equation. Arrowsmith [69] and Marathe [70] studied the stability regions for the Mathieu equation.
Reference [71] is an extensive compilation of the various fields of study in which the characteristics
of the Mathieu equation are found and employed.
In studying the Mathieu, we omit effects of other terms specific to the rotating beam model,
as our attention is drawn to this more fundamental equation in dynamical systems. Indeed, full
partial differential equations (PDEs) of the rotating beam can be analyzed and modally reduced.
The analysis of the full PDE is pursued through simulations to demonstrate that the underlying
dynamic phenomenon is that of the forced nonlinear Mathieu equation.
Analytical solutions of the Mathieu equation is performed using perturbation techniques, such
as the method of multiple scales (MMS) [72–75]. We analyze resonances of the forced nonlinear
Mathieu equation, as a phenomenological representative of blade motion. The analysis reveals the
existence of subharmonic and superharmonic resonances. We unfold the super- and sub-harmonic
resonance cases present in the system and identify those that may be critical to wind turbines during

19

normal operation.

1.5

Wind Turbine Drivetrain

The rotor hub is the mechanism that is used to attach the wind turbine blades to the drivetrain.
The wind turbine drive train is a huge single speed gearbox with speed multiplication in the range
of 100:1. Aerodynamic forces from the wind are translated to rotational forces which amount to
torque at the generator. This is the ‘useful’ work done by the wind turbine system. All the other
loads and moments (other than the rotational torque) are transferred to the drivetrain and tower.
Typical wind turbine design codes develop the structural model with about 16 to 24 degrees
of freedom [76]. Most of the dof’s are associated with the bending and torsional modes of the
blades and tower and the entire drivetrain is represented by 1 dof. This simplification implies that
the individual elements that make up the drivetrain do not contribute to the overall dynamics of
the system under study. This is however not true and only proceeds to mask our understanding
of the interaction between the dynamics of the blade and the components in the gearbox. The
coupling between the external low frequency excitations from the blade dynamics with the internal
high frequency dynamics of the gearbox (such as gear mesh and bearing frequency) could produce
undesirable loads. Also, as wind turbine blades and gearboxes increase in size the frequency gap
between the internal resonances in the gearbox and external system level resonance decreases,
necessitating a more detailed formulation of gearbox and drivetrain in the modeling environment
for wind turbine systems. There has been extensive efforts by National Renewable Energy Labs
and Romax Technology through the Gearbox Reliability Collaborative to improve industries understanding of the complexity of the loading conditions in the gearbox. This effort treats loading from
blades as external. Our focus is to integrate the blade and gearbox simulations as a path way for
future studies that will recognize the need for detailed system level models for design and analysis
of large wind turbine systems. Some preliminary modeling efforts performed as a part of this thesis
include assessment of gearbox design. In particular, detailed analysis has been performed to assess
the impact of bearing clearances on gear and bearing misalignments and stresses. A simulation

20

toolset is also used to predict the vibration of the gearbox housing under various operating conditions. Analysis from the models developed for the studies was used to perform a comprehensive
assessment of factors that influence the reliability of a wind turbine gearbox. Some critical parameters include the choice of bearings, manufacturing fits, and other detailed geometry modification of
gearbox internals. RomaxWIND provides solutions during steady state operation of wind turbine to
provide damage and life assessment of the drivetrain. We also explored quasi-transient approaches
to capture the effects of dynamic loads. However, a simple torsional model was also developed to
perform some estimates of cyclic variations of torque and speed under various operating conditions.
A summary of this model is also presented.

21

CHAPTER 2
EQUATIONS OF MOTION FOR BLADE VIBRATION

2.1

Formulation of In-Plane EOM

We formulate the equations of motion of in-plane deflection of a beam attached to a hub which
rotates at a constant rate. The equation of motion of the rotating beam is obtained by applying the
extended Hamilton’s principle [72]. The energy formulation for an inextensible beam with large
deflections, rotating about a horizontal axis with flexure is discussed in this section.
The gravitational potential energy of the rotating beam is integrated over the beam elements
with the help of Figure 2.1. The height h(x, t) of an element dx, located at a point x of the
beam, rotated by an angle φ, is given by the angular inclination x(1 − cos φ) plus the inclined
foreshortening s(y, x, t), plus the inclined transverse beam displacement, y(x, t), such that h(x, t) =
∍ x 02
y 04
y
+
)dx and the primes indicate
x(1−cos φ)+ s(x, t) cos φ+ y(x, t) sin φ, where s(y, x, t) =
(
2
8
0
partial derivatives with respect to x. Then the gravitational potential energy is

Vg =

∍ L
0

m(x)gh(x, t)dx =

∍ L
0

m(x)g[x(1 − cos φ) + s(x, t) cos φ + y(x, t) sin φ]dx

(2.1)

where L is the length and m(x) is the mass per unit length.
If the beam is rotating at a fixed angular speed Ω, then φ = Ωt. In a wind turbine, the speed
will naturally be a slowly changing dynamic variable, effected by the aerodynamic forces on the
blades, and the moment on the hub including effects such as generator dynamics, gearbox dynamics,
friction and control. In our resonance analyses, we restrict motion to a fixed angular speed. Any
variation about the mean speed can be dealt with as an extension of this analysis.
The flexural potential energy will include geometric nonlinearity. The nonlinear curvature is
given as k =

y00
.
(1+y02 )3/2

The square of the curvature is expanded into the form k 2 ≈ y002 (1 − 3y02 ),

and inserted into the bending potential energy as
22

g
+

ф
x

.

.

x(1-cos ф)

,t
y(x

.

) P

s(x, t)

h(x, t)

P

Figure 2.1: The arbitrary configuration of the position of a point P on the rotated and deflected
beam. The point P corresponds to location x on the undeflected beam.

23

∍ L
∍ L
1
1
2
Vb =
E I(x)k dx ≈
E I(x)y002 (1 − 3y02 )dx,
(2.2)
0 2
0 2
where E is the Young’s modulus and I is the area moment of inertia of a cross-section of the beam.
The kinetic energy of the beam is formulated from the velocity of each element. The position
of an element is r = (x − s(x, t))er + y(x, t)eφ in terms of the unit vectors in the direction of
Û φ and
the undeformed, rotated beam (er ) and transverse in the φ direction (eφ ). Using eÛ r = φe
Û r , and considering the case where φÛ = Ω, the velocity is v = (−sÛ −Ωy)er +[Ω(x −s)+ yÛ ]eφ .
eÛ φ = −φe
Û 2, where J(x) is the
The contribution to the integral of the rotational energy density, 21 J(x)(y0 + φ)
mass moment of inertia per unit length, can also be added. Then
∍ L
∍ L
1
1
Û 2 dx
T=
m(x)v ¡ vdx +
J(x)(y0 + φ)
0 2
0 2

(2.3)

is the total kinetic energy.
The extended Hamilton’s principle is now applied, such that
∍ t
2
(δT − δV + δW)dt = 0

(2.4)

t1

under the constraint that the varied path coincides with the true path at t = t1 and t = t2 . Constructing
δT, δV, and δW, and integrating by parts to obtain a common δy term in the integrands, yields the
partial differential equation and boundary conditions.
∍ t ∍ L
∍ t
2
2
δVg dt =
m(x)g cos φ(δs)dxdt,
We sketch the treatment of one of these terms,
t
0
t
1
1
∍ x
where δs =
(y0 + y03 /2)δy0 dz, written then as
0

∍ t ∍ L
2
t1

0

m(x)g cos φ

∍ x
0

(y0 + y03 /2)δy0 dz


dxdt.

On the nested integral, we integrate by parts with u = (y0 + y03 /2) and dv = δy0 dz, to produce

∍ t ∍ L
x ∍ x

2

3
3
0
0
0
0
0
m(x)g cos φ (y + y /2)δy  −
(y + y /2) δydz dxdt,
t1

expanded
as,
∍ t ∍ L
2

t1

0

0

0

x 
m(x)g cos φ
 dx
0
0
∍ x
 
∍ L
3
0
0
0
−
m(x)g cos φ
(y + y /2) δydz dx dt.



(y0 + y03 /2)δy
0

0

24

The first term endpoint evaluations lead to an x-dependent term in the first integral, and an x = 0
evaluated term that can be separated out while the integration variable is renamed z. Integrating by
parts on the second integral, the nested integral is u, and the rest is dv. Focusing in the dx integral
inside the dt integral, we have
∍ L

m(x)g cos φ(y0 + y03 /2)δydx

∍ L

m(z)g cos φdz
 ∍ x
 L
x

3
0
0
0
−
m(z)g cos φdz
(y + y /2) δydz 
0
0
0
0

∍ L ∍ x
m(z)g cos φdz y0 + y03 /2 δydx.
+
0∍

0

−



0

(y0 + y03 /2)δy



0

Evaluating the third term, the contribution is zero when x = 0. Therefore, then renaming one
of the integration variables from z to x,
∍ L

∍ L



m(z)g cos φdz [(y0 + y03 /2)δy]
0∍
0
 ∍ L

L
3
0
0
0
−
m(z)g cos φdz
(y + y /2) δydx
0
0

∍ L ∍ x
+
m(z)g cos φdz (y0 + y03 /2)0 δydx.
0

m(x)g cos φ(y0 + y03 /2)δydx

−

0

The third term has the form
∍ L ∍ L
−

0

0



m(z)g cos φdz (y0 + y03 /2)0 δydx

and therefore can be combined with the
term. As such,
∍fourth
 our integral is
∍ L
L
m(x)g cos φ(y0 + y03 /2)δydx −
m(z)g cos φdz [(y0 + y03 /2)δy]
0
0

∍ L ∍ L
−
m(z)g cos φdz (y0 + y03 /2)0 δydx.
0

x

In summary, the term

∍ t ∍ L
2
t1

0

m(x)g cos φδsdxdt has been reworked such that the integral on

the x variable, expressed as the above integral, contains either integral terms with a δy term or
boundary terms (in this case at x = 0).

25

A similar analysis on all of the integrals in δH = 0 put the grand integral in the form such that
inside the t integral, there are x integrals from 0 to L with a common δy, and then boundary terms
evaluated at either x = 0 or x = L. Since δy is arbitrary, its coefficient (an integrand) must be zero,
leading to a partial differential equation of motion. Similarly, at the x = L boundary, δy and δy0
are considered arbitrary, and so their coefficients are both set to zero, producing natural boundary
conditions at x = L.
We assume the x = 0 boundary to be clamped, and thus impose y(0, t) = y0(0, t) = 0 as
geometric boundary conditions.
To this end, we obtain an integral-partial differential equation (IPDE) of motion and boundary
conditions. For the case when φ = Ωt,
−m(a(Üs, s, yÛ, x, t) + g cos Ωt)(y0 + y03 /2) +

∍ L

m(a(Üs, s, yÛ, z, t) + g cos Ωt)dz(y0 + y03 /2)0

x
0
0
00
+mb( yÜ, y, sÛ, x, t) + (J(x) yÜ ) − (E I y − 3E I y00 y02 )00 − (3E I y002 y0)0 +

f (y, yÛ, x, t) = mg sin Ωt,
(2.5)

where the dots are partial derivatives with respect to time, with boundary conditions y(0, t) =
y0(0, t) = 0 at x = 0 and
−J(L) yÜ0 + (E I y00 − 3E I y00 y02 + 3E I y002 y0) = 0

(2.6)

E I y00 − 3E I y00 y02 = 0
at x = L, where a(Üs, s, yÛ, x, t) = sÜ +2Ω yÛ +Ω2 (x −s), b( yÜ, y, sÛ, x, t) = − yÜ +Ω2 y +2ΩsÛ, and f (y, yÛ, x, t)
represents distributed aeroelastic loads.
Interpreting these equations of motion, the IPDE has nonlinear inertial terms due to the geometric nonlinear term m(a(Üs, s, yÛ, x, t) + g cos Ωt)(y0 + y03 /2), along with linear inertial terms,
including rectilinear, centripetal, and Coriolis effects. The IPDE has parametric excitation, with
the cos Ωt terms, as well as direct excitation mg sin Ωt at the same frequency. These come from
the gravitational potential energy, with parametric excitation involving large-deflection nonlinearity. The aerodynamic force f has cyclic components and will contribute to direct and parametric
excitation terms as well. The E I terms include nonlinear curvature terms from large deflection.
26

The aeroelastic f (y, yÛ, x, t) terms have further dependence on s and sÛ, and produce additional direct
or parametric excitation terms.

2.1.1

Modal Reduction of the EOM

To make the integral-partial differential equation amenable to the first level of analysis, we perform
a modal reduction to project the IPDE onto a finite number of ordinary differential equations
(ODEs). We use cantilevered beam modes as assumed modes and apply a Galerkin projection. To
this end, we approximate the transverse deflection as

y(x, t) ≈

N
Õ

qi (t)ψi (x),

(2.7)

i=1

where N is the number of retained modes, ψi (x) are the assumed modal functions, and qi (t) are the
assumed modal coordinates.
We substitute this expression into equation (2.5), multiply by ψ j (x), and integrate over the
length of the blade to obtain the j th second-order ODE. To simplify the analysis, we take N = 1
and neglect the effect of J(x). Neglecting the contributions of the rotational inertia is a common
approximation for the analysis of thin beams [72]. Also, we expect the dominant failure mode to be
in the first resonant mode of the beam. Including higher order modes will reveal additional details
but our aim is to start simple and analyze the representative system for dynamic instabilities.
Employing the aforementioned simplifications to the model, the differential equation takes the
form

Ü 2 + c qq
Û + d qÛ2 q + (e1 + e2 cos Ωt)q + ( f1 + f2 cos Ωt)q3
qÜ + bqq
∍ L
∍ L
= g sin Ωt
ψ(x) dx −
ψ(x) f (y, yÛ, s, sÛ, x, t) dx
0

0

where b, c, d, e1, e2, f1, f2 are constant coefficients and can be evaluated according to


∍ L
∍ L∍ z
∍ x
2
2
00
0
0
0
b=
ψ(x) ψ (x)
ψ (u) du dz − ψ (x)
ψ (v) dv dx
0

c=

∍ L
0

x

ψ(x)



0

−2Ωψ(x)ψ 0(x) + 2ψ 00(x)Ω

0

∍ L



ψ(z) dz dx

x

27

(2.8)

d=

∍ L
0

ψ(x)



ψ 0(x)

∍ x
0

2
ψ 0 (v) dv

+ ψ 00(x)

∍ L∍ L
x

2
ψ 0 (u) du dz


dx

z

!
2 − x2
L
)ψ(x) − (E Iψ 00(x))00 dx
e1 =
ψ(x) Ω2 xψ 0(x) + Ω2 ψ 00(x)(
2
0
∍ L

e2 =

∍ L
0



ψ(x) ψ 0(x)g + gψ 00(x)(L − x) dx

∍ x 02
3
3
ψ 0 (x) 2
Ω2 2
ψ 0 (x) 0
ψ (v)
0
2
2
f1 =
ψ(x) (−
)Ω x + ψ (x)Ω
dv +
(L − x )(
)
2
2
2
2
0
0
"∍ ∍
#

L
z ψ 02
2
2
00
2
00
0
0
00
0
0
−
du dz ψ (x)Ω + E Iψ (x)(3ψ (x)) − (3E Iψ (x)ψ (x)) dx
2
x
0
∍ L



!
3
(L − x)g 03
ψ 0 (x)
g+
(ψ (x))0 dx
f2 =
ψ(x) −
2
2
0
∍ L

These coefficients are dependent on the assumed modal function ψ(x) and the distributed parameters
m(x) and E I(x). The form of f (y, yÛ, s, sÛ, x, t) influences the last integral in equation (2.8), and the
q terms born from it.
For example, if we approximate the system as a uniform beam, such that m(x) and E I(x) are
constants, and use the first modal function of a uniform cantilevered Euler-Bernoulli beam [77],
given by
ψ = [(cosh(λx/L) − cos(λx/L)) − δ(sinh(λx/L) + sin(λx/L)],

(2.9)

where Ν = 1.87510407 and δ = 0.7341, then we obtain parameter values that produce

Ü 2 − 0.918qq
Û + 1.37qÛ2 q + (−2.66E I L −2 + 0.363L 2 Ω2 + 5.12L cos Ωt)q
qÜ − 0.94qq
+(5.352E I L −3 + 13.05E I L −4 − 0.226Ω2 − 8.534L −1 cos Ωt)q3
∍ L
= 1.65L sin Ωt −
ψ(x) f (y, yÛ, s, sÛ, x, t) dx.

(2.10)

0

Including a model for aeroelastic loading via f will contribute some terms (with dependence
on q, L, etc.) already present in equation (2.10) and also some new terms. Wind shear and
28

tower passing will introduce cyclic terms in f which will contribute to the direct forcing term, and
possibly other parametrically excited terms. Aeroelastic forces will also be significant on models
of out-of-plane deflection. It is apparent, from inspection, that this resulting ODE in q has linear
and cubic stiffness effects, with parametric excitation on both the linear and cubic terms, along with
direct excitation.
In addition we include a damping term in our model to account for operational damping in the
system due to resistance in blade motion in the plane of motion. Also the cross coupled terms in
Ü q,
Û q give rise to lot more interplay between the terms in a purely analytical case study. Since the
q,
mode shapes are fixed, the coefficients are only dependent on the geometry of the blade (beam), its
material properties and the rotational frequency.
The aerodynamic force on the beam can be modeled based on the lift and drag force at each
section along the length of the blade. We use the framework of the blade element momentum theory
(discussed in [43]) which is widely used in preliminary calculations in the industry and academics
(detailed discussion in Chapter 5).
Our intent with this formulation is to analyze this system, rebuild the approximate y(x, t) under
resonant cases, and determine the resulting loads applied to the low speed shaft and hub via y(x, t)
and its derivatives. In reality, blade motion is a combination of in-plane motion, out-of-plane
motion and twisting of the blade structure. In the previous section, we developed the EOM for
the in-plane motion. Blades are more pliable in the flap direction due to lesser rigidity in the
out-of-plane mode. This is due to the fact that the wind forces on the blade that cause motion act
predominantly in the out of place direction. The lift force on the airfoil causes the rotation of the
blade and the drag is perpendicular to it. The component of these forces can then be split into
component in the in-plane and out-of-plane directions. Most commercial design codes (including
those discussed in section 1.3.3) ignore the torsion component of the blade motions due to the fact
that its contribution to overall loading on the blade is minimal when compared to in-plane and
out-of-plane loads. As such, a comprehensive method would include all the degrees of freedom
in the modeling approach, however, we seek to isolate one mode of motion to study its effect on

29

the blade motion and consequently, their contribution to drivetrain loads. In essence, the effect of
coupled modes for the 2-dof system would not be included in the current analysis.

2.2

Formulation of Out-of-Plane EOM

We formulate the equations of motion of out-of-plane deflection of a beam attached to a hub
which rotates at a constant rate similar to the process laid out in the previous section. In this section,
we employ the Lagrange approach to arrive at the ordinary differential equations.
Meirovitch [72], compares the Galerkin method, which is the projection of the PDEs onto the
assumed modes, to the Rayleigh-Ritz Method, which is the application of the assumed modes to the
energies (for the linear case). We formulate the out-of-plane equation of motion using the RayleighRitz approach. It has been verified that the two approaches produced the same fundamental equation
of motion that govern the system. Due to the process of integration and differentiation of terms,
the expressions defining coefficients in the ODE for the modal coordinate q(t), appear differently
in both methods, but numerical computations of the coefficients reveal that they are essentially
equivalent.
Following the approach laid out in the previous section, the gravitational potential energy of
the rotating beam is integrated over the beam elements with the help of Figure 2.2. The height
h(x, t) of an element dx, located at a point x of the beam, rotated by an angle φ, is given by the
angular inclination x(1−cos φ) plus the inclined foreshortening s(z, x, t) due to bending in the out of
plane direction, such that h(x,
! t) = x(1 − cos φ) + s(z, x, t) cos φ, where the nonlinear foreshortening
∍ x 02
4
0
z
z
+
dx is evaluated from geometry (see figure 2.3). The primes indicate
s(z, x, t) =
2
8
0
partial derivatives with respect to x. z(x, t) is the out-of-plane deflection at each element along the
length of the blade at a given time. Then the gravitational potential energy is
Vg =

∍ L
0

m(x)gh(x, t)dx
(2.11)

=

∍ L
0

m(x)g[x(1 − cos φ) + s(z, x, t) cos φ]dx

30

Figure 2.2: The arbitrary configuration of the position of a point P on a beam rotating in-plane
and deflecting in the out-of-plane direction. The point P corresponds to location x on the
undeflected beam.
where L is the length and m(x) is the mass per unit length.
The flexural potential energy will include geometric nonlinearity, following [31]. The same
expression was also used in the formulation of the in-plane equation of motion in section 2.1 to
formulate the bending potential energy as
∍ L
∍ L
1
1
2
E I(x)k dx ≈
E I(x)z002 (1 − 3z02 )dx,
Vb =
0 2
0 2

(2.12)

where E is the Young’s modulus and I is the area moment of inertia of a cross-section of the beam.
The kinetic energy of the beam is formulated from the velocity of each element. The position
of an element is r = (x − s(z, x, t))er + z(x, t)ez in terms of the unit vectors in the direction of the
Û φ and eÛ z = 0,
undeformed, rotated beam (er ) and transverse in the φ direction (eφ ). Using eÛ r = φe
and considering the case where φÛ = Ω, the velocity is v = −sÛer +Ω(x − s)eφ + yÛ ez . The contribution
Û 2, where J(x) is the mass moment of inertia
to the integral of the rotational energy density, 12 J(x)(φ)
31

Figure 2.3: The nonlinear foreshortening due to out-of-plane blade motion at a point along the
blade at a distance x
per unit length, can also be added. Then
∍ L
∍ L
1
1
m(x)v ¡ vdx +
J(x)φÛ2 dx
T=
0 2
0 2

(2.13)

is the total kinetic energy.
We approximate the out-of-plane deflection as

z(x, t) ≈

N
Õ

qi (t)ψi (x),

(2.14)

i=1

where N is the number of retained modes, ψi (x) are the assumed modal functions, and qi (t) are
the assumed modal coordinates. We substitute this expression for z into the kinetic and potential
energy expressions, which are now T(qi, qÛi ) and V(qi ).
The Lagrangian function L is defined as L = T − V, where T is the total kinetic energy and V
is the total potential energy (gravitational and flexure) of the system respectively. The generalized
form of the Lagrangian equation of motion is


d ∂L
∂L
−
= Qi 0
dt ∂ qÛi
∂qi

(2.15)

where Qi 0 accounts for frictional forces and cyclic forcing functions (aerodynamic lift and drag
force on the airfoils).
32

To simplify the analysis, we take N = 1 and neglect the effect of J(x). Neglecting the
contributions of the rotational inertia is a common approximation for the analysis of thin beams [72].
Also, we expect the dominant failure mode to be in the first resonant mode of the beam.
Substituting the terms for kinetic and potential energy in equation (2.15) and constructing it as
an equation in the modal coordinate q, we arrive at the differential equation of motion for the first
out of plane mode for the wind turbine blade as

Ü 2 + d qÛ2 q + (e1 + e2 cos Ωt)q + ( f1 + f2 cos Ωt)q3 = Q
qÜ + bqq

(2.16)

where b, d, e1, e2, f1, f2 are constant coefficients and can be evaluated as
b =
d =
e1 =
e2 =

∍ L
∍0

∍ x
m(x)

L

∍0

m(x)

0

0

∍ L
m(x)Ωx

0

∍
0

0

∍ x
m(x)g

0

2
ψ0

∍ x

∍ L

f1 = −

x

2
ψ0

0
L m(x)Ω2

2

2
du

dx
2

du

2
ψ0

dx


du dx +

2



∍ x

2
ψ0

∍ L
0

2

E Iψ 00 dx

ψ 0 du dx

0

(2.17)

2
du

dx


∍ x
∍ L
∍ L
2 2
4
m(x)Ω2 x
0
+
E Iψ 00 ψ 0 dx
ψ du dx − 6
2
0
0
0
∍ x

∍ L
4
1
f2 =
m(x)g
ψ 0 du dx
2 0
o
Equation (2.16) has a very similar structure to the equation (2.8) which governs the in-plane
dynamics. Indeed, all the terms excepts the coupled q qÛ appear in both expressions. The out-ofplane blade ODE has parametric excitation, with the cos Ωt terms, as well as direct excitation Q
which consists of aerodynamic forces. These come from the gravitational potential energy, with
parametric excitation involving large-deflection nonlinearity. The aerodynamic force has cyclic
components and will contribute to direct, and parametric, excitation terms as well. The E I terms
are nonlinear curvature terms from large deflection. The external force terms may have further
33

dependence on s and sÛ, and may produce additional direct or parametric excitation terms. An
aeroelastic formulation may also be accompanied by the appropriate fluid mechanics equation,
coupled through boundary conditions.

34

CHAPTER 3
FORCED NONLINEAR MATHIEU EQUATION

Modeling the in-plane blade and out-of-plane motion in its full depth, i.e. including all the nonlinearities and gravitational loading, leads to a nonlinear parametrically excited equation. Our
goal is to simplify the system of equations and look for analytical descriptions of representative
phenomena. For this we apply reduced-order modeling to obtain analytical expressions of response
characteristics as functions of parameters. We make assumed modal reductions using the cantilevered beam modes to obtain a single mode model of blade motion in the lead-lag direction. The
resulting second order differential equation has elements of a parametrically and directly excited
Mathieu/Duffing equation. This begs us to study it. The simplified equation thus becomes the
focus of our analysis, which is performed in search of possible phenomenon related to wind turbine
systems. The simplifications are relaxed in work where the analysis would be more numerical.

3.1

Perturbation Analysis: First Order

Embedded in the single mode ODE’s (2.8) and (2.16) are terms of a Mathieu-Duffing equation:
qÜ +  µ̂qÛ + (ω2 +  γ cos Ωt)q + αq3 = F sin Ωt,

(3.1)

Since this is a fundamental equation of dynamic systems, equation (3.1) takes our attention and
becomes the focus of an initial study. In this chapter we only consider the influence of sinusoidally
varying direct excitation of the system. For a wind turbine, both sinusoidal and constant direct
excitation forces will occur. This will be treated in detail in Chapter 4.
In equation (3.1), ω plays the role of the undamped natural frequency based on the mean
stiffness-to-mass ratio,  Îł is the magnitude of the parametric excitation corresponding to the cyclic
variation in stiffness, F is the magnitude of direct cyclic forcing due to transverse gravitational
and aerodynamic loading, Ω is the fundamental frequency of cyclic loadings, based on the rotation
rate of the rotor, and  µ̂ represents linear viscous damping. Parameter α is the coefficient of cubic
35

nonlinearity expected due to large deformation effects, and  is a small bookkeeping parameter.
We seek approximate solutions to equation (3.1) by using the method of multiple scales (MMS)
[74,75]. The analysis reveals the existence of various sub-harmonic and super-harmonic resonances.
We unfold these resonance cases and identify the critical ones for wind turbines during normal
operation.
When the forcing in equation (3.1) is of order , the analysis will indicate a primary resonance.
But, this resonance is the same as in the Duffing equation, and is not treated here. The forcing in the
above expression is of order one — also known as hard forcing. This will help us unfold secondary
resonances. Using MMS, we incorporate fast and slow time scales (T0, T1 ) and also variations in
amplitude. This allows for a dominant solution q0 and a variation of that solution q1 , i.e. such that

q = q0 (T0, T1 ) +  q1 (T0, T1 ) + ...

(3.2)

d = D +  D and D = ∂ .
where Ti =  i T0 . Then dt
i
0
1
∂Ti

We substitute this formulation into equation (3.1) and then simplify and extract coefficients of
 0,  1 . The expression for co-efficient of  0 is
D0 2 q0 + ω2 q0 = F sin ΩT0 .

(3.3)

q0 = AeiωT0 − iΛeiΩT0 + c.c.

(3.4)

The solution for this is

F
1
, and A = aei β .
2
2(ω2 − Ω2 )
Since equation (3.3) is a partial differential equation, the coefficient A, and hence a and β are

where Λ =

functions of T1 .
The expression for the co-efficient of  1 is
D0 2 q1 + ω2 q1 = − µ̂D0 q0 − 2D0 D1 q0 − γq0 cos ΩT0 − αq0 3 .
36

(3.5)

Substituting the solution for q0 from equation (3.4), we expand the terms on the right hand side.
We need to eliminate coefficients of eiωT0 that constitute the secular terms and would make the
solutions unbounded. The solvability condition, set by equating the coefficients of eiωT0 terms to
zero, leads to slow flow equations that describe the slowly varying variables a and β of the leading
order solution given by equation (3.4).

3.1.1

Non-Resonant Case

If there is no specific relation between Ω and the natural frequency (ω) of the system, then an
analysis of the solvability condition leads to the conclusion that a → 0 at steady-state solution.
Hence there is no effect of the nonlinear terms in the non-resonant case.

3.1.2
3.1.2.1

Superharmonic Resonances
2Ω ≈ ω

If 2Ω ≈ ω, then the α cos Ωt term and the nonlinearity from equation (3.5) contribute to the
secular terms. We detune the frequency of excitation such that 2Ω = ω +  σ. In this case the
solvability condition can be written as
−2A0iω − µ̂Aiω − α(3A2 Ā + 6Λ2 A) +

iγΛ
exp iσT1 = 0
2

(3.6)

where 0 denotes differentiation with respect to T1 . Letting A = 21 aei β and φ = σT1 − β, the real
and imaginary parts of equation (3.6) lead to
a2
γΛ
aφ0ω − σωa + 3αa(Λ2 + ) +
sin φ = 0,
8
2
µ̂aω γΛ
a0 ω +
+
cos φ = 0.
2
2
At steady-state, a0 = 0, φ0 = 0 and the relationship between the response amplitude and the
detuning parameter σ is

37

0.9

Response Amplitude of the system

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0

0

0.5

1

1.5

Frequency Ratio (Ω/ω)

Figure 3.1: Simulated response of equation (3.1) showing superharmonic resonances at orders 1/2
and 1/3; µ̂ = 0.5,  = 0.1, α = 0.5, F = 0.1, γ = 3. The frequency ratio sweeps up.

σ=

3αΛ2

 2 2

γ Λ
3 2
µ̂2 1/2
+ Îąa Âą
−
.
8
4
4a2 ω2

(3.7)

The peak amplitude would be
ap =

γΛ
.
ω µ̂

(3.8)

The corresponding value of σp is


3αΛ2
Îł2
1+
.
σp =
ω
8ω2 µ̂2
Hence we can conclude the following:
• the peak value is independent of α,
• α, Λ and γ/ µ̂ affect the peak location and as they increase, |σ2p | increases,
• the sign of α determines the sign of σp .

Vertical tangencies in the curve of equation (3.7) can be found and shown to coincide with
stability transitions, such that the unstable branch lies between the vertical tangencies.
38

A sample response of equation (3.1) is numerically simulated. The response curve shown in
Figure 3.1 has the primary frequency and two superharmonics shown at 1/3 and 1/2 the natural
frequency. Superharmonic response at this order is also noticed for a linear system (see Figure 3.2).
This agrees with our calculation of peak amplitude (see equation (3.8)). When nonlinearity Îą = 0,
the Îł term still shows up in the solvability condition (3.6), and this superharmonic resonance peak
persists.
Wind turbine blades are generally designed such that the natural frequency of the rotating blade
in lead-lag (in-plane) motion is below the rotational frequency. This analysis implies potential
existence of superharmonic resonances which would also provide additional critical frequencies
where the response of the blade would be dominant. Such responses would increase loading on the
gearbox and other components and increased bending on the blades. The superharmonic resonance
has also been observed in simulations and experiments on out-of-plane motion [54, 55]. The sign
of the nonlinear term also has a significant effect on the response of the system. Based on the
expression for Îą given in [78], the value can be negative for certain wind turbine blades that would
lead to a softening type behavior (the resonance curve bends to the left). This would be consistent
with the slightly lower frequencies associated with large oscillation of the blades.

3.1.2.2

3Ω ≈ ω

If 3Ω ≈ ω, the cubic nonlinearity contributes to the secular terms. We detune the frequency of
excitation such that 3Ω = ω +  σ. The solvability condition will have additional terms that have
3ΩT0 as an exponential argument in equation (3.5), and it takes the form
−2A0iω − µ̂Aiω − α(3A2 Ā + 6Λ2 A) − iαΛ3 eiσT1 = 0.

(3.9)

Following the analysis done in the previous section for the 2Ω superharmonic resonance, we
substitute a polar form for A and separate the equation into real and imaginary parts. Using
φ = σT1 − β, we arrive at a homogenous set of equations

39

Response Amplitude of the system

15

10

5

0

0

0.5

1

1.5

Frequency Ratio (Ω/ω)

Figure 3.2: Simulated response of the linear case of equation (3.1) showing superharmonic
resonances at orders 1/2 and 1/3; µ̂ = 0.5,  = 0.1, α = 0, F = 0.5, γ = 3
a2
aφ0ω − σωa + 3αa(Λ2 + ) − αΛ3 sin φ = 0
8
µ̂aω
0
+ αΛ3 cos φ = 0
aω+
2
in amplitude a and phase φ.
These expressions are the same as the ones we would get in the synthesis of a standard Duffing
equation with hard excitation [74]. Steady state solutions are satisfied by
! 
 
2
2 2
2 6
 ¾ 2
3αΛ
3Îąa

 a2 = α Λ .
+
σ
−
−
 2

ω
8ω
ω2





This is a quadratic expression in the detuning parameter σ. Solving for σ, we get
! 1/2
3αΛ2 3αa2
α 2 Λ6  µ  2
σ=
+
Âą
−
.
ω
8ω
2
ω2 a2

(3.10)

The peak amplitude and corresponding value of σ are

ap =

2αΛ3
µω

40

(3.11)

and


α 2 Λ4
3αΛ2
1+
σp =
.
ω
2ω2 µ2
Therefore, the peak amplitude and frequency of this resonance are affected by the nonlinearity
via α and by the amplitude of excitation via Λ. The parametric excitation term does not contribute
to this resonance at first order. The value of peak amplitude at this resonance is directly proportional
to Îą (from equation (3.11)). However, numerical simulations for a linear version of our system
show the persistence of this superharmonic peak (see Figure 3.2). Higher order perturbation
analysis (discussed in chapter 4) reveals the existence of this resonant condition in the absence of
nonlinearity.

3.1.3
3.1.3.1

Subharmonic Resonances
Ω ≈ 2ω

If Ω ≈ 2ω then the terms with i(Ω − ω)T0 as the exponent add to the secular terms. The linear
parametric excitation term and the nonlinearity contribute to this resonance. The solvability
condition with these added terms is
−2A0iω − µ̂Aiω − α(3A2 Ā + 6Λ2 A) −

γ Ā
exp iσT1 = 0
2

(3.12)

We substitute the exponential form for A, separate the real and imaginary parts, and let
φ = σT1 − 2β, to get

−

a2
Îła
aωφ0 σωa
−
+ 3αa(Λ2 + ) +
cos φ = 0
2
2
8
4
µ̂aω γa
a0 ω +
+
sin φ = 0
2
4

Seeking steady state, we let φ0 = a0 = 0, which yields
σa ι 3a3
Îła
− (
+ 3Λ2 a) =
cos φ
2
ω 8
4ω
− µ̂a
Îła
=
sin φ
2
4ω
41

We eliminate φ dependence to obtain the frequency response equation as
"
!# 2
σa 3ιa a2
µ̂2 a2
Îł 2 a2
2
=
−
+Λ
+
2
ω
8
4
16ω2

(3.13)

For this equation we have the trivial solution a = 0 and another set of solutions of the form
a2 = p1 ± (p21 − q1 )1/2 , where
64ω2
4σω
− 8Λ2 and q1 =
p1 =
3Îą
9Îą2

!2


2 ω2 − γ 2 
 σ 3αΛ2
4
µ̂


+
 2 − ω
16ω2 




These are similar in form to the standard Duffing subharmonic (presented in section 3.1.3.2), but
the dependence on parameters differs. The nontrivial solutions for a are real only when p1 > 0 and
p1 ≥ q. These conditions imply that solutions will exist if
Λ2

Îł2
σω
2
and µ̂ ≤
<
6Îą
4ω2

(3.14)

These are two distinct conditions on the parameters of our problem. The parametric excitation
term must be sufficiently large, or the damping sufficiently small. Also, the product of the coefficients of the nonlinear cubic term and Λ2 must be sufficiently small or detuning must be sufficiently
large. However, the presence of the cubic nonlinearity is essential for this subharmonic to occur.

3.1.3.2

Ω ≈ 3ω

If Ω ≈ 3ω, the cubic nonlinearity contributes to the secular terms. We detune the frequency of
excitation such that Ω = 3ω +  σ. Then ΩT0 = (3ω +  σ)T0 = 3ωT0 + σT1 .
The solvability condition takes the form
−2A0iω − µ̂Aiω − α(3A2 Ā + 6Λ2 A) − i3α Ā2 ΛeiσT1 = 0

(3.15)

which is the same as that of the standard Duffing equation. We substitute the exponential
form for A, separate the real and imaginary parts, and let φ = σT1 − 3β to transform this into an
42

autonomous system. Seeking steady state, we let φ0 = a0 = 0. Eliminating φ leads to the nonlinear
frequency response equation
! 
 2
2
2 2
2 2
 9 µ̂
9αΛ
9Îą

 a2 = 81α Λ a4
+
σ
−
−
 4
ω
8ω 
16ω2




(3.16)

The solutions of equation (3.16) are either a = 0 or a2 = p ± (p2 − q)1/2 where
"
#
2 9 µ̂2 2
2
8ωσ
64ω
9αΛ
p=
− 6Λ2 and q =
+ (σ −
)2 .
9Îą
ω
81Îą2 4
Since q is always positive, we need p > 0 and p2 ≥ q. This requires that
!
2
2
αΛ
4ωσ
63αΛ
and
Λ2 <
σ−
− 2 µ̂2 ≥ 0
27Îą
ω
8ω
which can be compared to the order-two conditions (3.14). Then
 2
 1/2
4σ
63αΛ2 2σ
=
Âą
− 63
2ω µ̂
µ̂
µ̂2
defines the boundary of the region in the Λσ plane where non-trivial solutions can exist.

3.2

Parameter Dependence Discussion

We have shown the details of the first-order analysis of superharmonic and subharmonic resonance for the system of interest. At first order, the superharmonic resonance of order 1/3 is
the same phenomenon as in the Duffing equation. The nonlinear parameter, Îą, scales the peak
response, while both the nonlinear parameter and the direct excitation level affect the frequency
value of the peak response. The superharmonic response of order 1/2 involves interaction between
the parametric excitation and both the nonlinear parameter and the direct excitation. In fact, if the
nonlinearity is not present, i.e. if Îą = 0, this resonance persists. As such, a linear system excited
both parametrically and directly at the same frequency can exhibit a superharmonic resonance. This
complements linear primary resonance phenomena and linear subharmonic resonances exploited
in parametric amplification [61, 62].
43

Response Amplitude of the system

15

10

Increasing Îł

5

0

0

0.5

1

1.5

Frequency Ratio (Ω/ω)

Figure 3.3: Amplitudes of simulated response of equation (3.1) showing the effect of the
parametric forcing amplitude, with  = 0.1, µ̂ = 0.5, α = 0, F = 0.5. Different curves depict
Îł = 0.5, 1 and 3.
There also exists a primary resonance which is the same at first order as that of the Duffing
equation. The parametric excitation, when at the same frequency as the direct excitation, does not
affect primary resonance amplitude at first order.
Subharmonic resonances at orders two and three exist in the first-order approximations. Ordertwo subharmonics involve interactions with the parametric excitation term, while order-three subharmonic resonance is the same as that of the Duffing equation to first order. The subharmonic
resonance may not be critical to wind turbines, which motivate the study, since they operate below
the natural frequency of the rotating blade.
The variation in system responses for changes in some parameters have been observed in
simulations. Their behavior has been summarized below.
1. Effect of parametric forcing term: The system resonances increase with increased parametric
forcing (Figure 3.3). This amplification occurs at resonances, but not between them. Beyond
a certain value of Îł the system goes unstable at primary resonance. When approaching this

44

60

Response Amplitude of the system

50

40

30

Increasing F

20

10

0

0

0.5

1

1.5

Frequency Ration (Ω/ω)

Figure 3.4: Amplitude of simulated response of equation (3.1) showing the effect of the direct
forcing amplitude, with  = 0.1, µ̂ = 0.5, α = 0, γ = 3. Different curves depict F = 0.5, 1 and 2.
7

Response Amplitude of the system

6

5

4

3

2

Increasing
Nonlinearity
(Îą)

1

0
0

0.5

1

1.5

Frequency Ratio (Ω/ω)

Figure 3.5: Amplitude of simulated sweep-up response of equation (3.1) showing the effect of the
nonlinear term, with  = 0.1, µ̂ = 0.5, γ = 3, F = 0.5. Different curves depict α = 0.1, 0.5 and 1.

45

instability, the primary resonance increases abruptly, which is not captured in the first order
analysis. Higher order perturbation analysis performed in the following chapters details the
onset of instability as we increase parametric forcing term.
2. Effect of direct forcing term: As expected, an increase in the direct forcing term F increases
the overall magnitude of response evenly over the entire spectrum (Figure 3.4).
3. Effect of damping term: Increase in the damping decreases the overall magnitude of response over the entire spectrum, particularly at resonances. The superharmonic and primary
resonance peaks scale as 1/ µ̂.
4. Effect of nonlinearity: An increase in nonlinearity (Îą) causes the response curve to bend
over more significantly. This bending can induce jump instabilities as the frequency slowly
varies. Also, in the presence of strong enough nonlinearity, the response amplitude of the
superharmonic of order two is of a similar order of magnitude as the primary resonance, as
seen in Figure 3.5.
This chapter has focused on the resonances of the forced nonlinear Mathieu equation, and not
stability. The Mathieu equation is well known to have instabilities in the space of stiffness and
parametric excitation parameters. These instabilities are observed when the system is simulated
with certain parameter values at primary resonance. A higher-order analysis is necessary to describe
stability transitions, as well as higher-order effects on resonances. The presence of superharmonic
resonances may be significant for wind turbines, which are designed to operate below the lowest
natural frequency. There are also numerous instabilities at superharmonic frequencies in the
Mathieu parameter space. One of the main motivations to perform high order expansions is to
characterize resonance and instability of the system as a function of system parameters. These can
then be translated to the complete ODE for in-plane and out-of-plane motion i.e., equation (2.8) and
(2.16) to establish limits for material, mass and stiffness of blades in order to not have undeniable
response amplitudes at their operating ranges.

46

CHAPTER 4
SECOND-ORDER PERTURBATION ANALYSIS

In this chapter, we seek to explain the existence of multiple superharmonics for a linear Mathieu
system with a second-order perturbation expansion. Higher-order perturbation techniques have
been used extensively across various sciences to get more accurate expressions of the solutions
to problems, in order to derive an understanding of the underlying physics [79, 80]. Dynamical
systems have been described by higher-order perturbation expansions [81, 82] as well as other
scaling techniques [83]. We extend the second-order formulations to systematically construct
asymptotic solutions of the linear forced Mathieu equation at primary resonance.
The regular Mathieu equation is well known to have instability regions in the parametric forcing
amplitude-frequency parameter space, such that the fixed point at the origin can be stable or
unstable depending on system parameters. The stability wedges of the linear Mathieu equation can
be studied by applying Floquet theory [71] with harmonic balance solutions of periodic transitional
motions [75], and also by a higher-order perturbation expansion [74].
For lightly damped cases, these regions of instability originate near frequency ratios Ω/ω =
2/N, N being an integer (see for example reference [75]). The dependence of the instability region
on the strength of damping has been further explored in detail in [84]. There have been numerous
studies on systems with parametric excitation with nonlinear variations of the Mathieu equation,
for example with added Van der Pol, Rayleigh, and Duffing nonlinear terms, revealing interesting
dynamics (e.g. [64, 85]).
The introduction of direct forcing to the Mathieu equation turns it into a system that does not
directly align it with Floquet theory. However, for the linear Mathieu equation, the solution can be
expressed as a sum of homogenous and particular solutions, given as q = qh +qp , where, qh satisfies
qÜh +2 µqÛh +(ω2 + γ cos Ωt)qh = 0 and qp satisfies qÜp +2 µqÛp +(ω2 + γ cos Ωt)qp = F0 +F sin Ωt.
The qh term, from the Floquet theory, has the various instabilities that originate due to parametric
excitation. The qp part of the solution, approximated by perturbation techniques, can have resonance
47

conditions as well as instabilities. This implies that the instability properties known about the linear
Mathieu equation apply to the system under study in this work. Some of the resonance conditions
obtained through solving the particular solution coincide with the instability property of the linear
Mathieu as will be explained further.
In this chapter, we analyze equation (3.1) under hard forcing with two orders of multiple-scales
expansion in order to capture superharmonic resonance at one-third for the linear system with hard
forcing. We extend our second-order perturbation analysis for a weakly forced system in order
to describe a parametric amplification at primary resonance, and to capture information regarding
stability transition curves. Additional sections also provide commentary on the influence of direct
forcing terms and nonlinearity towards to response of the system for various parameter ranges.

4.1

Second-Order Perturbation Analysis – Hard Forcing

Figure 3.3 shows the simulated response of the system in equation (3.1) for a system with O(1)
excitation. Based on first-order perturbation analysis, the existence of superharmonic resonance at
order one-half is expected for a system with direct and parametric excitation at the same frequency.
Numerical simulations confirm this, but also show a peak at a third of the natural frequency. We
thus proceed to conduct the second-order perturbation analysis of equation (3.1) with the specific
goal of describing the superharmonic resonance of order 1/3.
Equation (3.1) is written with hard excitation and sinusoidally varying direct force as shown
below:
qÜ + 2 µqÛ + (ω2 +  γ cos Ωt)q = F sin(Ωt + θ).

(4.1)

Employing MMS, we incorporate three time scales (T0, T1, T2 ), and allow for a dominant solution
q0 and small variations q1 and q2 , such that
q = q0 (T0, T1, T2 ) +  q1 (T0, T1, T2 ) +  2 q2 (T0, T1, T2 ) + ...

48

(4.2)

where Ti =  i t. Then

d
∂
= D0 +  D1 +  2 D2 , where Di =
.
dt
∂Ti

We substitute this into equation (3.1) and then simplify and balance the coefficients of  0,  1
and  2 :
O(1) :

D0 2 q0 + ω2 q0 = F sin ΩT0

O() :

D0 2 q1 + ω2 q1 = −2µD0 q0 − 2D0 D1 q0 − γq0 cos ΩT0

O( 2 )

D0 2 q2

:

+ ω2 q2

2

= −2D0 D1 q1 − (D1 + 2D0 D2 )q0

(4.3)

−2µ(D0 q1 + D1 q0 ) − γq1 cos ΩT0
Solving the O(1) equation, we arrive at the solution of q0 as

q0 = AeiωT0 − iΛeiΩT0 + c.c.

(4.4)

Feiθ
1
, A = aei β , and c.c. means complex conjugate, and where ω differs from
2
2(ω2 − Ω2 )
Ω (no primary resonance).
where Λ =

The coefficient A, and hence a and β, are functions of T1 and T2 . Substituting this into the O()
of equation (4.3), we arrive at the expression

D0 2 q1 + ω2 q1 = −2µ(AiωeiωT0 + ΛΩeiΩT0 ) − 2D1 AiωeiωT0
(4.5)
Îł
− (Aei(ω+Ω)T0 + Āei(Ω−ω)T0 − iΛe2iΩT0 ) + c.c.
2
Here at first-order, Ω ≈ ω/2 and Ω ≈ 2ω lead to secular terms. This was studied in the
first-order expansions in [86]. But, Ω ≈ ω/3 is not yet revealed as a combination that would lead to
secular terms. As we are seeking that specific resonance condition, we proceed to find the solution
of q1 of the nonresonant case of equation (4.5).

49

4.1.1

Superharmonic Resonances at ω/3

The solvability condition at O() for the non-resonant condition is −2µAiω − 2iωD1 A = 0. Hence,
the particular solution for q1 is
ÎłA
γ Ā
2µΛΩ iΩT
e 0+
ei(Ω+ω)T0 +
ei(Ω−ω)T0
(ω2 − Ω2 )
2(Ω2 + 2ωΩ)
2(Ω2 − 2ωΩ)
iγΛ
e2iΩT0 + c.c.
+
2
2
2(ω − 4Ω )

q1 = −

(4.6)

We substitute the solutions for q0 and q1 into the O( 2 ) expression in equation (4.3) and arrive
at


iγ(Ω + ω)D1 A i(Ω+ω)T
iγ(Ω − ω)D1 Ā i(Ω−ω)T
0
0
e
+
e
= −2
2(Ω2 + 2ωΩ) 
2(Ω2 − 2ωΩ)
2iΩ2 µΛ iΩT
i(Ω + ω)γ A i(Ω+ω)T
iωT
iωT
2
0
0
0
−(D1 Ae
+ 2iωD2 Ae
) − 2µ −
e 0+
e
2
2
(ω − Ω )
2(Ω2 + 2ωΩ)


i(Ω − ω)γ Ā i(Ω−ω)T
2µΛΩ iΩT
ΩγΛ
2iΩT
iωT
0
0
0
+
e
−
e
+ D1 Ae
−γ −
e 0
2
2
2
2(Ω − 2ωΩ)
(ω − 4Ω )
(ω2 − Ω2 )

ÎłA
γ Ā
iγΛ
1 iΩT
i(Ω+ω)T
i(Ω−ω)T
2iΩT
0+
0+
0
+
e
e
e
e 0
2
2
2
2
2
2(Ω + 2ωΩ)
2(Ω − 2ωΩ)
2(ω − 4Ω )
D0 2 q2

+ ω2 q2



(4.7)

+c.c.
We recognize that the terms in the differential equation for q2 involve q0 and q1 differentiated
over different time scales. Only the γq1 cos ΩT0 term, i.e. the last term in equation (4.7), has
interaction between the response and excitation. Expanding that term produces the exponential
term e3iΩT0 . This is secular when 3Ω ≈ ω, giving rise to the 1/3 superharmonic.
While a superharmonic resonance of order 1/3 is apparent, we note that the analysis does not
reveal a subharmonic resonance of Ω ≈ 3ω in this linear forced Mathieu equation.
The solvability conditions up to O( 2 ) for 3Ω ≈ ω are

50

O() :

−2µωAi − 2D1 Aiω = 0

O( 2 ) : −D1 2 A − 2iωD2 A − 2µD1 A −

Îł2 A
Îł2 A
−
4(Ω2 + 2ωΩ) 4(Ω2 − 2ωΩ)

(4.8)

iγ 2 Λ
eiσT1 = 0
2
2
4(ω − 4Ω )
where σ is the detuning parameter defined by 3Ω = ω +  σ.
−

From this it is clear that the forced linear Mathieu equation has a one-third superharmonic,
which required a second-order analysis to be captured. Furthermore, as evidenced from figure
3.3 the response is  order lower than the peak at one-half. Analytical expressions for amplitude
derived in later sections also corroborate this. In order to get an expression for A we need to solve
the solvability conditions at O() and O( 2 ) together. We do this by recombining the time scales
in a process called reconstitution [74]. Eliminating D1 A and D12 A from the O( 2 ) equation, using
ω + σ
, we can show that the resulting solvability
our expression for Λ, and substituting Ω =
3
conditions come from a multiple-scales expansion of
!
2A
2F
dA
9Îł
81iÎł
−2iω
− 2iωµA +  2 µ2 A +
−
eiσt = 0
dt
70ω2
320ω4

(4.9)

We seek a solution in the form A = (X + iY )eiσt , with real X and Y . We enforce this solution
in equation (4.9), separate imaginary and real parts and cancel the common exponential term to
obtain

2ω

9Îł 2Y
81Îł 2 F
dX
+ 2 σωY − 2 µωX −  2 µ2Y −  2
+ 2
cos θ = 0
dt
70ω2
320ω4

(4.10)

dY
9Îł 2 X
+ 2 σωX − 2 µωY +  2 µ2 X +  2
=0
dt
70ω2

(4.11)

2ω

We have to solve for X and Y to get relations between γ, F, µ, , σ and ω. In a standard unforced
Mathieu analysis, we would then seek a solution of the form (X, Y ) = (x, y)e kt which, through
k, would describe the stability of the zero solution. However, this is not appropriate with direct

51

forcing [74], as the origin is not a solution for the above set of equations. Instead, we can seek
equilibrium solutions for X and Y equations (4.10) and (4.11).
1
Alternatively, we can write A in polar coordinates, i.e. we let A = aei β in equation (4.9), and
2
separate real and imaginary parts to get
aω βÛ +  2

Âľ2 a
81Îł 2 F
9Îł 2 a
+ 2
sin(σt + θ − β) = 0
+ 2
2
140ω2
320ω4

−ω aÛ − ωµa −  2

81Îł 2 F
cos(σt + θ − β) = 0
320ω4

(4.12)
(4.13)

We make the system of equations (4.12) and (4.13) autonomous by letting σt − β = φ and
Û For steady state solutions we set
correspondingly σ − βÛ = φÛ to get to expressions for aÛ and φ.
aÛ = φÛ = 0, and square and add the two equations to get

aσ +

 2 Âľ2 a
2ω

+

!2
2
2
 9Îł a
140ω3

+ ( Âľa)2

 2

 81Îł 2 F 2
=
320ω5

(4.14)

Equation (4.14) has the solution

a=
where

r

q
p

9Îł 2
 Âľ2
+
σ+
2ω 140ω3

p =

(4.15)
!2

+ Âľ2
(4.16)


81Îł 2 F 2
q =
320ω5
Hence, we have an expression for the amplitude of the system as a function of all its parameters.


Since p and q are positive, a solution exists over the entire parameter space. To find amax we
da2
differentiate a2 from expression (4.15) with respect to σ and set
= 0. This yields
dσ
!
 Âľ2
9Îł 2
2 σmax +
+
a2 = 0
2ω 140ω3
From this we get the value of σmax and the corresponding value of amax as

52

9Îł 2
 Âľ2
+
σmax = −
2ω 140ω3
(4.17)
amax =

81Îł 2 F
320µω5

√
The response amplitude is a = 2 X 2 + Y 2 , and its stability can be ascertained from equations
(4.10) and (4.11). Noting that these equations are in the form xÛ = Ax + F, the equilibrium solution
xe = −A−1 F can be determined from the determinant and trace of A. We find that if µ > 0, the
response at this resonance is stable.

4.1.2

Simulations and Discussions

Figure 4.1 shows a numerical plot of equation (4.14) showing the variation of amplitude with respect
to the detuning parameter. The maximum amplitude in the plot is âmax = 0.23 at σ̂max = −0.07,
which agrees with the values obtained by using equations (4.17). We can approximate the total
amplitude of oscillation by adding the amplitudes of the free and forced oscillation components,
i.e. compute |q| ≈ |a| + |2Λ|, which is approximate because of phase effects, and compare the
response to the numerical simulations shown in figure 3.3. For the parameters shown in figure 3.3
(noting that µ̂ = 2µ, and figure 3.3 involved forcing amplitude at F compared to the current  F),
âmax = 0.46 is calculated from equation (4.17), which corresponds to the rise above the primary
resonance curve. This is consistent with the peak seen at ω = 1/3 for γ = 3 in figure 3.3.
Further numerical simulations were carried out varying the parameters that would make the
system go unstable at primary resonance. However, if a system could operate below Ω ≈ ω
excitation, such as with wind turbines, this instability would not be encountered. Figure 4.2 shows
one such scenario. We notice the existence of multiple peaks. These can be correlated to either
a harmonic frequency or the seeds of instability wedges of the unforced Mathieu equation, which
occur near frequency ratios, Ω/ω = 2/N (N = 1, 2, . . . ) in the frequency response curve [75]. It is
well known that these instability wedges become more slender and weaker as we increase N in a
typical Mathieu plot, and especially weaker in the presence of damping. As we can see from figure
53

0.5

Amplitude of response |a|

0.45

0.4

0.35

0.3

0.25

0.2

0.15

0.1
−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

Detuning Parameter σ

Figure 4.1: Numerical evaluation of equation (4.14) shows the amplitude versus detuning curve.
Graph shown is generated using ¾ = 0.25,  = 0.1, F = 0.5, γ = 3, and θ = 0. amax = 0.2278
occurs at σ = −0.07.
4.2, the response near a ratio of 2/3 is distinct in the sense that the amplitude increases abruptly.
This can be attributed to the presence of an instability wedge, such that the amplitude values are
not steady-state, but are increasing during the time interval of the sweep. The lower-frequency
instability wedges (1/2, 1/3, 2/5, 2/7, etc.) are not penetrated in this simulation. The responses at
the other frequency ratios are due resonance phenomena (1/2, 1/3, etc.). In some cases, there exists
both a resonance curve and an instability wedge. Future work will explore parameters that will
make the system operate in such critical zones.
Figure 4.3 shows the response of the system when excited at frequencies higher than the natural
frequency. We can clearly see the subharmonic instability occurring at twice the natural frequency,
as there is an instability wedge at Ω = 2ω. The figure is plotted with weak parametric excitation to
show the existence of a superharmonic at order 2. The superharmonic resonance conditions exists
even for this case, but the magnitude of response in comparison to primary resonance and order 2
subharmonic is relatively small.

54

Response Amplitude of the system

10
9
8
7
6
5
4
3
2
1
0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Frequency Ratio (Ω/ω)

Figure 4.2: Simulated response of the linear case of equation (4.1) showing multiple
superharmonic resonances ¾ = 0.05,  = 0.1, F = 0.5, γ = 7, and θ = 0. For the isolated
superharmonic resonance at Ω/ω ≈ 1/3, amax = 2.477 by (4.17).

Response Amplitude of the system

5
4.5
4
3.5
3
2.5
2
1.5
1
0.5
0

0

0.5

1

1.5

Frequency Ratio (Ω/ω)

2

2.5

Figure 4.3: Simulated response of (4.1) showing subharmonic resonance at 2;
¾ = 0.05,  = 0.1, F = 0.5, θ = 0, and γ = 0.5. (Low parametric excitation is chosen so that the
system doesn’t show unstable response at primary resonance. For the given set of values, system
is unstable at subharmonic resonance.)

55

4.2

Weak Excitation—Primary Resonance Amplification

In this section we perform second order analysis to uncover resonance conditions and instability wedges under soft excitation (O() forcing). Our primary interest is in describing primary
resonance amplification observed in simulations. However, we will take note of other resonances
that may arise. We focus our attention back to the original equation (4.18) restated as

qÜ + 2 µqÛ + (ω2 +  γ cos Ωt)q =  F sin(Ωt + θ)

(4.18)

We follow the second-order expansion analysis done in the previous section to arrive at expressions for coefficients of  0,  1 and  2 as

O(1) :

D0 2 q0 + ω2 q0 = 0

O() :

D0 2 q1 + ω2 q1 = −2µD0 q0 − 2D0 D1 q0 − γq0 cos ΩT0 + F sin(ΩT0 + θ)

O( 2 ) : D0 2 q2 + ω2 q2 = −2D0 D1 q1 − (D1 2 + 2D0 D2 )q0 − 2µ(D0 q1 + D1 q0 )

(4.19)

−γq1 cos ΩT0
The solution of the O(1) equation is

q0 = A(T1, T2 )eiωT0 + c.c.

(4.20)

1
where A = aei β .
2
The differential equation for q1 is then rewritten by substituting the solution for q0 into the O()
equation as

D0 2 q1 + ω2 q1 = −2iωD1 AeiωT0 − 2µiωAeiωT0
Ā
A
iF
−γ( ei(ω+Ω)T0 + ei(Ω−ω)T0 ) − ei(ΩT0 +θ) + c.c.
2
2
2
56

(4.21)

4.2.1

Weakly Forced Linear Mathieu Equation: Resonance Cases

We have three cases for equation (4.21) which can contribute secular terms:
1. No specific relationship between Ω and ω
2. Ω ≈ ω
3. Ω ≈ 2ω
The treatment of these cases will be laid out in the following sections. In each case, the secondorder expansions will give us the slow time scale variations of the amplitude A of our fast time scale
solution q0 . A second-order expansion will involve solvability conditions obtained by eliminating
secular terms at both O() and O( 2 ).

4.2.2

Case 1: No Resonance at O()

When there is no specific relationship between the forcing frequency Ω and the natural frequency
ω, we remove secular terms by letting −2iωD1 A − 2µiωA = 0, and solve for the particular solution
using the remaining terms in differential equation (4.21), noting that we treat A as constant with
respect to the independent variable T0 . As such,

q1 =

γ Ā
iFeiθ
ÎłA
ei(Ω+ω)T0 +
ei(Ω−ω)T0 +
eiΩT0 + c.c.
2
2
2
2
2(Ω + 2ωΩ)
2(Ω − 2ωΩ)
2(Ω − ω )

(4.22)

We substitute the solutions for q0 and q1 into the equation (4.19) at O( 2 ). After expanding
the resulting q2 differential equation, we seek to identify terms that can contribute to a resonance
condition. Although there were no resonances specified at O(), new resonances can be identified
at O( 2 ). When there is not a specified relationship between Ω and ω at either O() or O( 2 ), i.e. no
resonances are specified, secular terms are removed and the resulting solvability conditions from

57

both O() and O( 2 ) are
O() :

−2iωD1 A − 2µiωA = 0

O( 2 ) : −D1 2 A − 2D2 Aiω − 2µD1 A −

Îł2 A
4(Ω2 + 2ωΩ)

−

Îł2 A
4(Ω2 − 2ωΩ)

(4.23)
=0

Following references [74, 75], these equations can be analyzed for their fixed points and stability.
The removal of the secular terms in equation (4.19) at O( 2 ) leaves a differential equation for
q2 as
D0 2 q2

+ ω2 q2



= −2 i(Ω + ω)D1 Kei(Ω+ω)T0 + i(Ω − ω)D1 Lei(Ω−ω)T0
 


+iΩD1 NeiΩT0 − D1 2 AeiωT0 + 2iωD2 AeiωT0 − 2µ i(Ω + ω)Kei(Ω+ω)T0


Îł
iωT
iΩT
i(Ω−ω)T
0 −
0 + D1 Ae
0 + iΩNe
+i(Ω − ω)Le
Kei(2Ω+ω)T0 + KeiωT0
2

i(2Ω−ω)T
−iωT
i2ΩT
0 + Le
0 + Ne
0 + N + c.c.
+Le
where K =

(4.24)

ÎłA
γ Ā
iFeiθ
,L=
, and N =
, for the nonresonant case.
2(Ω2 + 2ωΩ)
2(Ω2 − 2ωΩ)
2(Ω2 − ω2 )

But there are also specific combinations of Ω and ω appearing at second order that could lead
to resonance conditions in the system, and thereby add terms to the solvability condition at O( 2 ).
These can be identified by examining the remaining terms of the O( 2 ) equation given by equation
(4.24). We arrive at the following candidates:
a. Ω ≈ ω/2
b. Ω ≈ ω
c. Ω ≈ 2ω
The first case was treated with a first-order multiple scales analysis in [86], in which hard
excitation induced secular terms at the O() expansion. Similarly, in this case, the excitation,
this time at O(), induced secular terms one order down at O( 2 ). Hence, similar results may be
58

expected, and so this case is not emphasized in this paper. The latter two cases actually produce
secular terms at O(), and are therefore cannot be treated within Case 1, since it is a nonresonant
case at O().
Also, in this analysis, using weak forcing, we do not capture the linear superharmonic at order
1/3, as with weak forcing this superharmonic is buried “one-order lower". (If we were to continue
our analysis further, i.e., do a third-order MMS expansion, we will encounter that resonance term.)
The increasingly higher-order superharmonics show up at an increasingly lower order of . These
superharmonics will occur at frequency ratios of 1/n, n being an integer. If the Mathieu system
were weakly forced we would have to extend the perturbation analysis to nth order expansions to
capture corresponding effects. If a hard-forced Mathieu equation were considered, the order of
MMS expansion to get a 1/n, 2/n order superharmonic would be (n − 1)th . However, for higher
orders of , the effect of superharmonic resonances become increasingly negligible.

4.2.3

Case 2: Primary Resonance and Parametric Amplification

When Ω ≈ ω, i.e. letting Ω = ω +  σ in equation (4.21), the solvability condition is

O() :

iFeiθ iσT
e 1 =0
−2iωD1 A − 2µiωA −
2

(4.25)

The particular solution for q1 then becomes

q1 =

ÎłA
γ Ā
ei(Ω+ω)T0 +
ei(Ω−ω)T0 + c.c.
2(Ω2 + 2ωΩ)
2(Ω2 − 2ωΩ)

(4.26)

We substitute the solution for q0 and q1 into the equation (4.19) at O( 2 ). After expanding the
right-hand side of the resulting q2 differential equation, we seek to identify secular terms at this
resonance condition.
For the case of primary resonance, we let Ω = ω +  σ in equation (4.19) and eliminate secular
terms to get

59

O( 2 )

Îł2 A
: −D1
2 Aiω − 2µD1 A −
4(Ω2 + 2ωΩ)
2
2
Îł A
γ Ā
−
−
e2iσT1 = 0
2
2
4(Ω − 2ωΩ) 4(Ω − 2ωΩ)
2 A − 2D

(4.27)

Rewriting the O() equation in equation (4.25), we have
D1 A = −µA −

Feiθ eiσT1
4ω

We compute D1 2 A from the above expression as
D1

2A

=

Âľ2 A +

¾Feiθ iσT
iσFeiθ iσT
1
e
−
e 1.
4ω
4ω

We substitute these into the O( 2 ) expression in equation (4.27) to arrive at

iσFeiθ iσT
Feiθ eiσT1 

¾Feiθ iσT
1
1
+
− 2µ − µA −
−2D2 Aiω −
e
e
4ω
4ω
4ω
2
2
2
Îł A
Îł A
γ Ā
2iσT
1 =0
−
−
−
e
4(Ω2 + 2ωΩ) 4(Ω2 − 2ωΩ) 4(Ω2 − 2ωΩ)
µ2 A −

(4.28)

Following the approach in [74], the O() solvability condition from equation (4.27) and the
O( 2 ) solvability condition from equation (4.28) result from a multiple-scales expansion of


dA
iFeiθ iσt
−2iω
+  −2µiωA −
e
+
dt
2
(4.29)


γ 2 Ā 2iσt γ 2 A iσFeiθ iσt µFeiθ iσt
2
2
 Âľ A+
e
+
+
e
+
e
=0
4ω
4ω
4ω2
6ω2
Transforming the system to polar coordinates leads to steady-state equations in a and φ =
σT1 − β. The steady-state equations include sin(φ), cos(φ) and sin(2φ), cos(2φ) terms, which make
it more difficult to solve than other familiar examples.
Instead we follow the solution procedure given in [74]. We seek a solution in the form
A = (X + iY )eiσt , with real X and Y . We enforce this solution in equation (4.29), separate
imaginary and real parts and cancel the common exponential term to obtain
2ω

dX
 2Îł2
+ 2 µωX − (2ωσ +  2 µ2 −
)Y = Îą,
dt
12ω2
60

(4.30)

2ω

5 2 Îł 2
dY
+ (2ωσ +  2 µ2 +
)X + 2 µωY =  δ,
dt
12ω2

(4.31)

F
 2 Fσ
2F Âľ
F
 2 Fσ
2F Âľ
cos θ +
cos θ +
sin θ and  δ = −
sin θ +
sin θ −
cos θ.
2
4ω
4ω
2
4ω
4ω
At steady state XÛ = YÛ = 0. Hence from equations (4.31) and (4.30) we get

where α = −

+ −(2ωσ +  2 µ2 −

2 µωX
(2ωσ

4.2.3.1

+  2 Âľ2

5 2 Îł 2
+
)X +
12ω2

2 µωY

 2Îł2
)Y = Îą,
12ω2

(4.32)

=  δ.

Parametric Amplification

We solve (4.32) for X and Y , and then write the square of the total amplitude r 2 = (X 2 + Y 2 ) as
!

2 2

Îł
2ωµα + δ 2ωσ +  µ2 −
12ω2
!
!
r2 = 
2
2 2

Îł
5
Îł
(2ωµ)2 + 2ωσ +  µ2 −
2ωσ +  µ2 +
2
12ω
12ω2
!
2 2
5
Îł
2ωµδ − α 2ωσ +  µ2 +
12ω2
!
!

2 2
2

Îł
5
Îł
(2ωµ)2 + 2ωσ +  µ2 −
2ωσ +  µ2 +
12ω2
12ω2


+

(4.33)

Equation (4.33) provides information related to the amplitude of the primary resonance as a function
of parameters, including the detuning parameter σ and the phase angle θ between the parametric
1
and direct excitation terms. Note that we had A = aei β = X + iY = reiφ , and as such the response
2
amplitude is a = 2r. Since ι and δ depend on θ, the response amplitude is also dependent on θ.

61

Looking, for example, at the special case of θ = 0, we have
!


2
2  − F µ   F
 Fσ

Îł
2
−
−
(2ωµ)
2ωσ +  µ −
4ω
2
4ω
12ω2
!
!
r02 = 
2
2 2

Îł
5
Îł
(2ωµ)2 + 2ωσ +  µ2 −
2ωσ +  µ2 +
2
12ω
12ω2
!
 

2 2
− F µ
F  Fσ
5
Îł
(2ωµ)
+
−
2ωσ +  µ2 +
4ω
2
4ω
12ω2
!
!

2
2 2

Îł
5
Îł
(2ωµ)2 + 2ωσ +  µ2 −
2ωσ +  µ2 +
12ω2
12ω2


+



(4.34)

In the limit as γ → 0, σ → 0 and for small  (such that we can neglect second order  terms),
F
then in equation (4.34) r0 =
. With r0 = a/2,  µ = ζω and recognizing that the excitation
2ωµ
F
in equation (4.18) is  F, this amounts to a =
, which is consistent with the forced linear
2ζω2
oscillator at resonance.
Equation (4.34) provides a predicted response amplitude as a function of system parameters for
the 1:1 primary resonance of the linear forced Mathieu equation with damping. Figure 4.4 shows
the variation of amplitude as a function of parametric excitation Îł for a given set of parameter
values. The response amplitude predictions agree with numerical simulations of the linear Mathieu
equation given by equation (3.1).
Figure 4.4 shows how the response amplitude (a = 2r) at primary resonance varies with
increasing γ for various values of σ. The response amplitude is nearly uniform for a range
of small values of parametric excitation. In this range of Îł, there is no significant parametric
amplification, and the response levels at Îł = 0 in figure 4.4 are consistent with those of simple
forced linear oscillators. For larger values of Îł, the steady-state amplitude increases hyperbolically
with increasing Îł and becomes unbounded when Îł reaches a critical value, Îłcr . The extension or
stretch in the amplitude is known as parametric amplification.
Figure 4.6 shows the variation in the amplitude of the response versus the detuning parameter,
σ. In this case, the peak occurs at σ = −0.01, and the associated peak amplitude for the given
62

30

20

15

10

5

0

σ = 0.1
σ = 0.25
σ = 0.4

60

Primary Resonance Amplitude

25

Primary Resonance Amplitude

70

σ = −0.25
σ = −0.15
σ = −0.05

50

40

30

20

10

0

0.5

1

1.5

2

2.5

3

0

3.5

0

1

2

3

4

5

6

7

8

9

Parametric Excitation (Îł)

Parametric Resonance (Îł)

Figure 4.4: Variation of maximum amplitude of vibration at primary resonance as a function of
parametric excitation γ for various values of σ. (Left) Excitation frequency is detuned to be
slightly less than natural frequency. (Right) Excitation frequency is detuned to be slightly more
than natural frequency. Parameters are F = 5, Âľ = 0.1 and  = 0.1.
60
1000

50

Primary Resonance Amplitude

900

Primary Resonance Amplitude

Analytical
Numerical

Analytical
Numerical

800

700

600

500

400

300

200

40

30

20

10

100

0

0
0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

Parametric Excitation (Îł)

Parametric Excitation (Îł)

Figure 4.5: (Left) Numerical simulations show that the predicted parametric resonance can have
errors when the amplification becomes large as Îł reaches critical values. (Right) Graphed is
zoomed in to show that numerical and analytical results match for moderate levels of parametric
excitation.
parameter set is a = 2r = 10 for the case of Îł = 0 (no parametric excitation).
Inspection of equation (4.33), including the special case of equation (4.34) for θ = 0, reveals
that denominator D is a fourth-degree polynomial in σ, and the numerator N is quadratic in σ.
Therefore, to find the value of σmax for the maximum amplitude of the resonant response from
d(r 2 )/dσ = (N 0 D − D0 N)/D2 = 0 involves finding roots of a fifth-degree polynomial in the
63

20

Primary Resonance, . from 0 to 4

15

amplitude, a

amplitude, a

15

10

10

5

0
-1

-0.5

0

0.5

5
-0.4

1

<

-0.2

0

0.2

0.4

<

Figure 4.6: (Left) The numerical evaluation of equation (4.34) showing the amplitude versus the
detuning. Graph shown was generated using Âľ = 0.25,  = 0.1, F = 5, Îł ranges from 0 to 4 in 0.5
increments; as Îł increases a increases. (Right) The numerical evaluation of equation (4.34) for
Îł = 3 showing the amplitude versus the detuning at various orders of  expansion. Note that the
curves with O() (dotted line) and O( 2 ) (dashed line) expansions are visibly erroneous.
numerator. So we will not present an analytical expression for σmax or rmax .
However, Figures 4.4 and 4.6 suggest that the peak value may come close to σ = 0, at least for
some parameter sets with small values of Îł (this deviates as Îł increases). An expression for r 2 at
σ = 0 for this approximation of the peak value can be expressed. As such,
!
!
!
2 

2
2 FÂľ
2 2
−
F
Âľ

Îł
F
5
Îł
− ωµF +
−  µ2 −
+
 Âľ2 +
2 4ω
2
2
12ω
12ω2
2 =
!
!
r00


5 Îł 2 2
 Îł2
2
2
2
(2ωµ) +  µ −
Âľ +
12ω2
12ω2

(4.35)

The amplifier gain can be approximated by taking the ratio of amplitudes r00 as predicted by
F
equation (4.35) and the associated value for the simple harmonic oscillator, i.e. a h =
. In
2µω2
the gain expression, the force level F cancels out. For general phases, the expression r 2 /ah can
be evaluated at σ = 0 and plotted versus θ to show how the amplifier gain depends on the phasing
of direct and parametric excitations, as done with previous works on subharmonic parametric
amplification [61, 62].
Since this response is predicted with small-parameter asymptotic theory, it is likely that the
64

predicted amplitude values in the large-response range have some error. This is shown in the plots
of figure 4.5. Next, we will see how the unbounded response relates to stability.

4.2.3.2

Stability of Primary Resonance

Equations (4.31) and (4.30) can be rewritten in the form

2ωxÛ = Ax + F,

(4.36)

where, in the case of θ = 0,


  F  2 Fσ 

 X 
 −
+




x= 
 , F =  2 2 4ω 

 Y 

 FÂľ




−




4ω
2 2 


2 µ2 −  γ 
−2
µω
2ωσ
+



12ω2 
and A = 
2
2

 −2ωσ −  2 µ2 − 5 γ

−2 µω


2
12ω


The equilibrium x0 = A−1 F was found previously, and its amplitude r was studied. Since
equation (4.36) is linear, dynamics of X and Y , and hence the leading-order slowly varying response
amplitude A, about the equilibrium, and hence the stability, are governed by matrix A. In addition,
the homogeneous solution has its own equilibrium at the origin. The homogeneous solution is that
for which F = 0 in equation (4.18), and consequently F = 0 in equation (4.36). Thus, matrix
A simultaneously governs the stability of the homogeneous and particular solutions of equation
(4.18). The stability is determined by the trace (T) and determinant (D) of A, which are
T = −4 µω

(4.37)
2 3 γ 2 σ
5 4 Îł 4
+  4 µ4 −
3ω
144ω4
Since  > 0, µ > 0 and ω > 0 the trace T < 0. Therefore, the stability is governed by
D = 4 2 ω2 (µ2 + σ 2 ) + 4 3 µ2 ωσ +

the determinant D. The instability transition is given by parameters for which D = 0. Close
examination of equation (4.34) shows that the denominator of the r 2 expression is in fact D2 ,
65

Figure 4.7: Stability wedge for the linear forced Mathieu equation with damping at primary
resonance ((ω +  σ)2 = δ = N 2 /4, N = 2 for traditional Mathieu equation). Figure shows
comparison of stability of numerical system response/amplitude and the analytically calculated
stability boundary.
consistent with the denominator of A−1 . Thus, the instability criterion for the forced response
coincides with the approach of r 2 to infinity, and with the instability of xh of the unforced Mathieu
equation. Solving D = 0 for Îł yields the critical parametric excitation as
s
Îłcr =

q
24ω2
12ω2
2
9(2ωσ +  µ2 )2 + 5(2µω)2
(2ωσ +  µ ) ±
5
5

(4.38)

The calculated value of Îłcr is the value of Îł in figure 4.4 for which the amplitude becomes
unbounded for specific values of detuning parameter. The increase in amplitude until that point is a
result of parametric amplification, analytically approximated by equation (4.34). The analytically
predicted values of γcr are plotted versus σ in figure 4.7. Numerical simulations were also checked
for stability in this parameter space, and stable (x’s) and unstable (dots) events are plotted along
with the predicted stability boundary from equation (4.38) with good agreement.

66

4.2.3.3

Summary of Primary Parametric Resonance

We have seen that the linear forced Mathieu equation has a homogeneous and particular solution.
The homogeneous solution has the behavior previously studied traditionally by Floquet theory [71],
perturbation analysis [74], or Floquet-motivated harmonic balance [75]. The particular solution is
a forced response with the potential for a primary resonance, which was analyzed in this section by
multiple scales. A second-order analysis shows that the primary resonance is stretched, or amplified,
by the parametric excitation. Given the parameter set, as the level of parametric excitation Îł grows,
the resonance amplitude is at first hardly affected by small Îł, but then begins to grow quickly for
larger Îł, and approaches infinity at a critical Îłcr . The bounded resonance response is stable. As
the response becomes infinite at Îłcr , the response destabilizes. At the same time, the homogeneous
part of the response also becomes unstable, as the parameters cross into the unstable part of the
Arnold tongue.

4.2.3.4

Cautionary Note on Second-Order Expansions in 

It may be tempting to expand r 2 into powers of  consistent with the level of analysis for the
reconstituted differential equation, for example to powers of  2 in a second-order analysis, as
in [81]. This would lead us to an expression for r 2 from a geometric series expansion:




F 2 γ 2 σ − 24F 2 µ2 ω2 σ − 12F 2 ω2 σ 3
2
+
rapprox
=

 +


5 ¾2 + σ 2 2
16ω2 µ2 + σ 2
192ω

2
 35F 2 γ 4 µ2 − 120F 2 γ 2 µ4 ω2 + 3F 2 γ 4 σ 2 − 264F 2 γ 2 µ2 ω2 σ 2 + 1440F 2 µ4 ω4 σ 2
 

3


2
2
2
4
2
2
4
4
2
4
6
8
2
2
−48F γ ω σ + 1008F µ ω σ + 144F ω σ
9216ω µ + σ
+ O( 3 )
F2

(4.39)

where all terms up to O( 2 ) are kept in the analysis. However, terms of order  3 and higher are
dropped from the expansion of equation (4.34) in . In the limit as  → 0, σ → 0, equation (4.39)
is consistent with the common forced harmonic oscillator.

67

Quantitative predictions for r 2 values computed from equation (4.39) can be quite erroneous.
This can be seen by comparing the values for amplitude obtained by equation (4.34) and (4.39) to
the values obtained numerically in figure 3.3 (recalling that µ̂ = 2µ, and figure 3.3 involved forcing
amplitude at F compare to the current  F). For example, for Âľ = 0.25,  = 0.1, F = 5, Îł = 3
and σ = 0 the amplitude obtained from figure 3.3 at primary resonance (Ω ≈ ω) is |q| = 15.
The analytical calculations from equations (4.34) and (4.39) for amplitude are |q| = 13.95 and
|qapprox | = 13.3 respectively (insight from figure 4.6 suggests some error due to σ = 0 not
σmax ). Furthermore, the transition to unbounded responses as predicted by the denominator of
the simplified equation (4.39) deviates from that plotted in figure 4.7 from equation (4.38) and
the stability data from simulations. Thus the truncated expression including up to O( 2 ) gives an
erroneous result.
This indicates that the process of truncating to the 2nd power of  in effort to achieve -expansion
consistency with the multiple scales expansion is incorrect. Instead, it seems that the reconstituted
equation (4.29) has information of a second-order expansion, and the approximate steady-state
solution (4.34) uses all of this information. The truncated expression (4.39) loses information and
is more quickly divergent from the true solution. Thus, both the results of the amplitude expression
and the stability boundary indicate that we need to consider all powers of  while solving the
reconstituted equation (4.29). We follow the same approach while solving for response amplitude
and stability boundaries at other resonance conditions.

4.2.4

Case 3: Subharmonic Resonance of Order Two

This case was examined in the previous first-order multiple scales analysis [86]. It is possible that
a second-order analysis may improve the asymptotic expressions somewhat, and so we sketch its
analysis. When Ω ≈ 2ω, i.e. Ω = 2ω +  σ, the solvability condition for equation (4.21) is

O() : −2iωD1 A − 2µiωA −

γ Ā iσT
e 1 =0
2

The particular solution of equation (4.21) for q1 then becomes
68

(4.40)

q1 =

iF
ÎłA
ei(Ω+ω)T0 +
eiΩT0 + c.c.
2
2
2
2(Ω + 2ωΩ)
2(Ω − ω )

(4.41)

Applying to equation (4.19) at O( 2 ), and removing secular terms, leaves us with a solvability
condition as
O( 2 ) : −D1 2 A − 2D2 Aiω − 2µD1 A −

Îł2 A
= 0.
4(Ω2 + 2ωΩ)

(4.42)

We combine equations (4.40) and (4.42) and reconstitute to get an ODE





γ Ā iσt
γσ Ā iσt
Îł2 A
dA
2
2
+  −2µiωA −
e
+ Âľ A+
e
−
=0
−2iω
dt
2
4ω
32ω2

(4.43)

We follow the steps laid out in the previous section to arrive at the expressions for trace (T)
and determinant (D) for order 2 subharmonic, from the XÛ and YÛ that govern the stability of the
forced resonant response, and likewise for the homogeneous solution. The trace and determinant
are given by
T = −4ωµ
3 2 Îł 2
D = (2ωµ)2 + 2 σω +  2 µ2 −
32ω2

!2

 γ  2 σγ
−
−
2
4ω

!2

(4.44)

The expression for D is quadratic in Îł 2 which when equated to zero gives the stability boundaries
of the order-2 subharmonic under hard excitation. These agree with detailed analyses in the existing
literature [74,75] and we do not present the details of this analysis. The expressions for the resonance
under weak excitation also appear at an order lower than the subharmonic found using first-order
expansions with hard excitation in [86].

4.3

Effect of a Direct Forcing Term F0

The previous section we have considered a sinusoidal forcing term and parametric excitation at
the same frequency. The motivation for the study of this equation arises from the simplifications
of the wind turbine blade equation of motion formulation. The mean wind loading on a blade
69

represents a constant direct forcing term F0 whose magnitude with respect to the harmonic forcing
is dependent on the direction of motion of the blade (in plane versus out of plane). In this section,
we look at the fundamental dynamics of equation (3.1) with the constant force term. As such, the
Mathieu equation with constant loading is
qÜ + 2 µqÛ + (ω2 +  γ cos Ωt)q = F0 .

(4.45)

We follow the second-order procedure established in sections 4.1 and 4.2 to arrive at expressions
for amplitude of resonance at primary and other harmonic resonance conditions as a function of
system parameters. The multiple scales expansions of equation (4.45) are
O(1) :

D0 2 q0 + ω2 q0 = F0

O() :

D0 2 q1 + ω2 q1 = −2µD0 q0 − 2D0 D1 q0 − γq0 cos ΩT0

O( 2 ) : D0 2 q2 + ω2 q2 = −2D0 D1 q1 − (D1 2 + 2D0 D2 )q0

(4.46)

−2µ(D0 q1 + D1 q0 ) − γq1 cos ΩT0
The solution of the O(1) equation for q0 is

q0 = AeiωT0 + κ + c.c.

(4.47)

F0
.
2ω2
Substituting this into the O() of equation (4.46), we arrive at the expression

where Îş =

D0 2 q1 + ω2 q1 = −2D1 AiωeiωT0 − 2µ(AiωeiωT0 )
(4.48)
Îł
− (Aei(ω+Ω)T0 + Āei(Ω−ω)T0 − κeiΩT0 ) + c.c.
2
We can immediately identify two values of Ω which contribute secular terms (Ω ≈ ω and
Ω ≈ 2ω). For the case when there is not a specific relation between the frequencies of parametric
excitation and natural frequency, the Ω equation assumes the form

70



D0 2 q2 + ω2 q2 = −2 i(Ω + ω)D1 Kei(Ω+ω)T0 + i(Ω − ω)D1 Lei(Ω−ω)T0


 
iΩT
iωT
iωT
2
+iΩD1 Ne 0 − D1 Ae 0 + 2iωD2 Ae 0 − 2µ i(Ω + ω)Kei(Ω+ω)T0


Îł
i(Ω−ω)T
iΩT
iωT
0 + iΩNe
0 + D1 Ae
0 −
Kei(2Ω+ω)T0 + KeiωT0
+i(Ω − ω)Le
2

i(2Ω−ω)T
−iωT
i2ΩT
0 + Le
0 + Ne
0 + N + c.c.
+Le

(4.49)

ÎłA
γ Ā
ÎłF0
,L=
, and N =
, for the nonresonant case.
2(Ω2 + 2ωΩ)
2(Ω2 − 2ωΩ)
4(Ω2 − ω2 )
There is no contribution from F0 to the secular terms at O( 2 ) in the non-resonant case.

where K =

The resonant cases and the associated secular terms for two orders of multiple scale expansions
are given below (also referred to as solvability conditions).
Case 1: Ω = ω +  σ

O() :

−2iωD1 A − 2µiωA −

κγ iσT
e 1=0
2

O( 2 ) : −D1 2 A − 2D2 Aiω − 2µD1 A −

Îł2 A
4(Ω2 + 2ωΩ)

−

(4.50)

Îł2 A
4(Ω2 − 2ωΩ)

=0

Case 2: Ω = 2ω +  σ
O() :

−2iωD1 A − 2µiωA −

γ Ā iσT
e 1 =0
2

O( 2 ) : −D1 2 A − 2D2 Aiω − 2µD1 A −

Îł2 A
4(Ω2 + 2ωΩ)

(4.51)
=0

Case 3: Ω = ω/2 +  σ
O() :

−2iωD1 A − 2µiωA = 0

O( 2 ) : −D1 2 A − 2D2 Aiω − 2µD1 A −
−

κγ 2
eiσT1 = 0
4(Ω2 − 2ω2 )
71

Îł2 A
Îł2 A
−
4(Ω2 + 2ωΩ) 4(Ω2 − 2ωΩ)

(4.52)

6.5

Stable Respone at Ω/ω=2
Untable Respone at Ω/ω=2

Response Amplitude of the system

6

5.5

5

4.5

4

3.5

3

2.5

2

0

0.5

1

1.5

2

2.5

Frequency Ratio (Ω/ω)

Figure 4.8: Amplitudes of simulated responses of equation (4.45) showing unstable response at
subharmonic resonance due to increase of the parametric forcing amplitude;
 = 0.1, Âľ = 0.25, F0 = 2. Different curves depict Îł = 0.5 and 1.

70

Primary and 1/2 subharmonic resonance
Primary, 1/2 and 1/3 superharmonic resonance

Response Amplitude of the system

60

50

40

30

20

10

0

0

0.2

0.4

0.6

0.8

1

1.2

Frequency Ratio (Ω/ω)

Figure 4.9: Amplitudes of simulated responses of equation (4.45) showing the occurrence of
superharmonic resonances at high parametric forcing amplitude;  = 0.1, Âľ = 0.25, F0 = 2.
Different curves depict Îł = 2 and 4.

72

Response Amplitude of the system

30
Primary Resonance
1/2 Superharmonic Resonance
1/3 Superharmonic Resonance

25

20

15

10

5

0

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

Parametric excitation magnitude (.)

Figure 4.10: Trend of amplitudes of simulated responses of equation (4.45) at different critical
frequencies. Order 2 subharmonic goes unstable at significantly lower values of Îł and hence is
not shown.
Other superharmonic resonance conditions are revealed at higher orders of expansion. For example, the third-order expansion reveals resonance cases at Ω ≈ ω/3 and Ω ≈ 2ω/3. The existence
of resonances and instability zones is governed by the solutions to the solvability conditions for
specific combinations of Ω and ω. The set of equations for the various cases above are set up in the
same way as were the systems of equations for solvability conditions for the harmonically forced
equation (following the procedure detailed in sections 4.1 and 4.2). Next, we take an example case
and sketch out the solution to verify the predictions with numerical plots.

4.3.1

Primary Resonance

After reconstituting the solvability conditions at orders  and  2 in Equations (4.50), we have



γκ iσt 
κγσ iσt i ¾κγ iσt γ 2 A
dA
2
2
+  −2µiωA −
e
+ Âľ A+
e
−
e
+
= 0.
−2iω
dt
2
4ω
4ω
6ω2

73

(4.53)

We let A = (X + iY )ei σt and then split the equation (4.53) into its imaginary and real parts
following the process established in previous sections to arrive at
−2ω

γκ ¾ γ 2Y
dX
+ (−2µωX + 2ωσY ) +  2 (µ2Y −
+
) = 0,
dt
4ω
6ω2

γκ
γκσ Xγ 2
dY
+ (2µωY + 2ωσX −
) +  2 (Âľ2 X +
+
) = 0.
dt
2
4ω
6ω2
At steady state XÛ = YÛ = 0. Hence rearranging equations (4.31) and (4.30) we get
2ω

−2 µωX

+ (2ωσ +  2 µ2 +

 2Îł2
 2 κ ¾γ
)Y =
4ω
6ω2

(4.54)
(4.55)

(4.56)

 2Îł2
 γκ  2 κγσ
−
)X +
2 µωY
=
2
4ω
6ω2
We solve (4.56) for X and Y by using Cramer’s rule, and then write the square of the total amplitude
(2ωσ +  2 µ2 +

r 2 = (X 2 + Y 2 ) as


2

 
γκ  κγσ 
 Îł2
 κ ¾γ 
2
−
(2ωµ)
2ωσ +  µ +
−
2
4ω
4ω
6ω2
!

2 2 2

Îł
(2ωµ)2 + 2ωσ +  µ2 +
6ω2



! 

2
2
 κ ¾γ 

Îł
γκ  κγσ
2
−
2ωσ +  µ +
+
(2ωµ)
4ω
2
4ω
6ω2
!

2 2 2

Îł
(2ωµ)2 + 2ωσ +  µ2 +
6ω2

r2 =

+

(4.57)

From the expression for r 2 for the case of primary resonance (given by equation (4.57)) it is
clear that at two orders of expansion for the differential equation prescribed by (4.45) there is no
instability realized. This is due to the fact that the expression for the determinant as given by the
denominator of expression for r in equation (4.57) is always positive. The values for amplitude at
primary resonance as obtained from equation (4.57) are in agreement with the numerical simulations
of equation (4.45). The analytical prediction for amplitude at Îł = 1 from equation (4.57) is in
agreement with the values from numerical simulations as shown via figure 4.10.
74

4.3.2

Subharmonic Resonance

We follow the same procedure to show the existence of instability at the sub-harmonic resonance
for equation (4.45). The expression for amplitude r 2 and the corresponding trace and determinant
for the reconstituted equations from the solvability conditions at various orders of expansion are as
follows:



2
2    κ ¾γ 
Îł

Îł
 κγσ 
−
2ωσ + +  µ2 +
−
(2ωµ)
4ω
2
4ω
6ω2
!

2 2  Îł 2 2

Îł
(2ωµ)2 + 2ωσ +  µ2 +
−
2
6ω2



! 

2
2
 κ ¾γ 

Îł

κγσ
Îł
−
2ωσ − +  µ2 +
(2ωµ)
4ω
2
4ω
6ω2
!

2 2  Îł 2 2

Îł
(2ωµ)2 + 2ωσ +  µ2 +
−
2
6ω2

r2 =

+

(4.58)

T = −4 µω
 Îł2
D = (2ωµ)2 + 2ωσ +  µ2 +
6ω2

!2
−

 Îł 2

(4.59)

2

From equation (4.59) we can compute a critical value for parametric resonance, beyond which
the system will be unstable at the Ω ≈ 2ω subharmonic, as


18ω2 1

2
−
(2ωσ +  µ )
Îłcr =

4 3ω2
s
! 1/2
2
2



1

Âą
−
(2ωσ +  µ2 ) −
4µ2 ω2 + 2ωσ +  µ2
4 3ω2
9ω4

(4.60)

Beyond the critical value, the system response would be unstable.

4.3.3

Superharmonic Resonance

We cannot reveal instabilities with the case of Ω = ω/2 superharmonic at this level. The expressions
for amplitude at this level are in accord with two orders of expansion, but, to analyze the critical
75

values at which the system would become unstable, we would have to perform higher order
expansions. As such, the instability wedges at these frequencies are very slender and would not
be excited for “reasonable" values of system parameters. As we can see in figure 4.9, if we keep
increasing the value of parametric excitation, lower order superharmonics also show-up in the
response plots. However, if we were to continue increasing Îł at a certain point, the system could hit
unstable wedges corresponding to various frequencies and the response would become unstable. In
this work, we restrict analysis to the dominant harmonics that are observed for reasonable system
values.
The response predictions for systems with direct excitation in combination with the sinusoidal
excitation, exacerbate the response at the resonance conditions thus increasing the likelihood of
resonant depending on the relative magnitude of the terms. We showcased the comparison of
amplitudes at primary resonance for this condition to check the accuracy of predictions. By
superposition, the response of equation (3.1) would be a combination of the response of equations
(4.18) and (4.45), thus increasing the likelihood of large amplitude oscillations.

4.4

Effect of Nonlinearity on Mathieu Resonances and Instabilities

Typically while analyzing the Mathieu equation we look for stability transition curves in the
 − δ (i.e. magnitude of parametric forcing - frequency) space. The transition curves are the Arnold
tongues/Mathieu wedges. In our analysis of the weakly forced linear Mathieu equation, we sought
to reconstruct the Arnold tongue in a similar space and study the effect of direct forcing on the
system. Addition of nonlinearity will impact the stable and unstable solutions arrived at through
the linear analysis. From [75] we know that there are supercritical and subcritical bifurcations that
take place for an unforced nonlinear Mathieu equation as we vary parameters in the  − δ space.
Here we apply second-order multiple scales to the nonlinear Mathieu equation (3.1). Inserting
the expansion (4.2) into equation (3.1) where µ̂ = 2µ, and equating coefficients of like powers of 
yields

76

O(1) :

D0 2 q0 + ω2 q0 = 0

O() :

D0 2 q1 + ω2 q1 = −2µD0 q0 − 2D0 D1 q0 − γq0 cos ΩT0 − αq0 3 + F sin ΩT0

O( 2 )

D0 2 q2

:

+ ω2 q2

2

= −2D0 D1 q1 − (D1 + 2D0 D2 )q0 − 2µ(D0 q1 + D1 q0 )

(4.61)

−γq1 cos ΩT0 − 3αq0 2 q1
The solution of O(1) equation is

q0 = A(T1, T2 )eiωT0 + c.c.

(4.62)

We substitute this into the O() equation and identify other resonance conditions to eliminate
secular terms and seek the solution of q1 .
The equation for q1 obtained by substituting the solution for q0 into O() equation is

D0 2 q1 + ω2 q1 = −2iωD1 AeiωT0 − 2µiωAeiωT0 − α(A3 e3iωT0 + 3A2 ĀeiωT0 )
Ā
iF
A
−γ( ei(ω+Ω)T0 + ei(Ω−ω)T0 ) − eiΩT0 + c.c.
2
2
2

(4.63)

We have three cases for equation (4.63) which can contribute towards resonance condition.
1. No specific relation between Ω and ω
2. Ω ≈ ω
3. Ω ≈ 2ω
Case 1: When there is no specific relationship between the forcing frequency Ω and the natural
frequency ω we equate the secular terms to zero, such that
−2iωD1 A − 2µiωA − 3αA2 Ā = 0
and solve the remaining ODE in equation (4.63) to get the particular solution, by treating A as
constant with respect to the independent variable T0 , as
77

αA3 3iωT
ÎłA
0+
q1 =
e
ei(Ω+ω)T0
2
2
8ω
2(Ω + 2ωΩ)
(4.64)
iF
γ Ā
i(Ω−ω)T
iΩT
0+
e
e 0 + c.c.
+
2(Ω2 − 2ωΩ)
2(Ω2 − ω2 )
Case 2: For the second condition of when Ω ≈ ω i.e. Ω = ω +  σ1 (σ1 being the detuning
parameter for this case) from equation (4.63),
−2iωD1 A − 2µiωA − 3αA2 Ā −

iF iσ T
e 1 1 =0
2

forms the solvability condition. The particular solution for q1 then becomes

q1 =

αA3 3iωT
ÎłA
γ Ā
0+
e
ei(Ω+ω)T0 +
ei(Ω−ω)T0 + c.c.
8ω2
2(Ω2 + 2ωΩ)
2(Ω2 − 2ωΩ)

(4.65)

Case 3: And finally, when Ω ≈ 2ω i.e. Ω = 2ω +  σ2 (σ2 being the detuning parameter for this
case)
−2iωD1 A − 2µiωA − 3αA2 Ā −

γ Ā iσ T
e 2 1=0
2

forms the solvability condition. The particular solution for q1 then becomes
q1 =

iF
αA3 3iωT
ÎłA
0+
ei(Ω+ω)T0 +
eiΩT0 + c.c.
e
2
2
2
2
8ω
2(Ω + 2ωΩ)
2(Ω − ω )

(4.66)

We now substitute the appropriate solutions for q0 and q1 into the expression (4.61) at O( 2 ).
After expanding the terms in the expression for the q2 ODE, we seek to identify terms that can
contribute to a resonance condition.
For Case 1 when there is no relation between Ω and ω, the expression for q2 is:

78

D0 2 q2

+ ω2 q2



= −2 i(Ω + ω)D1 Kei(Ω+ω)T0 + i(Ω − ω)D1 Lei(Ω−ω)T0
 

3iωT
iΩT
iωT
iωT
2
0 + iΩD1 Ne
0 − D1 Ae
0 + 2iωD2 Ae
0
+3iωD1 Me

−2µ i(Ω + ω)Kei(Ω+ω)T0 + i(Ω − ω)Lei(Ω−ω)T0 + 3iωMe3iωT0 + iΩNeiΩT0


Îł
iωT
0
+D1 Ae
− Kei(2Ω+ω)T0 + KeiωT0 + Lei(2Ω−ω)T0 + Le−iωT0
2


i(Ω+3ω)T
i(−Ω+3ω)T
i2ΩT
0 + Me
0 + Ne
0 + N − 3α A2 Kei(Ω+3ω)T0
+Me

(4.67)

+ Ā2 Kei(Ω−ω)T0 + 2A ĀKei(Ω+ω)T0 + A2 Lei(Ω+ω)T0 + Ā2 Lei(Ω−3ω)T0
+2A ĀLei(Ω−ω)T0 + A2 Me5iωT0 + Ā2 MeiωT0 + 2A ĀMe3iωT0 + A2 Nei(Ω+2ω)T0

i(Ω−2ω)T
iΩT
2
0 + 2A ĀNe
0 + c.c.
+ Ā Ne
where K =

γ Ā
ÎąA3
iF
ÎłA
,L=
,M=
.
,N=
2
2
2
2
2(Ω + 2ωΩ)
2(Ω − 2ωΩ)
8ω
2(Ω − ω2 )

Following the analysis done in [74] for the unforced, linear case, we write the solvability
conditions together for various resonance cases. Once we enlist the solvability conditions for each
of the resonance case identified, we look for the fixed points and stability for each case separately.
For the case of no specific relationship between ω and Ω

O() :

−2iωD1 A − 2µiωA − 3αA2 Ā = 0

O( 2 ) : −D1 2 A − 2D2 Aiω − 2µD1 A −

3α2 A3 Ā2
Îł2 A
Îł2 A
−
−
=0
8ω2
4(Ω2 + 2ωΩ) 4(Ω2 − 2ωΩ)
(4.68)

There are various other combinations of Ω and ω appearing at second-order that could lead to
resonance conditions in the system. Examining the O( 2 ) equation given by equation (4.67) for the
general case we arrive at the following combinations:
1. Ω ≈ 3ω
2. Ω ≈ ω/2
79

3. Ω ≈ 4ω
The other resonance combinations identified at O() i.e.

4. Ω ≈ ω
5. Ω ≈ 2ω

lead to slightly different expressions for an ODE in q2 based on the expression of q1 computed
for each case and given in equations (4.65) and (4.66). We examine the resulting expression for
solvability conditions at O().
For each of these conditions, the corresponding equations at O() and O( 2 ) are given below.

Ω ≈ 3ω
O() :

−2iωD1 A − 2µiωA − 3αA2 Ā = 0

O( 2 ) : −D1 2 A − 2D2 Aiω − 2µD1 A −
−

Îł2 A
3α2 A3 Ā2
−
8ω2
4(Ω2 + 2ωΩ)

(4.69)

Îł2 A
i3αF Ā2 iσT
−
e 1 =0
4(Ω2 − 2ωΩ) 2(Ω2 − ω2 )

Ω ≈ ω/2
O() :

−2iωD1 A − 2µiωA − 3αA2 Ā = 0

O( 2 ) : −D1 2 A − 2D2 Aiω − 2µD1 A −
−

3α2 A3 Ā2
Îł2 A
−
8ω2
4(Ω2 + 2ωΩ)

ÎłFi
Îł2 A
−
eiσT1 = 0
2
4(Ω − 2ωΩ) 4(Ω2 − ω2 )

Ω ≈ 4ω

80

(4.70)

O() :

−2iωD1 A − 2µiωA − 3αA2 Ā = 0

O( 2 ) : −D1 2 A − 2D2 Aiω − 2µD1 A −
−

3α2 A3 Ā2
Îł2 A
−
8ω2
4(Ω2 + 2ωΩ)

(4.71)

Îł2 A
3αγ Ā3
γα Ā3 iσT
−
eiσT1 −
e 1 =0
4(Ω2 − 2ωΩ) 2(Ω2 − 2ωΩ)
16ω2

Ω≈ω
O() :

−2iωD1 A − 2µiωA − 3αA2 Ā −

O( 2 ) : −D1 2 A − 2D2 Aiω − 2µD1 A −
−

iF iσT
e 1 =0
2

3α2 A3 Ā2
Îł2 A
−
8ω2
4(Ω2 + 2ωΩ)

(4.72)

Îł2 A
γ 2 Ā
−
e2iσT1 = 0
2
2
4(Ω − 2ωΩ) 4(Ω − 2ωΩ)

Ω ≈ 2ω
O() :

−2iωD1 A − 2µiωA − 3αA2 Ā −

O( 2 ) : −D1 2 A − 2D2 Aiω − 2µD1 A −

γ Ā iσT
e 1 =0
2

3α2 A3 Ā2
Îł2 A
−
8ω2
4(Ω2 + 2ωΩ)

(4.73)

3αγ A Ā2
eiσT1 = 0
2(Ω2 + 2ωΩ)
The σ ’s above are the detuning parameter in each case respectively.
−

In this analysis, using weak forcing, we do not capture the linear superharmonic at order 1/3.
This happens because of the way the terms and the bookkeeping parameter  appear in the equations.
The terms are multiplied at varying orders, thus burying the effect of that superharmonic “one-order
lower". If we were to continue our analysis further we will encounter that resonance term.
The solutions for the set of equations given for each resonance case above would give us the
slow time scale variations of the amplitude A of our fast time scale solution q. In order to arrive
at a solution for A we re-combine the terms at O(), O( 2 ) into a single ODE and seek its solution.
To arrive at this juncture would require elimination of D1 A terms from O( 2 ) equation.

81

We sketch the treatment of one of the resonant cases listed. We consider the terms for the
primary resonance case Ω ≈ ω given in equation (4.72). From the O() equation we get,
3αi A2 Ā FeiσT1
−
D1 A = −µA +
2ω
4ω
and
D1 Ā = −µ Ā −

3αi Ā2 A Fe−iσT1
−
.
2ω
4ω

We compute D1 2 A from the above expressions as
Fiσ iσT 3αFi A Ā iσT 3αFi A2 −iσT
6αi µA2 Ā 9α2 A3 Ā2 F µ iσT
1.
+
e 1−
e
D1 2 A = µ2 A −
−
e 1−
e 1−
ω
4ω
4ω
4ω2
4ω2
8ω2
We substitute these in the O( 2 ) expressions in equation (4.72) to arrive at

6αi µA2 Ā 9α2 A3 Ā2 F µ iσT
Fiσ iσT
3αFi A Ā iσT
−2D2 Aiω + − µ2 A +
+
e 1+
e 1+
−
e 1
2
2
ω
4ω
4ω
4ω
4ω



3αFi A2 −iσT
3αi A2 Ā FeiσT1
3α2 A3 Ā2
Îł2 A
1
−
+
e
− 2µ − µA +
−
−
2ω
4ω
8ω2
8ω2
4(Ω2 + 2ωΩ)
2
2
γ Ā
Îł A
−
e2iσT1 = 0
−
2
2
4(Ω − 2ωΩ) 4(Ω − 2ωΩ)
(4.74)


It can easily be shown that the O() solvability condition from equation (4.72) and the O( 2 )
solvability condition from equation (4.74) result from a multiple-scales expansion of




dA
iF
3αFi A2 −iσt γ 2 Ā 2iσt
−2iω
+  −2µiωA − 3αA2 Ā − eiσt +  2 µ2 A +
e
+
e
dt
2
8ω2
4ω2

9αi µA2 Ā γ 2 A 15α2 A3 Ā2 3αFi A Ā iσt Fiσ iσt F µ iσt
−
+
+
+
e
+
e
+
e
=0
ω
4ω
4ω
6ω2
32ω2
4ω2

(4.75)

Transforming the system to polar coordinates, leads to steady state equations in a and φ =
σT1 − β. The steady state equations include sin(φ), cos(φ) and sin(2φ), cos(2φ) terms, which make
it more difficult to solve than other familiar examples.

82

We follow the solution procedure given in [74] and detailed in the previous chapter. We seek a
solution in the form A = (Br + iBi )eiσt , with real Br and Bi . We enforce this solution in equation
(4.75), separate real and imaginary parts and cancel the common exponential term to obtain,



3ÎąF
Îł2
dBi
3
2
2
Br Bi +
Br
2ω
+  2µωBi + 2ωσBr − 3α(Br + Br Bi ) +  µ2 Br −
dt
4ω2
4ω2
(4.76)

2
Îł2
9ι¾ 3
F
Âľ
15Îą
+
Br +
(Br 5 + Bi 4 Br + 2Bi 2 Br 3 ) +
(Bi + Br 2 Bi ) +
= 0,
ω
4ω
6ω2
32ω2

 F 
dBr
3Îą
3
2
+  −2µωBr + 2ωσBi − 3α(Bi + Bi Br ) +
1+
(Br 2 − Bi 2 )
−2ω
dt
2
4ω2


Îł 2 Bi
σ

Îł
9ι¾ 3
15Îą2 5
2
+
+ −
+ Âľ2 Bi +
Bi −
(Br + Bi 2 Br ) +
(Bi + Br 4 Bi
2
2
2
2ω
ω
4ω
6ω
32ω


(4.77)

+2Br 2 Bi 3 ) = 0.

As stated in the previous section, we have to solve for Br and Bi to get relations between
γ, F, µ, , σ and ω. The fifth degree terms still pose a challenge. If we proceed with enforcing
polar coordinate forms to the Bi, Br terms in the equation then we are faced with sin(φ), cos(φ) and
sin(2φ), cos(2φ) terms.
Similar analysis has been done for all of the resonant conditions that exist in the nonlinear
system. Here, as was in the case with the linear system, we are left with a differential equation
with higher order polynomial terms. The origin is not a solution for this set of equation as we have
direct forcing.
Having obtained these equations, we look for stability characteristics based on system parameters. This will require a lot more analysis before we can draw conclusions on stabilities and the
parameters that dictate behavior. The influence of parametric excitation on a linear version of
equations (4.76) and (4.77) has been demonstrated in the previous section. The completed form
of the expressions for the nonlinear system will show the explicit dependence on the amplitude on
parametric excitation and the nature of the stability of the numerous operating points that arise as
steady-state solutions of the system. The second order analysis will also provide expressions for
the local and global stability of the system. The current work, lays the foundation for constructing
the stability transition curves and analyzing parameter zones where in the solution would either
83

Table 4.1: Stability Wedge and Resonance Chart. R1 : Resonance identified at 1st order of MMS
expansion. R2 : Resonance identified at 2nd order of MMS expansion. W2 : Instability wedge
(Arnold tongue) expression can be found at second order of MMS expansion. −: Known
resonance case/ Instability not uncovered up to two orders of expansion
Forcing (Ω)
2ω
ω
2ω/3
ω/2
2ω/5
ω/3

O(1) Forcing
R1, W2
R1, W2
R1, W2
R2, W2

O() Forcing
R2
R1, W2
R1, W2
-

No Forcing
R1, W2
R2, W2
-

destabilize or be resonant. Both would lead to sustained oscillation in the system thereby causing
increased loading in our wind turbine system which has been modeled using these equations. Table
4.1 lists the frequency ratios at which resonances or wedges have been identified.

4.5

Summary

The main findings reported in this chapter show that a linear forced Mathieu equation has
a homogenous and particular solution. When direct and parametric excitations are of the same
frequency, the system can exhibit superharmonic resonances, subharmonic resonances, and primary
resonance. A second-order perturbation analysis uncovers superhamonics at orders 1/3 and 1/2, and
subharmonic instability at order 2. It also reveals a primary resonance parametric amplification,
which becomes unbounded at a critical value of parametric excitation level Îłcr , at which the
homogenous solution also becomes unstable. When the direct forcing is a constant load, the system
can have superhrmonic, subharmonic, and primary resonances.

84

CHAPTER 5
APPLICATION TO WIND TURBINE BLADES

5.1

Introduction

The most essential feature of the wind turbines are the blades, which harness the power from the
wind. The blades are expected to be in operation through the entire design life and often undergo
extreme loading. Environmental (wind and other), electrical (grid) and soil conditions (foundation)
affect loading, durability and operation of wind turbines and hence should be taken into account
in the design process [12]. Each external environmental condition is typically classified into a
normal or extreme operating condition with typical occurrence intervals stated for design purpose.
For example, certifications of wind turbine blades consider operation of the equipment not only in
normal and turbulent wind conditions but also in extreme conditions such as gusts and direction
changes. Despite these efforts wind turbine blade failures are common and their maintenance
requires significant capital expenditure. As these blades grow larger in size, the dynamic loads
experienced is not well understood. In chapters 2, 3 and 4 we developed the equations of motion
and perturbation analysis necessary to capture and predict the behavior of the system of equations
that define the dynamic motion of these systems. In this chapter, we extend the application to actual
data from wind turbine blades.

5.2

Numerical Simulations of ODE’s

In this section, we address the response of the system described by equations (2.8) and (2.16).
These describe the vibration motion of the wind turbine blade and account for nonlinearities and
variations in geometry along the length of the blade. We dealt with a few terms in these equations
that would dominate the response of the system. These were obtained by assuming small modal
coordinates (i.e. q =  q̂) and eliminating terms that have higher orders of . Since the cross-coupled
terms contribute at least a 2nd or 3r d power of  they can be safely ignored in the calculations of

85

Figure 5.1: Response amplitude of system described by equation (2.8) that described in-plane
motion of the blade. Values used in simulation
b = 0.1, d = 0.1, e1 = 1, e2 = 0.3, f1 = 0.1, f2 = 0.1, F = 1. To account for damping we add a 2µqÛ
term with Âľ = 0.25
the resonances and stabilities for the forced nonlinear Mathieu representation of this system.
Figure 5.1 shows the response of the system described by the ODE in equation (2.8). The
superharmonic responses at orders 1/3 and 1/2 are those that have been captured in the simplification of this ODE. Thus, the fundamental behavior of the system is governed by the forced
nonlinear Mathieu equation. This also holds true for the out-of-plane resonance case. We then
proceed to quantify motion for actual wind turbine blade parameters. Before we do this, we need
a representative aerodynamic force that accounts for external loads due to the wind. We apply the
BEM theory to arrive at an expression for aerodynamic load. Commercial codes like FAST have
also been used to assess blade vibration, tip deflection, etc.

5.3

Aerodynamic Forces on the Blade

We estimate the aerodynamic loads under the following assumptions:
• Constant angle of attack
86

• Steady wind
• Constant rate of rotation
• Linear wind shear
• Uniform wind direction
Blade element momentum theory gives a reasonably good first estimates for aerodynamic
loading on the blades for analysis from a dynamic system perspective. The magnitude of forces
computed using this method are only used to get a coarse estimate as it relies on two key assumptions;
(i) there are non interactions between the different blade elements and (ii) forces on the blade
elements are determined by the lift and drag forces on the undeformed blade. More detailed models
accounting for dynamic stall or wake vortices may contribute additional terms and would have to
be considered depending on the problem definition. However, for resonance analysis, the BEM
theory provides reasonably good lift and drag forces along the length of the blade.
The lift force (L) and drag force (D) per unit length is:
1
ρcW 2CL
2
1
D = ρcW 2CD
2
L=

(5.1)

where c is chord length, ρ is air density, W is relative velocity at blade, CL and CD are coefficients
of lift and drag. Since we are interested in the forces acting normal to and tangential to the plane of
rotation (out-of-plane and in-plane respectively) the lift and drag forces can be projected in these
directions (shown as FN and FT in figure 5.2). The aerodynamic loads for out-of-plane and in-plane
blade motion can thus be written as:
FN = L cos φ + D sin φ
FT = L sin φ − Dcosφ
where L and D are from equation (5.1).

87

(5.2)

decrease of lift. The lift and drag coecients need to be projected onto the NTdirection.
CN = CL cos φ + CD sin φ
(3.7)
and
CT = CL sin φ − CD cos φ

(3.8)

L

Vrel

φ

FT
φ

R

90◦ − φ

FN

Rotor plane

D

Figure 3.1: The local forces on the blade. Redrawn from [31].

Figure 5.2: The local forced on an aerosol structure of the blade, reproduced from [3]
23

The information on the coefficients of lift and drag are required to calculate the in-plane and
out-of-plane loading and these depend on the angle of attack, Îą. Further, from figure 5.3, we can
write the relation between mean free stream wind velocity U∞ and angle of attack as

tan φ =

(1 − a)U∞
(1 + a0)ωr

(5.3)

where a and a0 are axial and tangential induction factors. From equations (5.1), (5.2) and (5.3)
we can write expression for aerodynamic load at every section of the blade and compute the
aerodynamic loading over the entire blade. Typically, estimation of a and a0 requires an initial
guess to calculate the aerodynamic forces and iterating through the solution at every step until a
and a0 do not change beyond a defined tolerance. This standard procedure is described in detail
in [43]. We follow the procedure laid out there for an example blade in the next section.
The initial estimate for the aerodynamic forces from the BEM theory can be improved by adding
blade solidity ratio, wind shear, tip loss correction factors, etc. We incorporate wind shear into
this formulation. Wind shear is the variation of mean wind speed as we move vertically along the
height of the wind turbine. The mean wind speed can vary significantly along the length of the
blade as the blade rotates through the rotor plane. The mean wind speed at certain height is given
88

CHAPTER 3. WIND TURBINE DESIGN CALCULATIONS

U∞ (1 − a)
Îą
Vrel

θ

ωr(1 + a′ )

Rotor plane

φ

Figure 3.2: Velocities at the rotorplane. Redrawn from [31].

Figure 5.3: Velocities at the rotor plane, reproduced from [3]
Further, the solidify σ is defined as the fraction of the annular area in the control
volume, which is covered by the blades

with respect to the mean wind speed at height 10c(r)N
m above ground by
σ(r) =

2πr



(3.9)

Îş

h
U = U10
The normal force and the torque on the controlh10
volume of thickness dr, is since FN
where N denotes the number of blades.
and FT are forces per length
2 from 0.01
where Îş is based on the terrain and typically1varies
to 0.7 depending on the location of
U∞
(1 − a)2

dT = N FN dr = ρN
2

sin2 φ

cCN dr

(3.10)

wind turbine installations. This curve can also be approximated by simpler function for each blade
and

as

1
U∞ (1 − a)ωr(1 + a′ )
dQ = rN FT dr = ρN
cCt rdr
2
sin φ cos φ
Finally, the two induction factors are defined as
1
U(h)
+ h ∗ w+
a == Ubase
2
4 sin φ
+1
σCN

(3.11)

(3.12)

where w + isand
the approximate rise in wind speed for an increase in height h. Tower shadowing also
1
(3.13)
4
sin
φ
cos
causes a cyclic excitation for each tower pass. There φare
− 1numerous models that quantify the effect of
σCT
a′ =

tower shadow
the wind
turbinehave
blade
forbeen
every
pass (for
e.g.,BEM
Powles’
model
that
All on
necessary
equations
now
derived
for the
model.
Since
theis implemented
different control volumes are assumed to be independent, each strip may be treated

in GH Bladed,
Blevin’s
andthe
theresults
JET wake
The
look at modeling the
separately
andmodel
therefore
for onemodel).
radius can
besimplest
computedmodels
before solving
for another one. For each control volume, the algorithm can be divided into eight

loss of windsteps:
speed and these can be written as a summation of sinusoidal functions at the frequency
of rotation and multiples of frequency. Since we are modeling only one blade, the frequency of
24

tower pass is the same as the frequency of rotation. This will be altered appropriately once we have
∍ L
0
coupled blades. The total aerodynamic load on the blade is compute as Q =
Q aero dx, where
0

Q aero is the load per unit length.
89

solution. Because the full load was applied statically, satisfaction of the buckling requirement is met
when buckling mode frequencies are computed to be greater than the required safety factor. This is
because the buckling frequencies are equal to the scaling factor of the applied static load at which a
particular buckling mode will occur.

Figure 5.4: Screenshot of NuMAD blade geometry, reproduced from [4]

Figure 17. Screenshot of NuMAD Blade Geometry
5.3.1 Sandia 100 m Blade
47
In order to compute the total aerodynamic load we must know the exact geometry of the blade. We
go through an example calculation of the loads on Sandia National Lab’s 100 m blade. These blades
are currently in development to study the scope of usage of large offshore wind turbine blades and
the data is used extensively for research and development. A NuMAD model of the blade used is
shown in figure 5.4. NuMAD (Numerical Manufacturing And Design) is an open-source software
tool written in Matlab that enables users to create a three-dimensional model of a wind turbine
blade. The tool manages all blade information including airfoils, materials, and material placement.

As we can see from figure 5.4, the airfoil shape and geometry changes significantly as we go
along the length of the blade. There are 34 sections at which information about blade geometry,
material, profile, etc., is provided for this blade as tabulated in tables 5.1 and 5.2.
For a given given wind speed and blade velocity, the angle of attack/incidence on an airfoil
changes. The Sandia 100 m blade has various airfoil profiles along the length of the blade as
indicated in table 5.1. Figure 5.5 shows the coefficients of lift and drag for two sample airfoils in

90

Table 5.1: Blade fraction align the station, chord length,twist and airfoil description for the Sandia
100m blade
Station Number
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

Blade Fraction
0
0.005
0.007
0.009
0.011
0.013
0.024
0.026
0.047
0.068
0.089
0.114
0.146
0.163
0.179
0.195
0.222
0.249
0.276
0.358
0.439
0.52
0.602
0.667
0.683
0.732
0.764
0.846
0.894
0.943
0.957
0.972
0.986
1

Chord (m)
5.694
5.694
5.694
5.694
5.694
5.694
5.792
5.811
6.058
6.304
6.551
6.835
7.215
7.404
7.552
7.628
7.585
7.488
7.347
6.923
6.429
5.915
5.417
5.019
4.92
4.621
4.422
3.925
3.619
2.824
2.375
1.836
1.208
0.1

91

Twist (deg)
13.308
13.308
13.308
13.308
13.308
13.308
13.308
13.308
13.308
13.308
13.308
13.308
13.308
13.177
13.046
12.915
12.133
11.35
10.568
9.166
7.688
6.18
4.743
3.633
3.383
2.735
2.348
1.38
0.799
0.28
0.21
0.14
0.07
0

Airfoil Description
Cylinder
Cylinder
Transition (99.25%)
Transition (98.5%)
Transition (97.75%)
Ellipse (97%)
Ellipse (93.1%)
Ellipse (92.5%)
Transition (84%)
Transition (76%)
Transition (68%)
Transition (60%)
Transition (51%)
Transition (47%)
Transition (43.5%)
DU99-W-405
DU99-W-405 (38%)
DU99-W-350 (36%)
DU99-W-350 (34%)
DU97-W-300
DU91-W2-250 (26%)
DU93-W-210 (23%)
DU93-W-210
NACA-64-618 (19%)
NACA-64-618 (18.5%)
NACA-64-618
NACA-64-618
NACA-64-618
NACA-64-618
NACA-64-618
NACA-64-618
NACA-64-618
NACA-64-618
NACA-64-618

Figure 5.5: Variation of coefficient of lift and drag for the NACA 64 airfoil and DU30 airfoil
the blade. The blade mass is 114,172 kg, maximum rotational speed is set to 7.44 rpm and the
rated wind speed is 11.3 m/s. Using the aeroelastic code FAST (models available for download
at [6]), the tip deflection, flap-wise and edge-wise moments at the rotor hub can be estimated for
various operating conditions. It is typical to check for loading under most naturally occurring
wind conditions as defined by IEC-61400. These are EWM (extreme wind model), EOG (extreme
operating gust), ETM (extreme turbulent model), EDC (extreme direction change), EWS (extreme
wind shear), etc. We use the loading conditions derived from the FAST model as aerodynamic
loading for our ODE in the out-of-plane direction. From the values predicted for the moment at the
root, we calculate the forces acting on the blade necessary to produce the given moment.
This will give rise to a distributed loading information on the wind turbine blade due to wind.
This in combination with gravity are the most dominant loads on the wind turbine blades. We only
focus on routine operational ranges, as it is in our interest to investigate potential resonances away
from the natural frequency of the blade (viz., super- and sub-harmonic resonances) for standard
loading. This would lead to abnormal loads on the entire structure of the wind turbine and all the
machinery components leading to premature failure.
We use the information from the tables and curves prescribed for the Sandia 100 m blade and
compute the lift and drag forces at various sections of the blade.

92

5.4

Response of Wind Turbine Blades

The response of different wind turbine blades are studied. Data from the Sandia 9 m CX blade,
NREL 23.5 m GRC blades, a commercial 40 m and 63 m blade, and a concept 100 m blade are used
to study the influence of direct and parametric excitation on out-of-plane blade response. Based on
analysis in [54, 87], the out-of-plane motion has significant response to cyclic load variations.
We compute the coefficients of equation (2.16) restated below based on the expression given in
equation (2.17).
Ü q+b
Û qq
Ü 2 +d qÛ2 q+(e1 +e2 cos Ωt)q+( f1 + f2 cos Ωt)q3
q+Âľ
• Linear co-efficient of q: e1 =

=

1
(CL cos φ+CD sin φ) ρW 2 c dr
2
0

∍ L

12.36E I
+ 1.193Ω2
3
mL

121.32E I .926Ω2
+
• Nonlinearity: f1 = −
L2
mL 5
f2 =

39.66
L3

• Cross-coupled terms: b = d =
• Parametric excitation: e2 =

4.59
L2

15.41
L

• Direct excitation: Q aero
The values of the coefficients of the equation of motion that describe the out-of-plane motion of
the blade are dependent on blade elasticity and moment of inertia E I, mass M and length L. The
computation of the aerodynamic terms is dependent on the evaluation of lift and drag coefficients
for each blade. The forces are then integrated along the length of the blade to obtain the net force.
The aerodynamic term essentially bottles down to three terms: a constant term, a cyclically varying
term.
We then simulate the response of the 23.5 m blade shown in figure 5.6. The response of the
23.5 m blade does not show any superharmonic resonances. Commercial blade data for other
wind turbine systems in operation were collected (for e.g. the 40 m blade) and the responses of
93

Figure 5.6: Amplitude response of the NREL GRC blade, blade properties from [5]
those blades also showed no harmonic response. We investigated the condition under which the
resonance and instability conditions are exhibited from the simplified analysis in Chapter 4. We
notice that as the magnitude of the the parametric excitation increases in relation to the direct
forcing the system starts to exhibit the dynamic behavior inherent in the Mathieu equation analyzed
with “choice parameters". The relative magnitude of the terms in the ODE for an example system
can be tuned to demonstrate this behavior, but, for a wind turbine blade, the blade properties govern
the relative magnitude.
Ishida et al. [55] have shown the effect of varying the magnitudes direct and sinusoidal forcing
for for a rotating beam. The constant component and the cyclic varying term of Q have equal
contributions to the magnitude of the response at the superharmonic resonances. Figure 5.7 shows
a sample response of a 60 m blade with constant and cyclic variation in loads. We can see
that as the blades get larger, the superharmonic resonances continue to persist. Wind turbine
blades are typically designed to operated well below the natural frequency. Let us consider two
examples of blades we have studied and qualitatively study the likelihood of blades operating under
superharmonic resonance conditions. For the NREL GRC 23 m blade, the lowest natural frequency
94

Figure 5.7: Response of a linear system with cyclic loading (left) and constant loading (right),
based on parameter from a 60 m blade.
is 1.72 Hz. The rotational speed of the blade is limited to 22 rpm for rated power. This translates to
about 1/5 of the natural frequency of the rotating blades. The Sandia 100 m blade has the lowest
natural frequency of 0.49 Hz and a the rotation is set to 7.44 rpm for rated power. This translates
to about 1/4 of the natural frequency of the blade. If the rotor speed were to slightly increase due
to other phenomenon, these could operate near the 1/3 superharmonic.
In the formulations of out-of-plane equations of motion, we have only considered one blade with
one projected mode. If we were to include higher modes then the system becomes less stiff. Also,
we could include other blades to account for blade-to-blade dynamic interaction. We see from the
preliminary calculations that the larger blades tend to operate close to the superharmonic captured
in our analysis. There are other resonances that occur for example at much lower frequencies.
However, the response amplitude increase due to them, is orders of magnitude lower than those
captured in the second-order multiple scales expansions. Also, blades tend to be operated at speeds
above their design recommendation due to various reasons (stall regulated blades, slow pitch
control, unexpected wind flow, system requirements, etc.). The work in this thesis should guide
blade designers and wind turbine manufacturers towards avoiding operation under superharmonic
resonance conditions to protect the blades. The increased vibratory loads on blades have a direct
95

impact on the reliability of components in the gearbox. We will address this in the next chapter.
The relative magnitude of the terms in the differential equation of motion are critical towards
triggering responses that lead to instability of increased amplitudes at harmonic frequencies. The
amplitudes of direct and parametric responses are especially of critical importance as has been seen
with analysis of the in-plane equation of motion. Hence, in our ODE given by equation (2.16), the
magnitude of e2 relative to e1 and Q determine the onset or presence of instabilities and resonances.
We plot the ratio of the magnitude of parametric excitation to the square of the natural frequency for
three wind turbine blades of lengths 23.5 m, 63 m and 100 m (see figure 5.8). The blade material
properties are obtained from data sheets for each wind turbine blade (an example data sheet for the
100m blade is shown in Table 5.2, where, flap stiffness is E Izz , edge stiffness is E I yy , and torsional
stiffness is GI x x and x, y, z are axes in the axial, flapwise, and edgewise directions, E and G are
Young’s and shear moduli, and I is second moments of area. As the ratio of e2 /e1 increases, the
response of the system at the harmonics steadily increases. After a certain value this will breach
the stability boundary at these resonances and will start to exhibit unstable response (similar to
the stability transition in the traditional Mathieu equation). The magnitude of the parametric to
direct forcing is also critical to ensure that the system response is stable. An increase in the ratio
of the parametric to direct forcing will also result in an increase in the amplitude of blade motion
at these harmonic frequencies leading to possibility of amplification of amplitude of response and
eventually unstable behavior. Designers, with appropriate choice of airfoil that has a direct impact
on the lift and drag forces on the system, can tune the response to be away from these resonances,
provided they are aware of these phenomena.
Careful attention needs to be paid to the dynamics characteristics of the large wind turbine
blades. As we see from figure 5.8, the response characteristics of the wind turbine system are
altered as the blades become larger. The 100 m blade under study at Sandia National Laboratory is
an “upscaled" version of the 23.5 m blade. If we define the scaling factor k as the ratio of the lengths
of the blades, then the other critical properties are scaled as follows: mass → k 3 , aerodynamic loads
→ k 3 , power → k 2 and gravitational force → k 4 . This implies that we pay a huge penalty in terms

96

Table 5.2: 100m blade structural properties
Span Fraction
(-)
0
0.005
0.007
0.009
0.011
0.013
0.024
0.026
0.047
0.068
0.089
0.114
0.146
0.163
0.179
0.195
0.222
0.249
0.276
0.358
0.439
0.52
0.602
0.667
0.683
0.732
0.764
0.846
0.894
0.943
0.957
0.972
0.986
1

Flap Stiffness
(N-m2 )
3.22E+11
2.88E+11
2.50E+11
2.12E+11
1.75E+11
1.62E+11
1.45E+11
1.31E+11
1.04E+11
7.99E+10
7.05E+10
5.58E+10
4.97E+10
4.76E+10
4.17E+10
3.88E+10
3.36E+10
2.85E+10
2.31E+10
1.45E+10
8.96E+09
5.38E+09
3.11E+09
1.80E+09
1.55E+09
1.01E+09
7.36E+08
3.55E+08
2.08E+08
8.10E+07
4.29E+07
1.92E+07
5.13E+06
5.84E+03

Edge Stiffness
(N-m2 )
3.22E+11
2.88E+11
2.52E+11
2.15E+11
1.80E+11
1.63E+11
1.54E+11
1.41E+11
1.24E+11
1.04E+11
9.37E+10
8.34E+10
8.76E+10
9.85E+10
1.05E+11
1.05E+11
1.03E+11
9.91E+10
6.65E+10
5.48E+10
3.37E+10
2.18E+10
1.44E+10
1.14E+10
1.07E+10
8.72E+09
7.41E+09
4.90E+09
3.66E+09
1.74E+09
9.07E+08
4.19E+08
1.19E+08
6.22E+04

97

Torsional Stiffness
(N-m2 )
1.68E+11
1.50E+11
1.29E+11
1.09E+11
8.97E+10
7.98E+10
7.22E+10
6.47E+10
5.01E+10
3.53E+10
2.54E+10
1.50E+10
9.27E+09
7.66E+09
5.89E+09
4.41E+09
3.85E+09
3.22E+09
2.67E+09
1.72E+09
1.10E+09
6.90E+08
4.32E+08
2.94E+08
2.64E+08
2.10E+08
1.86E+08
1.29E+08
9.78E+07
4.51E+07
2.58E+07
1.15E+07
3.03E+06
2.57E+03

Mass per Length
(kg/m)
5.74E+03
5.09E+03
4.44E+03
3.78E+03
3.15E+03
2.88E+03
2.91E+03
2.68E+03
2.33E+03
1.98E+03
1.84E+03
1.63E+03
1.68E+03
1.74E+03
1.76E+03
1.81E+03
1.79E+03
1.76E+03
1.58E+03
1.43E+03
1.26E+03
1.11E+03
9.35E+02
7.79E+02
7.40E+02
6.07E+02
5.04E+02
3.59E+02
2.90E+02
2.10E+02
1.56E+02
1.20E+02
7.91E+01
6.55E+00

Ratio of coecients of q term e2/e1

0.07

0.06

0.05

0.04

0.03

0.02

0.01
20

30

40

50

60

70

80

90

100

Length of Wind Turbine Blade (m)

Figure 5.8: Ratio of parametric excitation to the square of natural frequency as a function of
length of the blade
of the loads and mass for increase in power of the order of two. The increased mass and loading
coupled with large diameters of blades exposed to different wind velocities can lead to significant
challenges in capturing the physics that govern the flow of wind around the turbine blades and
consequently, our estimates of loads transferred to the wind turbine blades will be incorrect. Also,
as blades get larger there is a reducing factor between the blade rotation speed and the flutter speed
(as seen in figure 5.9).
Hence, even though there are significant benefits towards designing larger wind turbine structures, care must be taken to perform a detailed analysis during the concept and design phase to
account for all the loading conditions. As blades are designed to be longer, the ratios become more
conducive towards increased response at superharmonic resonances and this will have a direct
impact on the reliability of wind turbine systems. Flutter onset also has increased possibilities.
Repair and maintenance for these mega structures will become increasing expensive. Accurate

98

Figure 5.9: Flutter speed as a function of blade length. Data obtained from [6]
modeling of external and internal loads of the system will enable us to set constraints on system
parameters to prevent wind turbines from having undesirable response and potential failures.

99

CHAPTER 6
DRIVETRAIN DYNAMICS

Vibrations in the blades are transferred to the gearbox via the rotor hub. Any drastic changes
in loading of the blades either due to wind gusts or due to increased vibration response due
to resonance phenomenon is thus transferred directly to the main shaft of the gearbox. This
necessitates accommodation of blade vibration loads while designing and certifying gearboxes for
wind turbines. The common practice is to model these elements for dynamic responses separately as
the critical frequencies are magnitudes apart. However, increasing failures in gearbox components
has ushered the industry towards gaining a complete understanding of all the dynamic loads the
drivetrain would experience during operation. In this chapter, we look at two modeling approaches
for wind turbine drivetrains.

6.1

Modeling of the Gearbox Reliability Collaborative (GRC) Gearbox

First, we look at a detailed drivetrain model built in RomaxWind, a commercial computer-aided
engineering software that is used for comprehensive analysis, optimization and certification of wind
turbine drivetrains. Some of the early validation work of the RomaxWind model of the NREL GRC
750 kW wind turbine drivetrain is included in this chapter as well.
The GRC is a program run by the NREL National Wind Technology Centre (NWTC) to
determine the reasons for premature failure that are common in gearboxes in wind turbines. The
GRC focuses on three areas of research: drivetrain modeling and analysis, full scale dynamometer
testing and field testing. A model of the GRC gearbox is shown in Figure 6.1, with some salient
features labeled. Figure 6.2 shows the schematic of the GRC gearbox with various speed stages and
components labeled. The gearbox has been modeled in the RomaxWind software, a virtual product
development and simulation environment for the design and analysis of wind turbine gearboxes,
bearings and drivetrains. The gearbox is a single stage planetary with two parallel helical gear
stages - a typical configuration for medium sized gearboxes in wind turbines. The gearbox is a
100

Figure 6.1: RomaxWind model of GRC gearbox
Figure 3: GRC gearbox internal components view

High Speed Stage
HS-ST

High Speed
Shaft
HS-SH

The figure below shows the nomenclature used to describe the locations of the bearings

Pinion

PLC

Planet

Pinion

Intermediate Speed
Shaft
IMS-SH

Annulus
Gear

Annulus HS-SH-A

Pinion

HS-SH-B HS-SH-C
Gear

Planet
Low speed Shaft
LS-SH

IMS-SH-A

Pinion

IMS-SH-B IMS-SH-C

PL-A PL-B
Sun

Gear

Low Speed Stage
LS-ST

PLC-A

Sun

PLC-B

Gear

LS-SH-A

Intermediate Speed Stage
IMS-ST

LS-SH-B LS-SH-C

Figure 5: Bearing nomenclature and location

Figure 4: GRC gearbox internal nomenclature and abbreviations

Figure 6.2: Schematic showing GRC gearbox layout and component nomenclatures
GRC Analysis Round Robin

6

The GRC Analysis Round Robin seeks to evaluate the different analytical tools and
approaches utilized in the development and design of horizontal axis wind turbine
drivetrains. The round robin utilizes the GRC drivetrain to provide its participants with a
comprehensive characterization of the gearbox and other components. From there, the
participants can develop models with different levels of complexity and fidelity.

speed increaser, with approximately 1:81 ratio, and is designed for high speed shaft speeds of 1200
A series of loading conditions, such as the calibration load case (CLC), are specific to this
rpm and 1800 rpm depending on the running configuration
of the two speed (4/6 pole) generator.
effort, though they are not a guideline according to the IEC standard. As the round robin
progresses, the loading conditions will increase in complexity ranging from fully static to
dynamic and transient events. This progressive approach assures that the discrepancies in
the analysis results are due to the different tools employed and not the result of an
analyst’s assumptions.

Two 750 kW gearboxes that have had some redesign from the original arrangement are considered for modeling and testing. They have been modified to allow a comprehensive instrumentation

101

8

Table 6.1: Validation of gearbox housing natural frequencies
Mode Number
1
2
3
4
5
5

Test Frequency (Hz)
317
341
359
394
415
445

Simulation Frequency (Hz)
338
359
387
422
464
494

% Error
6.21
5.01
7.23
6.64
10.56
9.92

package (over 150 channels). Measurements of tooth root strains on the ring gear, raceway strains
on the planet bearings, deflection of the planet carrier, temperatures near the bearings, input and
output torque and speed, amongst other operating values are recorded. We use these measurements
to validate the model of the drivetrain in RomaxWind. There are two aspects to having a reliable
model for predictive computation; (a) have a comprehensive model to capture all the necessary
physics to predict response of a system and (b) validate the model at component and system level to
ensure that the behavior of the model emulates the real drivetrain it represents. We used measurements from accelerometers of a vibration tests with impulse excitations to check if the frequency
and modal response of the gearbox housing matches those predicted by the simulation. Figure 6.3
shows a snapshot of the test set up and the corresponding simulation in RomaxWind. The gearbox
housing was instrumented with ∟300 accelerometers to capture the mode shape and we use the
measurements to qualitatively compare the mode shapes. From table 6.1, we can see that for the
first six modes the natural frequencies are at about 10% error. This is a sufficient level of validation
as the exact material properties, joints in the structure represented by finite element structures and
boundary conditions are not exactly captured.
The RomaxWind drivetrain model is built from the engineering drawings provided by NREL
GRC. Key features of the model include beam finite element representation of shafts, solid finite
element representation of gearbox housing, gear blanks, planet carrier and torque arms and 6DOF
spring connections for (elastomeric) trunnion mounts. The gears and bearings are modeled with
semi-analytical formulations that take account of important factors such as misalignment, area of
contact under load, micro-geometry, radial and axial clearances or preload and material properties.
102

Figure 6.3: Comparing the 1st natural frequency and modal response of the gearbox housing.
Figure shows the test set-up and test result (319 Hz) on the left and RomaxWind model and
simulation result (338 Hz) on the right
Static non-linear analysis can be performed for prescribed loading conditions and the global deflections are solved simultaneously with the contact mechanics for the gears and bearings. Static
analysis is considered sufficient as the inertia terms are orders of magnitude lower than the forces
due to the loads carried by the gearbox. Thus the effect of the whole system behavior on contact
elements is captured. We use some of these metric as measures to ensure our model is validated.
Designers typically use these models for achieving good alignment of the system under the
loads, calculating gear and bearing contact stress and life, optimizing the micro-geometry of the
gears for increased life and transmission error and predicting the gear vibration magnitudes, as well
as for many other purposes.
103

Figure 6.4: RomaxWind model of GRC gearbox installed in wind turbine with a load case
example: Input Torque (Mx) = 325 kNm, Off Axis Moment = -100 kNm, Rotor Weight = -1500 kN

6.2

Case Study - Influence of Carrier Bearing Clearances on Gear and
Bearing Alignment and Stresses

Figure 6.4 shows a model used for the purpose of simulating the GRC gearbox as installed
in the dyno. The load case shown with the figure includes a significant off-axis moment as well
as loading from rotor mass (wind turbine blades and hub). Large off-axis moments are common
in wind turbines, under conditions such as wind direction change, yaw error (putting the blades
askew to the wind), yawing of the machine and wind shear. Most dynamometer tests only load
the drivetrain with torque and it is important to understand that the loading and deflections of the
gearbox may be markedly different in the two installations; particularly for three point mounting
arrangement (main bearing and two torque arm supports).
For the given load case a non-linear static analysis has been performed and the contact and
stress in the planetary stage is calculated at a certain carrier rotation. Planet positions are chosen
at 0◦ , 120◦ and 240◦ and subsequent results focus on the planet/sun gear mesh and the planet and
carrier bearings are shown in table 6.2.
The load as well as clearance in the main bearings can misalign the carrier to the housing and

104

Table 6.2: Metrics on the Quality of Gear Contact
Planet Gear
Number
1
2
3

Maximum contact
stress (MPa)
1223
1146
1197

Misalignment
FβX (¾m)
173
103
57.2

Face load distribution
factor KHβ
1.72
1.25
1.35

ring gear and affect the gear contact (note the main bearing is a spherical roller and provides little
tilt stiffness).

6.2.1

Time-varying Misalignment and Stress

When the planet carrier rotates, the planets move with the carrier and subsequently the planet
load share and planetary gear misalignment vary with time. This results in increased vibration
and reduced life, increasing the probability of failure for the planetary gear set. By performing
the analysis at incremental rotations of the carrier the time varying behavior can be examined
and figures 6.5 and 6.6 chart the planet load share and sun/planet gear misalignment respectively.
The time varying misalignment and load sharing has significant effects on the planet bearings, a
component that is known to have high failure rates in wind turbines and are expensive to replace.

6.2.2

Optimizing the Design to Improve Contact, Stress and Life

A design analysis needs to be performed to understand the influences of misalignments on the
planet gears durability. If the misalignment can be reduced then the time varying contact stress and
face load distribution can be reduced. This will in turn decrease vibration and reduce fatigue on
the planet gears and planet bearings. For many wind turbine gearbox designs, the alignment of the
planetary stage is sensitive to the magnitude of the clearances in the planet carrier bearings. For
this model, with the gearbox installed in the dyno, a full factorial study has been performed where
the radial internal clearance of the upwind and downwind planet carrier bearings are varied and the
misalignment and face load distribution factors (FβX and Khβ [88]) are calculated for each planet.
The full factorial study (see figure 6.7) demonstrates how sensitive the misalignment of the sun and
105

Figure 6.5: Planet load share in kNm vs. carrier rotation in degrees

Figure 6.6: Sun/Planet gear misalignment in microns vs. carrier rotation in degrees

106

	
  

Figure 6.7: Results of full factorial study for upwind and downwind planet carrier bearing’s radial
internal clearance (shown as PLCA RIC and PLCB RIC in figure): (left) Sun/Planet
misalignment; (right) Sun/Planet load distribution.
planet is to these bearing settings and designers can use such models to optimize the parameters of
the bearings for improved performance and life.
The best clearance specification would be the subject of a full design review. However, a
brief example is provided here where the carrier clearances are set to an optimal value so that the
planet/sun misalignments are consistent with position. For this design and load case the optimum
settings are 50 Âľm and 250 Âľm on the upwind and downwind bearings respectively (indicated
by markers on Figure 6.7). With the optimized clearances a similar mean stress is achieved; the
maximum stress is reduced by approximately 5% and the time varying amplitude is reduced by a
factor of four. The gears are still misaligned, yet this misalignment is now consistent with rotation
(plots not shown) and can be corrected with lead slope adjustments.

6.2.3

Redesigning to Reduce Sensitivity

The improved carrier bearing clearance setting leads to lower vibrations levels in the components of
the gearbox. But, typical manufacturing tolerances on bore of the carrier and planetary alignments
are of the order of 100 microns. This would mean that even though we can specify clearances

107

Table 6.3: Metrics for the quality of gear contact for cases shown in figure 6.9
Case
(a)
(c)

Maximum load
per unit length (N/mm)
1131
921

Face load distribution
factor Khβ
1.72
1.11

while designing the gearbox, in practice, we may not achieve optimal settings. A change of a few
microns affects the load sharing significantly. This prompts the use of other design decisions to
reduce the sensitivity of the model to bearing clearances. For this system, we study the effect of
change in the bearing package along with a corresponding micro-geometry optimization of the
planet gear set. We replace the cylindrical roller bearing of the carrier with tapered roller bearings.
Micro-geometry corrections on the planet gear set are also made to correct gear misalignments
and improve the face load distribution. Figure 6.9 shows the improvement in the gear contact at
the face of a planet gear (which had extreme loading pattern in the original case) with a changed
bearing package and changes in the lead slope. Similar improvements are noticed for the contact
patches on other planet gears as well. Thus, an overall reduction in planet gear misalignments
and a corresponding improvement in face loading can be obtained. The improved design is quite
robust and small errors in planet position due to assembly don’t affect the load sharing patterns.
The proposed changes in clearance and bearing configuration also improves the loading patterns
on the bearing elements. We show the variation in maximum contact stress of the downwind planet
bearing for the two cases shown in figure 6.9. The variations in maximum stress per revolution is
reduced by about 70% and the maximum stress is also reduced by about 6%, both of which would
lead to lesser fatigue loading of the bearings.
Thus, we have used the RomaxWind model of the drivetrain, to assess loading of the drivetrain
under different input loads at the rotor hub, quantified sensitivities of the model to different
parameters in the gearbox and drivetrain design, and suggested methods to improve/optimize the
design of components in a wind turbine drivetrain to improve reliability. We continue to explore
the usage of these models for improving the predictive capability.

108

Figure 6.8: Sun/Planet load distribution for (top) original case; (bottom) case with tapered roller
bearing with +100 Âľm preload and after micro-geometry optimization

6.3

Prediction of Wind Turbine Gearbox Vibration

In this section, we use the RomaxWind model to predict the vibration levels at various location
on the gearbox housing. This was specifically used to improve our confidence in the model, but, at
the same time was a necessary step to understand the usage of this model to benchmark standard
failure modes and design gearboxes to meet noise and vibration targets throughout their warranty.
Impact tests and dyno running tests were performed on the gearbox mounted on the dyno at NREL.
They are necessary to tune housing material properties and bushing stiffness to attenuate any
undesirable frequency response. We seek to use this data to characterize the vibration spectrum of
a healthy gearbox and for fault diagnostics. There are numerous condition monitoring techniques
109

Figure 6.9: Downwind planet bearing maximum contact stress with carrier rotation for (top)
original case; (bottom) case with tapered roller bearing with +100 Âľm preload and after
micro-geometry optimization

110

that use vibration signatures measured realtime from operating machinery to diagnose or detect
potential issues or defects with machine components [89]. Advanced signal processing techniques
are employed to determine the root cause of certain frequencies in the spectrum and trace it back
to component malfunction.
In figure 6.10, we compare measurements for vibration response of sensors at various locations
between test and RomaxWind model. It is clear that the trend for sensitivity of sensor position
towards the response amplitude correlates well, i.e., we are able to accurately identify the sensor
that detects the maximum amplitude for the two cases shown. The location of the sensor is crucial
to the response it records due to structural excitations. The precise location should be known and
a finite element node is the exact same location should be “tagged" in the RomaxWind model for
accurate analysis.
Capability to predict vibration response accurately in the RomaxWind model serves as an
important assessment tool in the design phase of the gearbox development process. This model
will be used to optimize the sensor location for diagnostics and also to design the gearbox for
low vibration. For design improvements, requirements such as limitations on housing deflection at
certain locations can be checked. While for condition monitoring, depending on the whole system’s
dynamic response, the choice of sensor location for vibration based condition monitoring can be
decided to get the best output. However, vibration amplitudes are difficult to predict from simulation
tools across different frequencies and care should be taken to model the structural flexibility of the
entire wind turbine system as accurately as possible to get a representative frequency response.

6.4

Redesigning the Wind Turbine Gearbox: For Improved Reliability

The NREL GRC gearbox phase 3 redesign was performed to assess an existing design of
the gearbox and suggest candidate design changes from the conceptual review amongst other
requirements. The specific task performed as a part of this dissertation, performed in collaboration
with Romax Technology Inc., was to provide an exploration of the effects of the design details on
the loads and deflections within the gearbox. This involved use of the RomaxWind model of the

111

Figure 6.10: Normalized amplitude of vibration at different sensor location for High speed shaft
(HSS) and Intermediate speed shaft (IMS) gear excitations
gearbox, which has been validated against many measurement on the gearboxes from phase 1 and
2. Candidate design changes arising from the concept design review relating to the planetary stage
were evaluated. Some specific recommendations made are discussed in this section.

6.4.1

Tapered Roller Bearing for Carrier and Planet

Tapered roller bearing offer more load-carrying and load-distributing ability. The bearing suggested
by NREL GRC supplier was assessed in “X" (force direction on pins make an X) and “O" (force
direction on pins make an O) arrangements and even though the “O" offered more stiffness,
112

Figure 6.11: RomaxWind model with carrier and planet gear tapered roller bearing in “X"
(marked at location of bearing by rectangles) and “O" (marked at location of bearings by circle)
configurations
it is not recommended due to difficulty in assembly for the carrier (figure 6.11). The driving
factor for choosing a taper roller bearing is to eliminate clearance and increase stiffness, hence
improving planetary stage gear mesh load distributions. As the preload is increased the ring
gear mesh misalignments become closer together (see figure 6.12). As clearance is increased the
misalignments diverge. If the misalignments are all close together then it is possible to correct with
lead slope modification. A set of misalignments that are different at each planet position indicates
some significant misalignment between the carrier and the housing or the sun.
The carrier cylindrical roller bearings should be replaced with taper roller bearings. This
decreases the misalignment between the carrier and the housing. The planet cylindrical roller
bearings should be replaced with integral gear/outer race taper roller bearings. This reduces
the misalignment due to the moment from the axial component of the mesh forces on the planet,
improves the load share between the upwind and downwind bearing rows and removes the possibility
of fretting between the outer race and the gear blank. The selection of a bearing with significantly
greater load capacity increases the system life calculation significantly as the planet bearings are
113

Figure 6.12: Mesh misalignment vs. preload for initial tapered roller bearing selection
the lowest life component.

6.4.2

Interference Fit Planet Pins

If the planet pin are designed to have a tight fit, then the carrier will also have a stiffening effect.
A comparison between the baseline, and models with carrier tapered roller bearings, ring gear
stiffening, carrier stiffening, planet pin interference fits is shown in figure 6.13. It can be seen
that increasing the ring gear stiffness has not had a significant effect on mesh misalignment. The
increase in carrier stiffness has had a moderate effect because this reduces the misalignment due to
torsional wind-up of the carrier. However, the most dramatic effect is found by interference fitting
the pins into the carrier. This has a greater stiffening effect in the carrier and consequently would
be easier to modify in the design of the gearbox.

114

Figure 6.13: Mesh misalignment for increasing ring and carrier stiffness with interference fit pins
for planetary gear contacts
6.4.3

Decreased Spline Crowing

The sun spline showed damage around the center contact region once the above design changes
were implemented. Reducing the crowning on the external spline would increase the contact area
and should help negate the damage in the center contact region. A potential knock-on effect is an
increase in the tilt stiffness of the spline and hence a worsening of planetary gear misalignment
and load share. However, in the current model for the load cases studied, the misalignments on the
planetary gear set were not affected due to changes in spline crown. We simulate changes in spline
crown from the baseline value of 173 µm to 50%, 75% and 87.5% of it’s value. Figure 6.14 shows
the deflection distribution of the spline teeth. The purpose of having the flexible spline joint here
is to allow the sun to float in the planetary stage. Decreasing the spline crown will increase the tilt
stiffness of the spline and reduce the amount of sun float. However, this is deemed to be acceptable
tradeoff for the reduction in peak center load on the spline.

115

Figure 6.14: Deflection distribution on tooth on horizontal axis for decreasing crown of (top left)
173Âľm, (top right) 86.5Âľm, (bottom left) 43.25Âľm, (bottom right) 21.625Âľm
6.4.4

Other Changes and Recommendations

There are other more detailed design changes that will have an impact on the reliability of the
gearbox design. In this section we explored a few design modification in an existing gearbox with
constrains on overall geometry and parts being used. However, future work in gearbox redesign for
the NREL GRC project or otherwise can explore other design changes listed below (though some
recommendations are specific to the layout of the GRC gearbox) and its implications on loading
and misalignments.
• The ability to modify the carrier taper roller bearing preload
• The ability to modify the planet taper roller bearing preload
• A method for adjusting the stiffness of the ring gear
• An optimization of gear micro-geometry
• An optimization of bearing preload settings
• The ability to change out the high speed pinion
• An improved design of the high speed shaft bearing system
116

• An improved design of the hollow shaft taper roller bearing support and lock-nut system
• An improved design of the oil feed rings
• A review of the assembly and disassembly procedure and proposed methods for improvements
• Improvements to the gearbox lubrication system and sealing, and the including of a instrumentation methodology to allow measurement of lubrication effectiveness

6.5

Torsional Model of Wind Turbine Drivetrain

Since RomaxWind is geared towards detailed component static analysis, in order to perform
transient analysis a simple matlab model was built in collaboration with Romax Technology Inc. in
order to capture the dynamic interaction between the wind turbine blades, drivetrain and generator.
A schematic of the torsional model is shown in figure 6.15. The model has been used to characterize
and study different drivetrain configurations. Some preliminary studies with this model involved
comparing the loads on drivetrain components of a medium speed drivetrain (new concept) vs.
a high speed drivetrain (existing concept). The moments transmitted and the associate drivetrain
deflections in the rotational axis are also captured.

6.6

Extension to Loading from Blade Vibration

In the previous sections, it has been shown that the loading in the gearbox and pressure at the
roller elements of bearings are highly sensitive to the design of the system. Moreover, the deflection
during operation of the machinery causes variations in the loads experienced and the quality of
contact in the gear tooth. It is hence clear that any blade vibration beyond designed levels will cause
changes in the loading experienced in the gearbox via coupling through the rotor hub. However, the
increased loads at superharmonic resonances cannot be easily implemented in RomaxWind model.
For starters, we would need the detailed geometry of multiple gearboxes to assess the impact of
increased vibration loads for different sizes of blades. As this is not feasible within the scope of

117

Figure 6.15: 1-D torsional model of a simple wind turbine gearbox, courtesy of Romax
Technology
this thesis, the torsional drivetrain model developed can be used for such studies. The aerodynamic
loads computed using the BEM theory in chapter 5 can also be used as inputs to this model.
The solution of the IPDE that are used to model the blade vibration would be considered.
We would seek solution to the moments that are transferred at the root of the blade. Existing
commercial codes such as FAST are already able to do this. Our analysis would account for
significant nonlinearity in the blades (which tend to be more longer and flexible) and associated
harmonic resonances. The dominant blade frequency and the associated harmonics typically occur
at much lower frequencies than the gear and bearing related frequencies. From torsional models
(as shown in figure 6.15), simulations that relate to the influence of gear mesh stiffness, component
inertias and mount damping to overall system behavior can be studied. Additionally, we can include
torsional fluctuations at the rotor hub that are of the order of blade non linear dynamics and assess
its impact.
118

CHAPTER 7
CONCLUSION AND FUTURE WORK

7.1

Conclusion

This work considered the cyclic loads on a wind turbine blade and how they affect blade
vibration and drivetrain dynamics. The work first provided integral-partial differential equations
and reduced ordinary differential equations of motion that model the in-plane and out-of-plane
dynamics of a wind turbine blade. Nonlinear effects (due to curvature and foreshortening) and
cyclic loading (due to gravity and aerodynamics) have been considered.
The single mode reduction led to the parametrically and directly forced nonlinear ODE equation
(2.8). This motivated the investigation of a forced nonlinear Mathieu equation as a paradigm in
dynamical systems research. Furthermore, we found that when typical wind turbine parameters are
included, equation (2.8) may not have the convenience of a small parameter. So as a preliminary
study, we examined the small parameter nonlinear forced Mathieu equation (in which the direct and
parametric excitations occur at the same frequency, out of phase by 90 degrees, and which includes
a cubic nonlinearity), with the expectation that the phenomena uncovered may carry over to wind
turbines. Perturbation analysis of equation (3.1) revealed the presence of multiple super- and subharmonic resonances. The superharmonic resonance of order 1/2 persists for the linear Mathieu
equation with direct forcing. Some resonance and instability phenomena were uncovered using
simulations. Second-order perturbation methods provided the analytical framework to understand
system behavior better. Primary resonance and multiple superharmonic and subharmonic resonance
cases for system with hard and weak excitations, direct forcing and nonlinearity were studied. The
presence of superharmonic resonances may be significant for wind turbines, which are designed to
operate below the lowest natural frequency. Resonances can be amplified beyond linear predictions
if there is simultaneous superharmonic and primary excitation as well.
We further extended the analytical framework to real-world example wind turbine systems.
119

Turbines with blade diameters ranging from 9 m to 100 m were simulated. The supporting
aerodynamic theories to assess aero loads were also benchmarked and developed. The phenomenon
of harmonic resonance and instability is likely to occur for larger wind turbines and this is shown
via a systematic analysis.
To understand the effect of unknown dynamic loading conditions on a gearbox, we have analyzed
the NREL GRC gearbox model. As a case study, the sensitivity to design and tolerances was
performed. Detailed models within the RomaxWind software tool have been developed with focus
on determining the quality of the function of the planetary gear stages, and predictions of planet
load distributions have been validated against the test data. Key design drivers are discussed, such
as the quality of alignment at the gears and bearings and the loads and stresses seen on these
components. It is illustrated that small clearances, such as in the carrier bearings, can have a large
effect on the performance of the design and a study shows how to identify and reduce time varying
misalignment and contact stresses which leads to lower vibration, lower fatigue and a more reliable
product. This model will be used to demonstrate the effect of superharmonic and subharmonic
resonances on gearbox dynamics and reliability of the wind turbine system.
Improving the reliability of wind turbines is important for the industry and for reducing the
cost of energy. The NREL GRC program is one example of industry, university and government
research labs working together to improve the understanding of the gearbox and exploring and
disseminating the reasons for premature failure. Other such projects and government initiatives
around the globe point to a pressing need towards understanding issues that are critical to the proper
functioning and reliability of wind turbine structures.

120

7.2

Future Work

The foundations laid in this thesis can be extended in multiple directions. Some of them listed
below:
• Develop model for coupled in-plane, out-of-plane and twist motion of the blade.
• Second order perturbation analysis for complete nonlinear partial differential equation of
wind turbine blades.
• Use the model of the beam to validate stable/unstable operating zones by experimenting with
a cantilever beam attached to a shaker. The attachment will allow for inclined orientation to
mimic parametric forcing.
• Use semi-empirical aerodynamic theory to get a more representative expression for forcing
and incorporate it into the calculations.
• Use calculated moments at the root of the blade as input to the NREL GRC gearbox model
in Romax and arrive at more realistic loading scenarios within the gearbox.

121

BIBLIOGRAPHY

122

BIBLIOGRAPHY

[1] Global Wind Report Annual Market Update, Global Wind Energy Council. Online, 2013.
[2] U.S. Wind Industry Annual Market Report, American Wind Energy Association. Online,
2013.
[3] Anders Ahlstrom. Aeroelastic simulation of wind turbine dynamics. PhD thesis, Royal
Institute of Technology, Sweden, 2005.
[4] Todd Griffith, Thomas Ashwill, and Brian Resor. Large offshore rotor development: Design
and analysis of the Sandia 100-m wind turbine blade. In 53rd AIAA Structures, Structural
Dynamics and Materials Conference, volume 1499. AIAA, April 2012.
[5] H. Link, W. LaCava, J. van Dam, B. McNiff, S. Sheng, R. Wallen, M. McDade, S. Lambert,
S. Butterfield, and F. Oyague. Gearbox reliability collaborative project report: Findings from
phase 1 and phase 2 testing. Technical Report NREL/TP-5000-51885, NREL, June 2011.
[6] http://energy.sandia.gov/energy/renewable-energy/wind-power/offshore-wind/offshorewind-sandia-large-rotor-development/.
[7] 20% wind energy by 2030 – increasing wind energy’s contribution to the us energy supply.
US Department of Energy, DOE/GO-102008-2567, July 2008.
[8] Wind turbine supply train strategies 2009-2020. Emerging Energy Research, July 2009.
[9] Landwirtschaftskammer Schleswig-Holstein (LWK), Schleswig-Holstein, Germany.
[10] P.J. Tavner, G.J.W. van Bussel, and E. Koutoulakos. Reliability of different wind turbine
concepts with relevance to offshore application. Belgium, March - April 2008. European
Wind Energy Conference (EWEC).
[11] F. Oyague, C.P. Butterfield, and S. Shuangwen. Gearbox reliability collaboration analysis
round robin. Technical Report CP-500-45325, May 2009.
[12] Germanischer Lloyd WindEnergie. Guideline for the Certification of Wind Turbines, 2003
edition, (with 2004 supplement).
[13] P Marzocca, L Librescu, and W A Silva. Aeroelastic response and flutter of swept aircraft
wings. AIAA Journal, 40(5):801–812, 2002.
[14] P J Attar, E H Dowell, and D M Tang. A theoretical and experimental investigation of the
efects of a steady angle of attack on the nonlinear flutter angle of a delta wing plate mode.
Journal of Fluids and Structures, 17(2):243–259, 2003.
[15] R A Ibrahim. Friction induced vibration, chatter, squeal, and chaos. part II: Dynamics and
modeling. Applied Mechanics Reviews, 47:227, 1994.

123

[16] Brian F Feeny. The nonlinear dynamics of stick-slip oscillators. Dynamics and Friction, 1996.
[17] Gunjit Bir and Jason Jonkman. Aeroelastic instabilities of large offshore and onshore wind
turbines. Journal of Physics: Conference Series, 75(012069), 2007.
[18] Ali H. Nayfeh. Nonlinear interactions: analytical, computational and experimental methods.
Wiley, New York, 2000.
[19] M Golubitsky and D G Schaeffer. Singularities and groups in bifurcation theory. SpringerVerlag, New York, 1985.
[20] C Pierre and P D Cha. Strong mode localization in nearly periodic disordered structures.
AIAA Journal, 27(2):227–241, 1989.
[21] C Pierre and D V Murthy. Aeroelastic modal characteristics of mistunes blade assemblies Mode localization and loss of eigenstructure. AIAA Journal, 30(10):2483–2496, 1992.
[22] D W Lobitz. Parameter sensitivities affecting the flutter speed of a MW-sized blade. Journal
of Solar Energy Engineering - Transaction of the ASME, 127(4):538–543, 2005.
[23] D W Lobitz and P S Veers. Aeroelastic behavior of twist-coupled HAWT blades. In Proc. of
the 1998 ASME/AIAA Wind Energy Symposium, pages 75–83, Reno, 1998.
[24] D W Lobitz. Aeroelastic stability prediction of a MW-sized blade. Wind Energy, 7(3):211–
224, 2004.
[25] M H Hansen. Stability analysis of three-bladed turbines using an eigenvalue approach. In
Proc. of the 2004 ASME/AIAA Wind Energy Symposium, pages 192–202, Reno, 2004.
[26] E H Dowell. A modern course in aeroelasticity. Kluwer Academic Publishers, Boston, 4th
revised and enlarged edition, 2005.
[27] G Bir and J Jonkman. Aeroelastic instabilities of large offshore and onshore wind turbines.
In TEAWE Torque from Wind Conference, number NREL/CP-500-41804, pages 28–31, Denmark, August 2007. Technical University of Denmark.
[28] S J Schreck. Rotationally augmented flow structures and time varying loads on turbine blades.
Number NREL/CP-500-40982 in 45th AIAA Aerospace Sciences Meeting and Exhibit, pages
8–11, Reno, Nevada, January 2007.
[29] J. Wendell. Simplified aeroelastic modeling of horizontal-axis wind turbines. Technical
Report DOE/NASA/3303-3; NASA-CR-168109, Massachusetts Inst. of Tech., Cambridge,
USA, September 1982.
[30] B S Kallesøe. Equations of motion for a rotor blade, including gravity and pitch action. Wind
Energy, 10(3):209–230, February 2007.
[31] D.I. Caruntu. On nonlinear forced response of nonuniform blades. In Proceeding of the 2008
ASME Dynamic Systems and Control Conference, number DSCC2008-2157. ASME, October
2008.
124

[32] Inderjit Chopra and J Dugundji. Nonlinear dynamic response of a wind turbine blade. Journal
of Sound and Vibration, 63(2):265–286, 1979.
[33] Jason Jonkman. NWTC design codes (FAST), Novermber 2010.
[34] Jason Jonkman. Modeling of the UAE wind turbine for refinement of FAST AD. Technical
Report TP-500-34755, NREL, Golden, CO, December 2003.
[35] Guido Weinzierl. A BEM based simulation-tool for wind turbine blades with active flow
control elements. Master’s thesis, Technische Universität Berlin, April 2011.
[36] Tonio Sant. Improving BEM-based aerodynamic models in wind turbine design codes. PhD
thesis, Delft University, January 2007.
[37] Y R Wang and D A Peters. The lifting rotor inflow mode shapes and blade flapping vibration
system eigen-analysis. Computer Methods in Applied Mechanics and Engineering, 134(12):91–105, 1996.
[38] S R Bhat and R Ganguli. Validation of comprehensive helicopter aeroelastic analysis with
experimental data. Defence Science Journal, 54(4):419–427, 2004.
[39] O Rand and S M Barkai. A refined nonlinear analysis of pre-twisted composite blades.
Composite structures, 39(1-2):39–54, 1997.
[40] E C Smith and Inderjit Chopra. Aeroelastic response, loads and stability of a composite rotor
in flight forward. AIAA Journal, 31(7):1265–1273, 1993.
[41] D.H. Hodges and E.H. Dowell. Nonlinear equations of motion for elastic bending and torsion
of twisted nonuniform rotor blades. Technical Note D-7818, NASA, 1974.
[42] Theodore Theodorsen. General theory of aerodynamic instability and the mechanism of
flutter. Technical report, National Advisory Committee for Aeronautics, 1935.
[43] Grant Ingram. Wind turbine blade analysis using the blade element momentum method.
Durham University, December 2005.
[44] Juan Mendez and David Greiner. Wind blade chord and twist angle optimization by using genetic algorithms. Institute of Intelligent Systems and Numerical Applications in Engineering.
[45] Jack Wasey. Optimizing wind turbine blade shape and pitch control strategy using blade
element theory, Athens 2004.
[46] R Lanzafame and M Messina. Fluid dynamics wind turbine design: Critical analysis, optimization and application of bem theory. Renewable Energy, 32(14):2291–2305, November
2007.
[47] Rachid Younsi, Ismail El-Batanony, Jeur-Bernard Tritsch, Hassan Naji, and Bernard Landjerit.
Dynamic study of a wind turbine blade with horizontal axis. Eur. J. Mech. A/Solids, 20:241–
252, September 2001.

125

[48] J Gordon Leishman. Challenges in modeling the unsteady aerodynamics of wind turbines.
40th AIAA Aerospace Sciences Meeting and 21st ASME Wind Energy Symposium. AIAA
2002-0037.
[49] A Baumgart. A mathematical model for wind turbine blades. Journal of Sound and Vibration,
251(1):1–12, March 2002.
[50] Henrik Broen Pedersen and Ole Jesper Dahl Kristensen. Applied modal analysis of wind
turbine blades. ISBN 87-550-3170-6 (Internet), Risø National Laboratory, Denmark, February
2003.
[51] W A A M Bierbooms. A comparison between unsteady aerodynamic models. Journal of
Wind Engineering and Industrial Aerodynamics, 39:23–33, 1992.
[52] J.W. Larsen, S.R.K. Nielsena, and S. Krenkb. Dynamic stall model for wind turbine airfoils.
Journal of Fluids and Structures, 23:959–982, May 2007.
[53] A C Hansen and C P Butterfield. Aerodynamics of horizontal axis wind turbines. Annuak
Reviews of Fluid Mechanics, 25:115–149, 1993.
[54] Tsuyoshi Inoue, Yukiko Ishida, and Takashi Kiyohara. Nonlinear vibration analysis of the
wind turbine blade (occurrence of the superharmonic resonance in the out of plane vibration
of the elastic blade). Journal of Vibration and Acoustics, (accepted).
[55] Yukio Ishida, Tsuyoshi Inoue, and Kohei Nakamura. Vibration of a wind turbine blade
(theoretical analysis and experiment using a single rigid blade model). Journal of Environment
and Engineering, 4(2):443–454, July 2009.
[56] Pramod Malatkar. Nonlinear Vibrations of Cantilever Beams and Plates. PhD thesis, Virginia
Polytechnic Institute and State University, Blacksburg, Virginia, July 2003.
[57] Mustafa Sabuncu and Kaan Evran. The dynamic stability of a rotating pre-twisted asymmetric
cross-section blade subjected to lateral parametric excitation. Finite Elements in Analysis and
Design, 42:1113–1122, 2006.
[58] C Venkatesan and V.T. Natarajan. Non-linear flapping vibrations of rotating blades. Journal
of Sound and Vibration, 84(4):549–556, 1982.
[59] H. H. Yoo and S. H. Shin. Vibration analysis of rotating cantilever beams. Journal of Sound
and Vibration, 212(5):807 – 828, 1998.
[60] Y N Al-Nassar, M Kalyon, M Pakdemirli, and B O Al-Bedoor. Stability analysis of rotating
blade vibration due to torsional excitation. Journal of Vibration and Control, 12(9-10):1379
– 1391, 2007.
[61] Jeffrey Rhoads and Stever Shaw. The impact of nonlinearity on degenerate parametric amplifiers. Applied Physics Letters, 96(23), June 2010.
[62] Jeffrey Rhoads, Nicholas Miller, Steve Shaw, and Brian Feeny. Mechanical domain parametric
amplification. Journal of Vibration and Acoustics, 130(6):061006(7 pages), December 2008.
126

[63] Manoj Pandey, Richard Rand, and Alan T. Zehnder. Frequency locking in a forced Mathieuvan der Pol-Duffing system. Nonlinear Dynamics, 54(1-2):3–12, February 2007.
[64] L.A. Month and R.H Rand. Bifurcation of 4-1 subharmonics in the non-linear Mathieu
equation. Mechanics Research Communication, 9(4):233–240, 1982.
[65] W I Newman, Richard Rand, and A L Newman. Dynamics of a nonlinear parametrically
excited partial differential equation. Chaos, 9(1):242–253, March 1999.
[66] Leslie Ng and Richard Rand. Bifurcations in a Mathieu equation with cubic nonlinearities.
Chaos Solitions and Fractals, 14(2):173–181, August 2002.
[67] M. Belhaq and M. Houssni. Quasi-periodic oscillations, chaos and suppression of chaos in
a nonlinear oscillator driven by parametric and external excitations. Nonlinear Dynamics,
18(1):1–24, June 1999.
[68] F Veerman and F Verhulst. Quasiperiodic phenomena in the van der Pol-Mathieu equation.
Journal of Sound and Vibration, 326(1-2):314–320, September 2009.
[69] D.K. Arrowsmith and R.J. MondragĂłn. Stability region control for a parametrically forced
Mathieu equation. Mecannica, 34:401–410, December 1999.
[70] Amol Marathe and Anindya Chatterjee. Asymmetric Mathieu equations. Proceedings of The
Royal Society A, (462):1643–1659, February 2006.
[71] N.W. McLachlan. Theory and Application of Mathieu Functions. Dover Publications, New
York, 1964.
[72] Leonard Meirovitch. Principles and Techniques of Vibration. Upper Saddle River, New Jersey,
1997.
[73] J.J. Stoker. Nonlinear Vibrations. Macmillan, New York, 1967.
[74] Ali H. Nayfeh and Dean T. Mook. Nonlinear Oscillations. Wiley Interscience Publications.
John Wiley and Sons, New York, 1979.
[75] Richard Rand.
Lecture notes on nonlinear vibration; available
http://audiophile.tam.cornell.edu/randdocs/nlvibe52.pdf. Ithaca, NY, 2005.

online

at

[76] Joris PEETERS. Simulation of dynamic drive train loads in a wind turbine. PhD thesis,
KATHOLIEKE UNIVERSITEIT LEUVEN, June 2006.
[77] Robert D. Blevins. Formulas for natural frequency and mode shape. Krieger Pub Co, 1979.
[78] Venkatanarayanan Ramakrishnan and Brian F Feeny. In-plane nonlinear dynamics of wind
turbine blades. In ASME International Design Engineering Technical Conferences, 23th
Biennial Conference on Vibration and Noise, pages no. DETC2011–48219, on CD–ROM,
Washington, D.C., 2011.
[79] M. H. Holmes. Introduction to Perturbation Methods. Springer-Verlag, 1998.
127

[80] James A. Murdock. Perturbations: Theory and Methods. A Wiley-Interscience publication,
1991.
[81] Ali H. Nayfeh. Perturbation Methods in Nonlinear Dynamics. Lecture Notes in Physics (247),
pp. 238-314. Springer-Verlag, Berlin, 1986.
[82] M. Sayed and Y. S. Hamed. Stability and response of a nonlinear coupled pitch-roll shi model
under parametric and harmonic exitations. Nonlinear Dynamics, 64(3):207–220, May 2011.
[83] L. A. Romero, J. R. Torczynski, and A. M. Kraynik. A scaling law near the primary resonance
of the quasiperiodic mathieu equation. Nonlinear Dynamics, 64(4):395–408, June 2011.
[84] Howard J Susskind. A stability analysis of the Mathieu equation with order-one damping.
Master’s thesis, Cornell University, August 1991.
[85] Randolph S Zounes and Richard H Rand. Subharmonic resonance in the non-linear Mathieu
equation. International Journal of Non-Linear Mechanics, 37:43–73, 2002.
[86] Venkatanarayanan Ramakrishnan and Brian F Feeny. Resonances of the forced Mathieu
equation with reference to wind turbine blades. Journal of Vibration and Acoustics, 134(6),
December 2012.
[87] Matthew Allen, Michael Sracic, Shashank Chauhan, and Motren Hartvig Hansen. Output-only
modal analysis of linear time periodic systems with appications to wind turbine simulation
data. In IMAC XXVIII, Jacksonville, FL, 2010.
[88] International Standard. Calculation of Load Capacity of Spur and Helical Gears – Part 1:
Basic Principles, Introduction and General Influence Factors, ISO 6336-1.
[89] Ian Howard. A review of rolling element bearing vibration “Detection, Diagnosis and Prognosis". Technical Report NAV 92/119, Aeronautical and Maritime Research Laboratory, October
1994.

128