VIBRATIONAL ANALYSIS OF VERTICAL AXIS WIND TURBINE BLADES
By
Onur Kapucu

A THESIS
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Mechanical Engineering - Master of Science
2014

ABSTRACT
VIBRATIONAL ANALYSIS OF VERTICAL AXIS WIND TURBINE BLADES
By
Onur Kapucu
The goal of this research is to derive a vibration model for a vertical axis wind turbine
blade. This model accommodates the affects of varying relative flow angle caused by rotating
the blade in the flow field, uses a simple aerodynamic model that assumes constant wind
speed and constant rotation rate, and neglects the disturbance of wind due to upstream blade
or post. The blade is modeled as elastic Euler-Bernoulli beam under transverse bending
and twist deflections. Kinetic and potential energy equations for a rotating blade under
deflections are obtained, expressed in terms of assumed modal coordinates and then plugged
into Lagrangian equations where the non-conservative forces are the lift and drag forces and
moments. An aeroelastic model for lift and drag forces, approximated with third degree
polynomials, on the blade are obtained assuming an airfoil under variable angle of attack
and airflow magnitudes. A simplified quasi-static airfoil theory is used, in which the lift and
drag coefficients are not dependent on the history of the changing angle of attack. Linear
terms on the resulting equations of motion will be used to conduct a numerical analysis and
simulation, where numeric specifications are modified from the Sandia-17m Darrieus wind
turbine by Sandia Laboratories.

Copyright by
ONUR KAPUCU
2014

I would like to dedicate this work to my mother Gonca Kapucu and my father Ali
Kapucu.

iv

ACKNOWLEDGMENTS

I would like to thank my advisor and my master Professor Brian Feeny for his endless
technical and mental support and sharing his wisdom generously to influence me to go
further. Without his dedication this research would have never existed.
I would also like to thank dear Katy Luchini Colbry, for always being there to listen and
help me.

v

TABLE OF CONTENTS

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

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

Chapter 1 Introduction . .
1.1 Background . . . . . .
1.1.1 Efficiency . . .
1.1.2 Blade Dynamics
1.1.3 Installation . .
1.2 Motivation . . . . . . .
1.3 Objective . . . . . . .
1.4 Contribution . . . . . .
1.5 Thesis Outline . . . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

. . . . . . . . . . . . . . . . . viii

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

1
1
4
5
5
6
8
8
8

Chapter 2 Formulation and Modeling
2.1 Elastic Blade Model for Bending . .
2.2 Reduced Order Modeling . . . . . .
2.3 Aeroelastic Modeling . . . . . . . .
2.4 Linearization . . . . . . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

10
10
14
17
25

Chapter 3 Numerical Analysis and Simulation .
3.1 Definition of the Giromill/H-Rotor Blade . . . .
3.2 Modal Frequencies . . . . . . . . . . . . . . . .
3.3 Blade with Quasistatic Aeroelastic Airfoil Model

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

27
27
31
32

Chapter 4 Conclusion and Comments
4.1 Summary of Results . . . . . . . .
4.2 Contribution . . . . . . . . . . . . .
4.3 Future Work . . . . . . . . . . . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

35
35
36
36

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

38

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

72

APPENDIX

.
.
.
.
.
.
.
.
.

vii

.
.
.
.

.
.
.
.

.
.
.
.

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

BIBLIOGRAPHY

vi

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

LIST OF TABLES

Table 3.1

Selected aluminium alloy’s material properties . . . . . . . . . . . .

30

Table 3.2

Specifications of the model created . . . . . . . . . . . . . . . . . . .

31

vii

LIST OF FIGURES

Figure 1.1

A sample HAWT sketch on left[2] and VAWT sketch on right [3] . .

2

Figure 1.2

A sample sketch for a Savonius WT on the left, Darrieus WT in the
middle and H-rotor configuration on the right [5] . . . . . . . . . . .

3

Comparison of aerodynamic efficiencies of different wind turbine types
[10] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

Basic VAWT Configurations: (a) Full Darrieus, (b)“H” , (c) “V”, (d)
“∆”, (e) “Diamond” and (f) “Giromill”[6] . . . . . . . . . . . . . . .

7

A sketch of blade and rotor. Undeformed on left and deformed on
right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

Top view representation of transverse deflection on top and twist on
bottom. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

Figure 2.3

Cross section of VAWT showing the coordinate axes . . . . . . . . .

18

Figure 2.4

Twist angle ψ and relative wind angle φ for the case when the pretwist
β=0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

Figure 2.5

Lift and drag force per unit length directions from top view . . . . .

20

Figure 2.6

NACA0012 CL and CD values with respect to angle of attack [21] .

21

Figure 3.1

Sandia 17m front view

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

28

Figure 3.2

Sandia 17m blade cross section showing aluminum spar and spline,
Nomex core and fiberglass skin. [22] . . . . . . . . . . . . . . . . . .

29

Figure 3.3

Blade dimensions for calculations

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

29

Figure 3.4

Variation of wnq with Ω on left and blade length a on right . . . . .

32

Figure 3.5

Variation of wnq and wnb with a when Ω = 0 . . . . . . . . . . . . .

33

Figure 3.6

Rotor dimensions used in calculations . . . . . . . . . . . . . . . . .

34

Figure 1.3

Figure 1.4

Figure 2.1

Figure 2.2

viii

Chapter 1
Introduction

1.1

Background

Our society’s increasing demand on energy has reached a critical level where the future
projections show us our fossil fuel sources will not be able to supply the demand in near
future. Being only dependent on several energy sources makes humankind vulnerable to energy shortage. Moreover, fossil fuel dependence in energy production unnaturally empowers
countries with rich fossil sources and increases the tension of international politics.
These facts show us the importance of the search for renewable energy sources. Being
aware of the aforementioned facts, most countries are now challenging a productive competition to develop and master the energy sources of the future. Facing the exhaustion of fossil
sources, we are spending more effort to exploit renewable sources such as wind, solar energy,
geothermal energy.
Among the renewable energy alternatives, wind power is a popular player. Energy Information Administration’s Annual Energy Outlook 2014[1] reference case projects 76% increase
in electricity production from wind power by 2040, where this increase is more than 300% in
the extended policies case. There exist several ways to absorb wind power and convert it into
a useful energy form. Wind turbines are the most popular setups to accomplish this purpose.
They basically convert the power of the wind into electricity via rotation in the generator.
Wind turbines can be distinguished into two main classes depending on the orientation of
1

the rotor: horizontal axis wind turbines (HAWTs) and vertical axis wind turbines (VAWTs).
In this notation “axis” represents the axis of rotation of rotor. Axes of rotation of VAWT
rotors are perpendicular to the ground, where they are parallel to ground in HAWT blades.
See Figure 1.1.

Figure 1.1 A sample HAWT sketch on left[2] and VAWT sketch on right [3]

Horizontal axis wind turbines dominate the current market; however, this might simply
be a coincidence [4]. The most common HAWT configuration has its blades facing the wind
to generate a rotation as a result of the lift forces on the blades. Being widely used, they
can be seen in various sizes from stand-alone installations to offshore fleets. As the energy
demand has increased over the time, HAWTs have become larger to produce more energy
from a single turbine. As of 2014, tallest commercial HAWT (MHI VESTAS V164) reaches
a tower height of 140 meters with a 220 meter tip height and 82 meter blade length and can
generate 8MW energy. Yet, as HAWT sizes have increased, problems with installation and
durability also have arisen. This might be a limiting factor for the future of the HAWTs [5].
The second most common wind turbine configuration, the VAWT, has its blades rotating
around a vertical axis. We can group VAWTs mainly into two: ones using lift-momentum
airfoils (i.e. Darrieus wind turbines) and ones using drag cups (Savonius wind turbine)(See
Figure 1.2). Among these there exist numerous configurations; such as, the Giromill/H2

rotor configurations, which are common types of Darrieus Turbine, due their simplicity.
Giromill/H-rotor style VAWTs consist of simply designed airfoils attached parallel to the
rotor via several struts or a single strut.

Figure 1.2 A sample sketch for a Savonius WT on the left, Darrieus WT in the middle and
H-rotor configuration on the right [5]

Our most significant experience with commercial VAWT installations dates back to the
early 80’s [5, 6]. Manufacturer FloWind initiated a large operation to install Darrieus type
VAWTs in the California Altamont and Tehachipi Wind Farms. The company installed more
than 500 turbines and reached a maximum total production of 105,000 megawatt/hour in
1987. However, shortly after this, the turbine fleet began collapsing due to failures in the
connection parts of the aluminum blades. By the late 90’s, their energy production decreased
down to 10,000 megawatt/hour and whole fleet was removed in early 2000’s. Nonetheless, a
thorough investigation of the blade dynamics was not conducted.
Next is a comparison of HAWTs and VAWTs in terms of efficiency, blade dynamics and
installation features.

3

1.1.1

Efficiency

Betz Law states that a turbine can capture up to 59.2% of the kinetic energy of the flow
[7]. When considered stand-alone, HAWTs are more efficient [8, 5, 9, 10] as can be seen
from Figure (1.3). However, when located in arrays, as in wind farms, HAWTs heavily affect
each other resulting a decrease in the total energy captured. On the other hand, several
findings [9, 11, 12] show us that VAWTs may be more effective in capturing wind energy
when positioned in wind farms. Moreover, since HAWTs have a fixed axis of rotation, their
performance is heavily dependent on wind direction. Yet, VAWT performance is not affected
with changing wind direction.

Figure 1.3 Comparison of aerodynamic efficiencies of different wind turbine types [10]

4

1.1.2

Blade Dynamics

HAWTs encounter cyclic gravitational loads on their blades. In the horizontal blade positions, gravity acts against either the leading or trailing edge. When a blade is in the bottom
horizontal position of its rotation, the gravitational load acts in tension, whereas when a
blade is in the top horizontal position of its rotation, gravitational load acts in compression. This cyclic change in the direction of load results periodic softening and stiffening
of the blade structure, as well as periodic direct forcing [13, 14]. As stated earlier, cyclic
gravitational loads may be a limiting factor for HAWTs as we keep increasing the turbine
scales. Ongoing research by Acar and Feeny [15] focuses on modeling bend-bend-twist type
deformations and resulting vibrations on HAWT blades. Also one can find more information
about equations of motion for a rotor blade in the works Kalsøe[16] and nonlinear vibration
analysis of the wind turbine blades in the works of Inoue and Ishida [14].
VAWTs don’t get affected by cyclic gravitational forces, but they also confront durability
problems due the cyclic lift and drag forces and moments. As VAWTs rotate, their blades
face the incoming flow with changing angles of attack and varying magnitudes of relative
flow velocity. Periodic changes of angle of attack create cyclic changes of lift and drag forces
and moments.
HAWT blades also experience cyclic aerodynamic loads due to rotation through a heightdependent wind profile, and due to a pulse effect during tower passes.

1.1.3

Installation

Due their geometry, the center of gravity of a HAWT is higher than the center of gravity
of a VAWT of the same size. This requires larger platforms for HAWTs, increasing the

5

material weight and cost. As stated in the Sandia Workshop of 2012, VAWTs might be more
advantageous in off-shore installations for this and a variety of other reasons [5].

1.2

Motivation

The following facts motivate the need of a comprehensive research in VAWTs [17]:
1. Not being the primary contributor of its market, our knowledge and experience with
VAWT dynamics is limited.
2. Despite all of the efforts and research, large scaled HAWTs suffer failures due cyclic
gravitational and aerodynamic effects. Therefore, being heavily dependent on HAWTs might
be a bottleneck for the future of wind power. VAWTs might reach larger scales and therefore
produce higher stand-alone energy rates in the future. Furthermore, VAWTs may perform
better in arrays[9].
3. Known durability problems for VAWT blades are outdated and might be overcome
with modern material and production technologies if investigated further.
4. A better understanding of VAWT dynamics might help us to prove and improve better
off-shore wind power generation capacities.
5. Efficiency of stand-alone VAWTs might be improved with different blade configurations
combining lift-momentum and drag cup.
However, as stated in a recent Sandia Report [5], VAWT technology is still open to
development and research in the following areas can be conducted to broaden our knowledge
and widen its usage.
1. Aerodynamics: Generally typical NACA blades are used in VAWTs. Improved aerodynamic analysis of turbines will help us find optimized blade solutions which can significantly

6

improve turbine efficiency. Similarly struts and joints can be optimized to improve performance.
2. Materials, drive train: Failures of VAWT blades in the 80’s in California are thought
to be related to the aluminum blade and joints. Further analysis of loads and vibrations
in VAWTs combined with the improvements in the material science may open the path for
durable VAWTs that can be scaled larger than HAWTs. Furthermore, once the loads and
vibrations are modeled, optimal drivetrain architecture can be developed.
3. Configurations: Different blade configurations (number of blades, locations, mixed
drag-cup, lift momentum usage etc.) can be investigated to improve the efficiency of VAWTs
[18]. See Figure 1.4.

Figure 1.4 Basic VAWT Configurations: (a) Full Darrieus, (b)“H” , (c) “V”, (d) “∆”, (e)
“Diamond” and (f) “Giromill”[6]
.
4. Off-shore designs: The geometric advantage of VAWT usage offshore is mainly due
to the lower center of gravity and zero cyclic gravitational forces. These advantages can be
investigated further to contribute widespread usage of VAWTs.

7

1.3

Objective

This research aims to derive a vibration model that accommodates the affects of varying
relative flow angle within a simple aerodynamic model. Throughout the research, a straight
vertical airfoil, such as the ones in H-Rotor and Giromill configured VAWTs, will be the
target of our interest. It will be supported by two struts located at the ends of the blade;
still, it can be modified easily for other locations on the blade.

1.4

Contribution

This research aims to contribute
1. A simple model of VAWT blade vibration, which currently does not seem to exist.
This two-mode model will consist of a pair of 2nd -order ODEs that are cyclically forced and
parametrically excited.
2. The model will be formulated to include up to cubic nonlinearity.
3. The ODEs will have many terms whose coefficients will be defined in a table of formulas
that will be provided.
4. This work serves to set up the follow-up research which will involve analysis and
numerical simulations.

1.5

Thesis Outline

In the upcoming chapters:
1. Energies for deriving nonlinear partial differential equations of the vibration of a single
blade will be obtained for one-dimensional transverse bending and twist.

8

2. A simplified aeroelastic model will be introduced based on semi-empirical quasi-static
airfoil theory.
3. Single assumed modes for transverse bending and twist will be used to reduce the order
of equations in the model. Two ordinary differential equations representing bend and twist
will be derived using Lagrange’s equations. The equations will be linearized relative to the
undeformed state.
4. A numerical evaluation of equation coefficients will be made for a specifically defined
turbine blade, resulting in equations that are dependent on few parameters, such as rotation frequency and blade length. Since data for an existing H-rotor/Giromill turbine is
not available, we will use data from a Sandia 17m Darrieus turbine as a guide for our Hrotor/Giromill model. The numerical evaluation will include an estimate of the first-mode
natural frequency, as a function of blade length, for the non rotating blade. The aeroelastic
terms will be based on a NACA0012 symmetric airfoil.

9

Chapter 2
Formulation and Modeling
We derive a vibration model for an H-rotor/Giromill blade. The blade is treated as a
uniform straight beam under bending and twisting deformations.
This chapter can be presented in several parts. First comes the derivation of the energy
equations for the bending and twisting blade. It will be followed by a simplified aerodynamic
model to obtain the lift and drag force and moment equations as functions of time and spatial
location on the blade.
Before going forward with modeling several assumptions were made:
1. Steady wind speed
2. Neglect the disturbance of wind due to upstream blade or post
3. Constant (steady-state) rotation rate
4. Inextensible beam
5. Shear center is at the geometric center of the blade

2.1

Elastic Blade Model for Bending

Following the ideas of Feeny [17] we will model our blades as an elastic Euler-Bernoulli beam.
A blade that we model could face bend-bend-twist deflections, yet only one-dimensional
transverse bending and twist will be taken into consideration for the rest of the formulation
and modeling. Here transverse bending is the deflection of the blade to the radial direction,

10

and twist occurs around the axis that passes through the centroid of the blade cross sections.
First, kinetic and potential energy equations for the particles on the beam will be expressed such that the extended Hamilton’s principle can be used to obtain the partial differential equations of motion. Alternatively, a Lagrange formulation based on assumed modes
can be applied to the energy formulations.
In order to do so, sections along the undeformed beam axis will be located from an origin
at the middle center of the rotor using the following position vector:

r(x, t) = (x − s(x, t))i + (R + y(x, t))ˆ
er

(2.1)

See Figure 2.1 and Figure 2.2.
In this equation, R is the distance between undeformed beam and the axis of rotation,
x is the vertical distance of the blade section to the origin, y(x, t) is the radially transverse
displacement of the selected section on beam at time t, s(x, t) is the cumulative foreshortening
affecting point x at time t, due deflection. The foreshortening term, s, is formulated to be
positive for x > 0. Also i is a unit vector in the direction in the undeformed x axis and ˆ
er
is a unit vector in the direction in the y axis, radial to the center of rotation.
We can also formulate the foreshortening for an Euler-Bernoulli beam due to one dimensional transverse bending in the radial direction y to retain up to 4th degree terms, as
[13]
s(x, t) ∼
=

y2 y4
+
)dz,
2
8
−a
+a

(

(2.2)

which will contribute linear and cubic terms to the equations of motion where a is half
the length of the blade. Assuming constant rotational speed and neglecting the effects of
transverse deflection in the circumferential direction on the energy, the kinetic energy for
11

Figure 2.1 A sketch of blade and rotor. Undeformed on left and deformed on right.
the particles on the beam can be formulated with
+a

+a 1
1
1
2
T =
m(x)v(x, t) dx +
Jxx (x)(ψ˙ + Ω)2 + I0 θ˙2 ,
2
−a 2
−a 2

(2.3)

where I0 is the mass moment of inertia of the rotor, m(x) is the mass per unit length of the
beam, and ψ is the angle of twist deformation. Furthermore, v is its velocity which can be
found by taking time derivative of the position function r(x, t) as

v = r˙ = −sˆ
˙ eθ + yˆ
˙ er + Ω(y(x, t) + R)k,

(2.4)

where Ω is the constant angular speed of the rotating parts.
Similarly, potential energy can also be formulated in the form of integrations. We will
express the total potential energy as V which is the summation of the gravitational potential

12

energy, Vmg , strain energy due bending and twist, Vs , and the elastic potential energy of the
boundaries of the beam, Vb .
As such, neglecting net vertical deflections of the blade, we have

+a

Vmg =

m(x)gs(x, t)dx,

(2.5)

−a

+a

Vs =

+a 1
1
EIz (x)[y 2 (1 − 3y 2 )]dx +
GIx (x)ψ 2 dx,
−a 2
−a 2

(2.6)

1
1
1
1
Vb = k1 (ug − s(−a, t))2 + kT1 (y (−a, t))2 + k2 (ug − s(a, t))2 + kT2 (y (a, t))2 , (2.7)
2
2
2
2
where E is the Young’s modulus, and Iz (x) is the area moment of inertia of the cross
section at x about the neutral axis. Vs includes effects of nonlinear bending curvature [19]
approximated to 4th degree. Also one should note that elastic modeling has two parts for
each strut. Struts are modeled as superposed linear and torsional springs where linear spring
constants are k1 and k2 , torsional spring constants are kT1 and kT2 for struts 1 and 2. In
order to simplify the model, k1 is assumed to be equal to k2 and kT1 is assumed to be equal
to kT2 .
The extended Hamilton’s principle can be used to obtain partial differential equations of
motion. Hamilton’s principle states that
t2

(δT − δV + δWnc )dt = 0

(2.8)

t1

or
t2
t1

L
0

ˆ nc )dx + δTb − δVb + δWnc
(δ Tˆ − δ Vˆ + δ W
b

dt = 0

(2.9)

where δ Tˆ, δ Vˆ are energy densities, defined as the integrands in eqns (2.3), (2.5) and (2.6),

13

and δTb , and δVb are boundary terms defined as the lumped terms in equations (2.3) and
(2.7).
Here Wnc is the nonconservative virtual work acting on the airfoil. Since the forces are
distributed over the length of the blade,
a

δWnc =

−a

ˆ nc dx
δW

(2.10)

Alternatively, we can apply Lagrange’s equations on assumed modal coordinates, as in
the next section.

2.2

Reduced Order Modeling

In the previous section, we have introduced our model. As mentioned, a transverse bending
and twist will be our sources of deflections on our rotating blade, which is modeled as an
Eulerian beam. We will go further to simplify our model by using reduced order modeling.
Deflections on the beam will be projected on assumed modes, based on hinged-hinged beam
modes to obtain second-order nonlinear ordinary differential equations. Using an assumed
mode formulation [20],

N

τk (x, t) ≈

ζi (t)ηi (x)

(2.11)

i=1

where τ is the deflection, N is the number of assumed modes, ηi (x) are the assumed modal
functions, and ζi (t) are the assumed modal coordinates. We can write

y(x, t) ∼
=

Nb

qi (t)ξi (x)
i=1

14

(2.12)

Figure 2.2 Top view representation of transverse deflection on top and twist on bottom.

ψ(x, t) ∼
=

Nt

bi (t)ρi (x)

(2.13)

i=1

where Nb and Nt are the numbers of bend and twist assumed modes. We will use one
assumed mode for both, so will have

y(x, t) ≈ q(t)ξ(x)

(2.14)

ψ(x, t) ≈ b(t)ρ(x)

(2.15)

˙ t) ≈ b(t)ρ(x),
˙
As such, y(x,
˙ t) ≈ q(t)ξ(x),
˙
ψ(x,
y (x, t) ≈ q(t)ξ (x), ψ (x, t) ≈ b(t)ρ (x).
Next we can express our energy functions

˙ t) =
T (q, q,
˙ b, b,

+a

1
2 + (Ω(y(x, t) + R))2 dx
m(x) (−s)
˙ 2 + q(t)ξ(x)
˙
2
−a
+a 1
2
1
˙
+
Jxx (x) b(t)ρ(x)
+ Ω + I0 θ˙2 ,
2
−a 2
15

(2.16)

+a

Vmg (q, q,
˙ t) =
˙ t) =
Vs (q, q,
˙ b, b,

m(x)gs(q, q,
˙ t)dx,

+a

1
EI(x)
−a 2

(2.17)

−a

q(t)ξ (x)

2

1 − 3 q(t)ξ (x)

2

dx

+a

1
2
+
GJxx (x) b (t)ρ(x) dx,
−a 2
˙ t) = 1 k1 ug − s(−a, t) 2 + 1 kT q (t)ξ(−a) 2
Vb (q, q,
˙ b, b,
2
2 1
1
1
2
2
+ k2 ug − s(a, t) + kT2 q (t)ξ(a)
2
2

(2.18)

(2.19)

Now, the energy equations on assumed modes can be plugged into Lagrange’s equation:

∂ ∂T
∂T
∂V
(
)−
+
= Qk
∂t ∂ q˙k
∂qk ∂qk

(2.20)

to obtain equations of motion. Here T stands for the kinetic energy, V stands for the potential
energy, and qk is the generalized coordinate (dependent variable) from k = 1, 2..., n. Qk is
the generalized force term, and

n

δWnc =

Qk δqk

(2.21)

k=1

Applying eqns. (2.14) and (2.15) for the case of a two assumed mode model

δWnc = Qq δq + Qb δb

16

(2.22)

In our case of bend and twist deflections, application of Lagrange’s equations leads to
a

x

2q q˙

m(x)
−a

−q q˙2

a

(ξ )2 dx dx + q 2 q˙

x

m(x)
−a

0
a

x

m(x)

(ξ

−a

)2 dx

a

dx − q

+q

m(x)Ω2 ξ 2 dx −

−a
a

gm(x)
−a

(ξ )2 dx dx + q 3

0
a

q
k + k2
+q 3 1
2

¨b

a
−a

m(x)
2
−a

0
a
−a

a

m(x)Ω2 Rξdx

−a
x

g

0
a

(ξ )4 dx dx+

(2.23)

6EI(x)ξ 2 ξ 2

−a
2

(ξ )2

Ixx ρ2 dx + b

a

EI(x)ξ 2 dx − q 3

−a
±a

m(x)ξ 2 dx

−a

0

0
a

a

(ξ )2 dx dx + q¨

+ q(kT1 + kT2 )(ξ |±a )2 = Qq

Jxx Gρ 2 dx = Qb

(2.24)

where Qq and Qb are non-conservative aeroelastic forces to be determined in the next section.

2.3

Aeroelastic Modeling

One of the main problems with VAWTs that requires further investigation is the cyclic forces
acting on the blades. Under a constant rotation rate of the rotor and blades, the relative
velocity and angle of the wind particles hitting the blade is changing periodically, which
creates cyclic lift and drag forces and moments on the blade. These forces are associated
with the nonconservative work, Wnc , in the Hamilton’s principle (2.8).
In order to simplify our aeroelastic model, our assumption is a flow field with constant
direction and velocity. This assumption also neglects the effects of the blades on the fluid
particles. Yet one should note that blades also have their own velocity due to rotation
and deflection. Therefore the relative wind speed has a cyclic and displacement dependent
behavior.
17

Figure 2.3 shows a sketch of the turbine rotating under the effect of constant wind speed
u, such that the wind velocity is u = −uj, where j is a unit vector in the y direction in Figure
2.3. Note that in Figure 2.3, the y axis is fixed, and is not to be confused with transverse
deflection, y(x, t), which is in the er direction that rotates with the turbine. er and eθ are
radial and tangential unit vectors. θ is the coordinate that represents the orientation of the
unit vector eθ with respect to the y-axis. Therefore j = er sin θ + eθ cos θ. Furthermore we
˙ θ and e˙ θ = −θe
˙ r.
let θ, θ˙ = Ω. Therefore, e˙ r = θe

Figure 2.3 Cross section of VAWT showing the coordinate axes

Considering the transverse deflection, the distance r between the axis of rotation and
blade’s center of gravity can be notated in radial and tangential coordinates as r = (R +
y(x, t))er and by differentiating r we will obtain the velocity of the blade as vb = r˙ =
˙ θ.
ye
˙ r + (R + y)θe
The flow velocity relative to the blade will be obtained as v = u − vb = −uj − vb , and
hence
v = −(y˙ + u sin θ)er − [(R + y)θ˙ + u cos θ]eθ

(2.25)

The angle between the relative flow velocity and tangential axis eθ is going to be denoted
18

Figure 2.4 Twist angle ψ and relative wind angle φ for the case when the pretwist β = 0
with φ, where
tan φ = vr /vθ =

y˙ + u sin θ
(R + y)θ˙ + u cos θ

(2.26)

One should note that angle of attack α is a summation of relative wind velocity angle, φ,
pre-twist angle, β, and twist deformation angle ψ, i.e.

α=ψ+β+φ

(2.27)

where Figure (2.4) shows a sketch for the case β = 0.
Using the Taylor series expansion for φ for small values of vr /vθ

φ = arctan

y˙ + u sin θ
y˙ + u sin θ
1
≈
−
(R + y)θ˙ + u cos θ
(R + y)θ˙ + u cos θ 3

3
y˙ + u sin θ
(R + y)θ˙ + u cos θ

(2.28)

Now, we will expand tan φ as a mutivariable power series for small variations of µk = y, y,
˙ z, z˙
and let θ = Ωt:
3

tan φ =

3

3

y˙ + u sin θ
u sin Ωt
=
+
µi
µj
µk dijk
RΩ + u cos Ωt
(R + y)θ˙ + u cos θ
i=0
j=0
k=0

19

(2.29)

where µ0 = 1 and d000 = 0. The lift and drag forces per unit length acting on the blade are

Figure 2.5 Lift and drag force per unit length directions from top view

1
L = CL (α)cρν 2 (− cos φer + sin φeθ ) = Lr er + Lθ eθ
2

(2.30)

1
D = CD (α)cρν 2 (− sin φer − cos φeθ ) = Dr er + Dθ eθ
2

(2.31)

in the directions that can be seen on Figure (2.5), where CL and CD are coefficients of lift and
drag respectively, ρ is the density of air and c is the chord length. CL and CD are empirically
determined functions of α, and for a symmetric airfoil we can assume CL (α) ≈ c1 α − c2 α3
and CD (α) ≈ d1 + d2 α2 , where c1 , c2 , d1 and d2 are positive constants to be empirically
determined case by case. For example, CL and CD for the case of NACA0012 are shown in
Figure 2.6 from reference [21].
Since α is dependent on displacement and velocity variables through equations (2.25)
and (2.26), CL and CD can be considered as functions of y, y,
˙ and ψ, as CL (α(y, y,
˙ ψ))
and CD (α(y, y,
˙ ψ)). In addition, the relative flow speed is a function of variables, such that
v(y, y)).
˙
Thus,
20

Figure 2.6 NACA0012 CL and CD values with respect to angle of attack [21]

CL ≈ c1 α(y, y,
˙ ψ)3 + c2 α(y, y,
˙ ψ)

(2.32)

CD ≈ c3 α(y, y,
˙ ψ)2 + c4

(2.33)

v 2 = (y˙ + u sin θ)2 + [(R + y)θ˙ + u cos θ]2

(2.34)

and

As such, the aerodynamic loads can be expressed in terms of y,y,
˙ and ψ,ψ˙ by inserting
equations (2.27) and (2.29) into equations (2.30) and (2.31), and also using equations (2.32),
(2.33), and (2.34). This results in
21

L ≈ Lr (y, y,
˙ ψ)er + Lθ (y, y,
˙ ψ)eθ

(2.35)

D ≈ Dr (y, y,
˙ ψ)er + Dθ (y, y,
˙ ψ)eθ

(2.36)

In our case of transverse deflection and twist, the lift and drag forces and moments
ˆ nc term in equation (2.10) such that
contribute to the δ W

˙ t)δy + M
˙ t)δψ
ˆ nc = fˆ(y, y,
ˆ (y, y,
δW
˙ ψ, ψ,
˙ ψ, ψ,

(2.37)

where radial forces will define function fˆ such as

fˆ(y, y,
˙ ψ, t) = Lr (y, y,
˙ ψ, t) + Dr (y, y,
˙ ψ, t)

(2.38)

ˆ as
and a lift moment function M

ˆ (y, y,
M
˙ ψ, t) = (L cos α + D sin α) lc

(2.39)

where lc is the distance between the aerodynamic center of the airfoil and its center of gravity.
One will find that lift and drag are cyclic functions, and are time periodic in the case
where Ω is constant. We will use the multi-variable Taylor series expansion (Appendix) up
ˆ (y, y,
to cubic terms on our functions fˆ(y, y,
˙ ψ, t) and M
˙ ψ, t) around (y,y,ψ)=(0,0,0).
˙
As

22

such,

fˆ(y, y,
˙ ψ, t) = fˆ000 + fˆ100 y + fˆ010 y˙ + fˆ001 ψ + fˆ200 y 2 + fˆ020 y˙ 2 + fˆ002 ψ 2
+fˆ110 y y˙ + fˆ101 yψ + fˆ011 yψ
˙ + fˆ300 y 3 + fˆ030 y˙ 3 + fˆ003 ψ 3
(2.40)
+fˆ210 y 2 y˙ + fˆ201 y 2 ψ + fˆ021 y˙ 2 ψ + fˆ120 y y˙ 2
+fˆ102 yψ 2 + fˆ012 yψ
˙ 2 + f111 y yψ
˙

ˆ (y, y,
ˆ 000 + M
ˆ 100 y + M
ˆ 010 y˙ + M
ˆ 001 ψ + M
ˆ 200 y 2 + M
ˆ 020 y˙ 2 + M
ˆ 002 ψ 2
M
˙ ψ, t) = M
ˆ 110 y y˙ + M
ˆ 101 yψ + M
ˆ 011 yψ
ˆ 300 y 3 + M
ˆ 030 y˙ 3 + M
ˆ 003 ψ 3
+M
˙ +M
(2.41)
ˆ 210
+M

y 2 y˙

ˆ 201
+M

y2ψ

ˆ 021
+M

y˙ 2 ψ

ˆ 120
+M

y y˙ 2

ˆ 102 yψ 2 + M
ˆ 012 yψ
ˆ 111 y yψ
+M
˙ 2+M
˙
ˆ ijk are time periodic.
where the coefficients fˆijk and M

Then we will substitute y =

q(t)ξ(x) and ψ = b(t)ρ(x) to write fˆ(y(q, ξ(x)), y(
˙ q,
˙ ξ(x)), ψ(b, ρ(x)), t) = f (q, q,
˙ b, t; x) and
ˆ (y(q, ξ(x)), y(
M
˙ q,
˙ ξ(x)), ψ(b, ρ(x)), t) = M (q, q,
˙ b, t; x), and integrate according to equation
(2.10) to write
a
−a

f (q, q,
˙ b, t; x)ξ(x)dx = f000 + f100 q + f010 q˙ + f001 b + f200 q 2 + f020 q˙2 + f002 b2
+f110 q q˙ + f101 qb + f011 qb
˙ + f300 q 3 + f030 q˙3 + f003 b3
+f210 q 2 q˙ + f201 q 2 b + f021 q˙2 b + f120 q q˙2
+f102 qb2 + f012 qb
˙ 2 + f111 q qb
˙

23

(2.42)

a
−a

M (q, q,
˙ b, t; x)ρ(x)dx = M000 + M100 q + M010 q˙ + M001 b + M200 q 2 + M020 q˙2 + M002 b2
+M110 q q˙ + M101 qb + M011 qb
˙ + M300 q 3 + M030 q˙3 + M003 b3
+M210 q 2 q˙ + M201 q 2 b + M021 q˙2 b + M120 q q˙2
+M102 qb2 + M012 qb
˙ 2 + M111 q qb
˙
(2.43)

where the coefficients fijk and Mijk are cyclic in time, and can be found in the Appendix.
Using equations (2.14) and (2.15), δy = ξ(x)δq and δψ = ρ(x)δb. Then the nonconservative work of equation (2.37) can be expressed in terms fˆ(y(q, ξ(x)), y(
˙ q,
˙ ξ(x)), ψ(b, ρ(x)), t)
ˆ (y(q, ξ(x)), y(
=f (q, q,
˙ b, t; x) and M
˙ q,
˙ ξ(x)), ψ(b, ρ(x)), t) = M (q, q,
˙ b, t; x) as

˙ t; x)ξ(x)δq + M (y, y,
˙ t; x)ρ(x)δb
ˆ nc = f (y, y,
δW
˙ ψ, ψ,
˙ ψ, ψ,

(2.44)

As such, considering equations (2.8) and (2.10) the generalized forces in equations (2.23)
and (2.24) have the form

a

Qq =

f (q, q,
˙ b, t; x)ξ(x)dx

(2.45)

M (q, q,
˙ b, t; x)ρ(x)dx

(2.46)

−a

a

Qb =

−a

which are expressed in terms of the expansions in equations (2.42) and (2.43).
Next we will linearize the equations that we obtained.

24

2.4

Linearization

In order to conduct a simple initial numerical analysis and simulate our model we will
decrease the complexity of our model. In order to do so, several linearization methods were
discussed.
Since the system has periodic excitation, equilibria do not exist. If we consider finding “unforced” equilibria, we can drop the direct excitation. The resulting equations are
nonlinear, with many terms, and with parametric excitation; so equilibria may not exist.
If we consider the equations when all cyclic time varying terms are omitted, the resulting
equations are still nonlinear with many terms, and so the equilibrium will be difficult to
express, and to use as a reference point. As such, we linearize about the equilibrium of the
non rotating system, which is zero. Thus, the “linearized” model we consider is obtained by
assuming y and ψ are small, hence quadratic and cubic terms in q and b and their derivatives
are dropped.
The resulting equations for small deflections of the rotating system within direct and
parametric excitations are

q¨ + ωq2 q + qa1 (t) + qa
˙ 2 (t) + ba3 (t) = a0 + f0 (t)

(2.47)

¨b + ω 2 b + qe1 (t) + qe
˙ 2 (t) + be3 (t) = e0 + g0 (t)
b

(2.48)

where
a mΩ2 ξRdx
a0 = −aa
2
−a mξ dx

25

a f
2
− −a
100 ξ (x)dx
a
2
−a mξ dx
a f
2
− −a
010 ξ (x)dx
a2 (t) =
a mξ 2 dx
−a
a f
− −a
001 ξ(x)ρ(x)
a3 (t) =
a mξ 2 dx
−a
a f
000 ξ(x)dx
−a
f0 (t) =
a mξ 2 dx
−a
a mΩ2 ξ 2 dx+ a gm x ξ 2 dx dx+ a EIξ 2 dx+k (ξ | )2 +(ξ | )2
− −a
a
−a
T
−a
−a
0
a mξ 2 dx
−a

a1 (t) =

ωq2 =
and,

e0 = 0
a M
2
− −a
100 ξ (x)dx
a
2
−a Jxx ρ dx
a M
2
− −a
010 ξ (x)dx
e2 (t) =
a
2
−a Jxx ρ dx
a
− −a M001 ξ(x)ρ(x)
e3 (t) =
a
−a Jxx
a M
ξ(x)dx
g0 (t) = −aa 000 2
−a Jxx ρ dx
a I G(ρ )2 dx
xx
ωb2 = −aa
2
−a Jxx ρ dx

e1 (t) =

26

Chapter 3
Numerical Analysis and Simulation

3.1

Definition of the Giromill/H-Rotor Blade

After formulating the energy equations and modeling the nonconservative forces to obtain
the equations of motion, now we want to make a vibration simulation and numerical analysis
for our blade. Instead of focusing on the effects of different variables on VAWT blade
vibration, we will focus simulating the behavior of a selected type of blade and accommodate
reliability of our models.
As already stated, our knowledge and experience with VAWT is limited. In order to
validate our model we need experimental data for reference. For that reason we will focus
on the blades of 17 meter diameter Darrieus wind turbine of Sandia (Sandia 17m) [22] which
were used by FloWind in California Altamont and Tehachipi wind farms.
Experimental data obtained from Sandia 17m shows us that maximum turbine performance is obtained around a tip ratio, RΩ/u, of 5 [22, 18]. Moreover, Sandia’s performance
tests are handled for rotor speeds varying from 29.6 rpm to 52.5 rpm, where one can see
further power output analysis on rotor speeds of 37 rpm and 48 rpm [22]. Moreover, their
addendum for 3 bladed turbine also focused on a rotor speed of 37 rpm [18]. Therefore, for
the parts that require a constant rate of rotation, we can assume a constant rate of rotation
of 37 revolution per minute or similarly Ω = 3.875 rad/s as a representative example for the
numerical study.
27

The Sandia Darrieus wind turbine uses NACA0012 symmetric airfoil blades of 21 inch,
or 0.53 meter, chord length that has a surface equation of [23]

yˆ(ˆ
x) = ±

0.12
.53 0.2969
0.2

xˆ
− 0.1260
.53

xˆ
.53

− .3516

xˆ 2
+ 0.2843
.53
−0.1015

xˆ 3
.53
xˆ 3
.53

(3.1)

where xˆ is the horizontal location along the chord line of the airfoil and yˆ defines the surface
of the airfoil as a function of xˆ. This equation will be used with a uniform thickness to
estimate the material cross sectional mass center, area moment of inertia and mass moment
of inertia per unit length.
Our energy and aeroelastic models are based on Giromill/H-rotor style VAWTs which
have straight blades that are parallel to axis of rotation. However, Darrieus VAWT blades
have a troposkein curved shape referred as an ”egg beater”[6]. Each Sandia 17m blade is
constructed from two straight parts and one circular part attached together. See Figure 3.1.

Figure 3.1 Sandia 17m front view
28

The cross section of a Sandia 17m blade has a joint construction of extruded aluminium,
Nomex core, and fiberglass skin. See Figure 3.2.

Figure 3.2 Sandia 17m blade cross section showing aluminum spar and spline, Nomex core
and fiberglass skin. [22]

We will focus on the straight sections of the Sandia 17m and approximate the blades to
be constructed only of extruded aluminum. Straight sections of Sandia 17m weigh 10.137
lb/ft and thus have a mass per unit length of m = 15.09 kg/m (Appendix) [22]. In order
to have same mass/length ratio, based on equation (3.1) the wall thickness of our extruded
aluminum airfoil with density of ρaluminium = 2700 kg/m3 should be 0.2 inches or rounded
to 0.005 meters. See Figure 3.3.

Figure 3.3 Blade dimensions for calculations

Using the NACA0012 with extruded aluminum of thickness 0.005 m and chord length of
0.53 m, see Figure 3.3, we can calculate the centroid of the blade cross section to be located
0.263 m from the front tip and the distance between centroid, and aerodynamic center, 1/4
of chord length for two dimensional incompressible flow [24], of the blade, lc in Equation
29

(2.39), is going to be 0.263 − 0.53/4 = 0.13 meters. Furthermore, we can calculate the mass
moment of inertia per unit length around centroid of blade cross section as Jzˆ = 0.37 kg.m,
second moment of area around the centroid of the blade as Izˆ = 0.000137 m4 and second
moment of area around chord line as Ixˆ = 3.1034 ∗ 10−6 m4 . Also we consider an aluminum
alloy with a density of 2700 kg/m3 , Young’s modulus of E = 70 GPa and shear modulus of
G = 26 GPa [25].
Table 3.1 Selected aluminium alloy’s material properties
Property
Value
Density
2700 kg/m3
Young’s Modulus
70 GPa
Shear Modulus
26 Gpa

In order to numerically analyze the equations of motion of an undeformed blade, the
radius R should be fixed. We don’t have configuration data for a mid-sized or large Hrotor/Giromill turbine, so we will use Sandia’s data on a 17m Darrieus blade as a reference
and guide for establishing a hypothetical H-rotor/Giromill of about the same size. Sandia
17m has a varying radius from 0 to 27.44 feet at maximum [22]. Yet Giromill/H-rotor type
VAWTs have constant undeformed blade radii. Approximating with a half circle shape for
Sandia 17m blade, one can calculate an average radius either for the lift force, L, or effective
torque by lift on rotor, T . Our calculations gave us an average radius of 19.40 feet in order
to obtain the same lift and an average radius of 20.62 feet in order to obtain the same torque
on rotor. Referring to these values obtained we will assume a constant radius R of 20 feet,
which is 6 meters, for our Giromill/H-rotor turbines. A summary of blade parameters is
given in Table (3.2).

30

Table 3.2 Specifications of the model created
Blade

NACA0012

Material and Construction

Extruded Aluminum Shell

Radius

6m

Cross sectional weight

15.09 kg/m

Mass moment of inertia around centroid per length

0.37 kg.m

Second moment of area around chord line

3.1034x10−6 m4

Second moment of area around centroid

0.000137 m4

Tip Speed Ratio

5

3.2

Modal Frequencies

Now we would like to see the effects of varying blade length and rotation rate on the natural
frequencies of the blade. For this purpose we assume that the blade is rotating under no
aeroelastic forces, such that Qq and Qb are equal to zero. Assuming zero torsional stiffness
on struts, i.e. the strut connections are ideal pin connections, Equations (2.47) and (2.48)
will be
a

q¨
−a

mξ 2 dx − q

a

mΩ2 ξ 2 dx −

−a

a

mΩ2 ξRdx + q

−a

a

x

gm
−a

ξ 2 dx dx+

0
a

q

(3.2)
EIξ

2 dx

=0

−a

¨b

a
−a

Jxx ρ2 dx − b

a
−a

31

Ixx G(ρ )2 dx = 0

(3.3)

where ξ and ρ are approximated by cos πx/2a for the first mode shapes. Note that the linear
conservative model is decoupled, since the shear center was assumed to be at the geometric
center of the cross section of the blade.
Plots of bending natural frequency for varying values of Ω and a are given in Figure 3.4
for varying half length a and rotor speed Ω.

Figure 3.4 Variation of wnq with Ω on left and blade length a on right

The approximate model frequencies for bend and twist for a stationary blade (rotor
speed Ω = 0)) are shown in Figure 3.5 as functions of half length a. Note that the torsional
frequency is much higher than the transverse bending frequency. This is consistent with
individual models of bending torsion [20].

3.3

Blade with Quasistatic Aeroelastic Airfoil Model

It is mentioned that lift and drag coefficients with varying angle of attack differ with
the type of blade used. Furthermore, we already stated that CL (α) ≈ c2 α + c1 α3 and
CD (α) ≈ c4 + c3 α2 . Using Figure 2.6 and curve fitting the data points we will obtain

32

Figure 3.5 Variation of wnq and wnb with a when Ω = 0
following equations for lift and drag coefficients:

CL (α) ≈ 4.4287α − 2.9916α3

(3.4)

CD (α) ≈ 0.0094 + 1.185α2

(3.5)

These values contribute to load coefficients fijk in equations (2.42) and (2.43) through equations (2.25), (2.27), and (2.28)-(2.34), and the formulas listed in the appendix. The offset of
the lift and drag forces is used to express the lift moment, based on the approximation that
the offset is constant.
As such, the coefficients shown in equations (2.47) and (2.48) are listed in the appendix
equations for specified values of parameters, but for arbitrary Ω, u, and a.
The key point to note is that the coefficients a1 (t) through a3 (t) and e1 (t) through
e3 (t) are complicated functions of time, each with fundamental frequency Ω and harmonics,
representing parametric excitation terms. The terms a0 and f0 (t) and e0 and g0 (t) represent
constant plus periodic direct excitations.

33

Figure 3.6 Rotor dimensions used in calculations
Thus, for the bend-twist VAWT blade vibration model, we derived an approximate lowerorder set of linearized equations that can now be verified, analyzed and simulated for resonances and instabilities. This sets up the model for future analysis.

34

Chapter 4
Conclusion and Comments

4.1

Summary of Results

As a result of the research, a blade vibration model for an H-rotor/Girmoill type VAWT
for bending and twisting deflections was formulated. In addition, energy equations for an
Eulerian beam under transverse bend and twist deflections was obtained. The system was
discretized using reduced order modeling on first assumed modes on each independent variable in energy equations. Lagrange’s method was applied to resulting equations to obtain
two equations of motion where generalized force terms were due aeroelastic forces on blades.
Next, an aeroelastic model was defined based on quasistatic airfoil theory. Lift and drag
forces and moments were formulated for an airfoil for changing angle of attack, where stall
effects were neglected. Resulting complicated formulas were simplified to cubic order using
Taylor series expansion. One should note that, resulting system had parametric and direct
excitation due to varying flow magnitude and direction relative to blade. Very complicated
terms were tabulated based on how to generate them with computer algebra.
In order to conduct a simple numerical analysis the system was linearized assuming small
deflections for bend and twist. Linearized equations of motion were derived for a specific
blade. Referring to Sandia 17m VAWT, hypothetical H-rotor/Giromill blade was defined
for numerical analysis where natural frequencies of the blade for a non rotating system were
found.
35

Furthermore, initial simulations showed complicated dynamics and propensity for instability.

4.2

Contribution

In a large absence of VAWT blade vibration models, this work has provided the initial steps
at modeling VAWT blade vibrations to the research literature.
The model presented here can be used as a foundation for subsequent vibration modeling
research on VAWTs. Research to be done while building this model is described in the next
section on Future Work.

4.3

Future Work

This work sets up follow-up work on a systematic study of vibrations and how responses
depend on parameters. To this end, future work will include
• verification of the terms which are very complicated
• perturbation analysis and numerical simulations to determine the potential resonances
and instabilities, and the effects of parameters on these phenomena
• higher degree of freedom reduced model, where more than one assumed mode will be
used in the calculations
• improvement of the aeroelastic model, i.e. including the effects of stall and the hysteretic effects caused by the influence of the changing angle of attack on the local flow
field.

36

• improvement of the structural model, for example in separating the shear center from
the geometric center, and including bend-bend-twist deflections.
• nonlinearity back into the model and perform nonlinear analyses and simulations
• modeling of the curved blades of Darrius-type turbines

37

APPENDIX

38

The following terms are obtained from Taylor Series expansion of radial force component
f and moment M in equations (2.42) and (2.43) where c1 , c2 , c3 and c4 are constants from
lift and drag coefficient formulations and dependent on the geometry of the airfoil.

f000 = 0.5ρair c(ΩR + u cos θ) (ΩR + u cos θ)2 + u2 sin2 θ
c1 β −

u3 sin3 θ
3(ΩR+u cos θ)3

+0.5ρair cu sin θ

f100 =
c1 β −

+

u sin θ
ΩR+u cos θ

(ΩR + u cos θ)2

3

+ c2 β −

+ u2 sin2 θ

u3 sin3 θ
3(ΩR+u cos θ)3

c3 β −

u sin θ
+ ΩR+u
cos θ

u3 sin3 θ
3(ΩR+u cos θ)3

+

u sin θ
ΩR+u cos θ

0.5ρair cΩ(ΩR+u cos θ)2
(ΩR+u cos θ)2 +u2 sin2 θ
u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

3

+ c2 β −

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

+ 0.5ρair c(ΩR + u cos θ) (ΩR + u cos θ)2 + u2 sin2 θ
3c1
+c2

Ωu3 sin3 θ
(ΩR+u cos θ)4
Ωu3 sin3 θ
(ΩR+u cos θ)4

c1 β −

+

−

−

Ωu sin θ
(ΩR+u cos θ)2

Ωu sin θ
(ΩR+u cos θ)2

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

β−

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

2

+ 0.5ρair cΩ (ΩR + u cos θ)2 + u2 sin2 θ
3

+ c2 β −

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

2
u3 sin3 θ
u sin θ
0.5ρair cΩu sin θ(ΩR+u cos θ) c3 β−
+
+c4
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)2 +u2 sin2 θ

+ ρair c3 cu sin θ
β−

(ΩR + u cos θ)2 + u2 sin2 θ

u3 sin3 θ
3(ΩR+u cos θ)3

f010 =
c1 β −

Ωu3 sin3 θ
(ΩR+u cos θ)4

−

Ωu sin θ
(ΩR+u cos θ)2

u sin θ
+ ΩR+u
cos θ

0.5ρair cu sin θ(ΩR+u cos θ)
(ΩR+u cos θ)2 +u2 sin2 θ
u3 sin3 θ
3(ΩR+u cos θ)3

+

u sin θ
ΩR+u cos θ

3

+ c2 β −

u3 sin3 θ
3(ΩR+u cos θ)3

+ 0.5ρair c(ΩR + u cos θ) (ΩR + u cos θ)2 + u2 sin2 θ
39

u sin θ
+ ΩR+u
cos θ

2

+ c4

3c1

1
ΩR+u cos θ

+c2

+

1
ΩR+u cos θ

−

−

u2 sin2 θ
(ΩR+u cos θ)3

u3 sin3 θ
3(ΩR+u cos θ)3

β−

u sin θ
ΩR+u cos θ

+

2

u2 sin2 θ
(ΩR+u cos θ)3

2
u3 sin3 θ
u sin θ
0.5ρair cu2 sin2 θ c3 β−
+
+c4
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)2 +u2 sin2 θ
3

3

u sin θ
u sin θ
+ 0.5ρair c (ΩR + u cos θ)2 + u2 sin2 θ c3 β −
+ ΩR+u
cos θ
3(ΩR+u cos θ)3

+ ρair c3 cu sin θ
β−

1
ΩR+u cos θ

(ΩR + u cos θ)2 + u2 sin2 θ

u3 sin3 θ
3(ΩR+u cos θ)3

−

2

+ c4

u2 sin2 θ
(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

f001 = 0.5ρair c(ΩR + u cos θ) (ΩR + u cos θ)2 + u2 sin2 θ
3c1 β −

u3 sin3 θ
3(ΩR+u cos θ)3

+ ρair c3 cu sin θ

+

2

u sin θ
ΩR+u cos θ

+ c2

(ΩR + u cos θ)2 + u2 sin2 θ β −

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

 f200 = 0.25ρair c

Ωu sin θ
Ωu3 sin3 θ −
u3 sin3 θ
+ u sin θ +β
4c3 Ωu(ΩR+u cos θ) sin θ
−
4
2

(ΩR+u cos θ)
(ΩR+u cos θ)
3(ΩR+u cos θ)3 ΩR+u cos θ

(ΩR+u cos θ)2 +u2 sin2 θ

+c3 u sin θ
2
−

(ΩR + u cos θ)2 + u2 sin2 θ

Ωu3 sin3 θ
(ΩR+u cos θ)4
u3 sin3 θ
3(ΩR+u cos θ)3

−

2

Ωu sin θ
(ΩR+u cos θ)2

2Ω2 u sin θ
(ΩR+u cos θ)3

+2

2 3 3
− 4Ω u sin θ5
(ΩR+u cos θ)

u sin θ
+ ΩR+u
cos θ + β


+u sin θ 
c3


Ω2
(ΩR+u cos θ)2 +u2 sin2 θ

3 3
− u sin θ 3
3(ΩR+u cos θ)

+ 0.5ρair c 2

+

−

u sin θ
ΩR+u cos θ

Ω2 (ΩR+u cos θ)2
3/2
(ΩR+u cos θ)2 +u2 sin2 θ



2

+β

Ω(ΩR+u cos θ)2
(ΩR+u cos θ)2 +u2 sin2 θ

+ c4
+Ω

(ΩR + u cos θ)2 + u2 sin2 θ

40

3c1

Ωu3 sin3 θ
(ΩR+u cos θ)4
Ωu3 sin3 θ
(ΩR+u cos θ)4

+c2

3

−

−

3 3
− u sin θ 3
3(ΩR+u cos θ)

Ωu sin θ
(ΩR+u cos θ)2

u sin θ
+ ΩR+u
cos θ + β

2(ΩR+u cos θ)Ω2
(ΩR+u cos θ)2 +u2 sin2 θ

+ (ΩR + u cos θ)





Ω2
(ΩR+u cos θ)2 +u2 sin2 θ

+c1 6 −

u3 sin3 θ
3(ΩR+u cos θ)3

2Ω2 u sin θ
(ΩR+u cos θ)3

−

3/2
(ΩR+u cos θ)2 +u2 sin2 θ

−

1
ΩR+u cos θ

+

+

−

2 3 3
− 4Ω u sin θ5

(ΩR+u cos θ)
2
Ωu sin θ
+
(ΩR+u cos θ)2

u sin θ
+ ΩR+u
cos θ + β

u2 sin2 θ
(ΩR+u cos θ)3

−

u3 sin3 θ
3(ΩR+u cos θ)3

2

2

(ΩR + u cos θ)2 + u2 sin2 θ 2 ΩR+u1 cos θ − u sin θ 3
(ΩR+u cos θ)


u3 sin3 θ
+ u sin θ +β
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)3


+

c3 −

u3 sin3 θ
3(ΩR+u cos θ)3


1
(ΩR+u cos θ)2 +u2 sin2 θ

−

u2 sin2 θ
 sin θ
3/2
2
2
2
(ΩR+u cos θ) +u sin θ

2u sin θ
(ΩR+u cos θ)2 +u2 sin2 θ

+ 0.5ρair c

2

u sin θ
+ ΩR+u
cos θ + β

(ΩR + u cos θ)2 + u2 sin2 θ

+c4
)

u 

2Ω2 u sin θ
(ΩR+u cos θ)3

u3 sin3 θ
3(ΩR+u cos θ)3

−



Ωu3 sin3 θ
(ΩR+u cos θ)4

u sin θ
+ ΩR+u
cos θ + β

(ΩR+u cos θ)

u2 sin2 θ
(ΩR+u cos θ)2 +u2 sin2 θ

4u sin θ −

Ω2 (ΩR+u cos θ)2

2 3 3
− 4Ω u sin θ5

f020 = 0.25ρair c 4c3

+c3 u sin θ

−

(ΩR + u cos θ)2 + u2 sin2 θ c2

+(ΩR + u cos θ)

3

+β

3

3

u3 sin3 θ
3(ΩR+u cos θ)3

2

Ωu sin θ
(ΩR+u cos θ)2

u sin θ
+ c1 − u sin θ 3 + ΩR+u
cos θ + β
3(ΩR+u cos θ)

+c2 −

+

u sin θ
ΩR+u cos θ

2u(ΩR+u cos θ) sin θ
(ΩR+u cos θ)2 +u2 sin2 θ

41

2

u sin θ
+ ΩR+u
cos θ + β

2

3c1

1
ΩR+u cos θ

1
ΩR+u cos θ

+c2

−

2

3 3
− u sin θ 3
3(ΩR+u cos θ)

+

u sin θ
ΩR+u cos θ

2

+β

2

− u sin θ 3
 (ΩR+u cos θ)

+(ΩR + u cos θ) 

+

u2 sin2 θ

3/2
(ΩR+u cos θ)2 +u2 sin2 θ

−

(ΩR+u cos θ)2 +u2 sin2 θ
3

u sin θ
ΩR+u cos θ

(ΩR + u cos θ)2

+(ΩR + u cos θ)



1

3 3
− u sin θ 3
3(ΩR+u cos θ)

c1

u2 sin2 θ
(ΩR+u cos θ)3

+ c2 −

+β

+ u2 sin2 θ

c1

u3 sin3 θ
3(ΩR+u cos θ)3

6

1
ΩR+u cos θ

u sin θ
+ ΩR+u
cos θ + β

−

u2 sin2 θ
(ΩR+u cos θ)3

2



−

−

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ + β

−

2
u3 sin3 θ
6u sin θ −
+ u sin θ +β 
3
ΩR+u cos θ
3(ΩR+u cos θ)


(ΩR+u cos θ)3

2c2 u sin θ
(ΩR+u cos θ)3

f002 = 0.5 3.ρair c1 c(ΩR + u cos θ) (ΩR + u cos θ)2 + u2 sin2 θ
β−

u3 sin3 θ
3(ΩR+u cos θ)3

+ρair c3 cu sin θ

f110 =
−
−

u sin θ
+ ΩR+u
cos θ

(ΩR + u cos θ)2 + u2 sin2 θ + 0.

0.5ρair cΩ
(ΩR+u cos θ)2 +u2 sin2 θ

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ + β

1
ΩR+u cos θ

3c1
2

+ c2

−

1
ΩR+u cos θ

u2 sin2 θ
(ΩR+u cos θ)3

−

u2 sin2 θ
(ΩR+u cos θ)3

(ΩR + u cos θ)2

0.5ρair cΩu sin θ
3/2
(ΩR+u cos θ)2 +u2 sin2 θ

c1

3 3
− u sin θ 3
3(ΩR+u cos θ)

+

u sin θ
ΩR+u cos θ

3

+β

+ c2 −

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ + β

(ΩR + u cos θ)2
ρair c3 cΩu sin θ

+

2 2
1
− u sin θ
ΩR+u cos θ (ΩR+u cos θ)3

u3 sin3 θ
+ u sin θ +β (ΩR+u cos θ)
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)2 +u2 sin2 θ
−

42

+

−

2
u sin θ +β +c (ΩR+u cos θ)
u3 sin3 θ
+
0.5ρair cΩ c3 −
4
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)2 +u2 sin2 θ
2
u sin θ +β +c (ΩR+u cos θ)
u3 sin3 θ
+
0.5ρair cΩu2 sin2 θ c3 −
4
3(ΩR+u cos θ)3 ΩR+u cos θ
3/2
(ΩR+u cos θ)2 +u2 sin2 θ
3Ωu2 sin2 θ
(ΩR+u cos θ)4

+ 0.5ρair c (ΩR + u cos θ)2 + u2 sin2 θ 3c1
−

u3 sin3 θ
3(ΩR+u cos θ)3

+6c1
−
+c2

1
ΩR+u cos θ

u3 sin3 θ
3(ΩR+u cos θ)3
3Ωu2 sin2 θ

u sin θ
+ ΩR+u
cos θ + β

−

u2 sin2 θ
(ΩR+u cos θ)3

−

(ΩR+u cos θ)4

+c2

u3 sin3 θ
3(ΩR+u cos θ)3
Ωu3 sin3 θ

Ωu3 sin3 θ
(ΩR+u cos θ)4

−

+ ρair c3 cu sin θ

2

Ωu3 sin3 θ
(ΩR+u cos θ)4

−

Ω
(ΩR+u cos θ)2

3c1

−

Ωu sin θ
(ΩR+u cos θ)2

−

(ΩR + u cos θ)

Ωu3 sin3 θ
(ΩR+u cos θ)4

u sin θ
+ ΩR+u
cos θ + β

(ΩR+u cos θ)4

+ ρair c3 cu sin θ

Ω
(ΩR+u cos θ)2

u sin θ
+ ΩR+u
cos θ + β

0.5ρair cu sin θ
(ΩR+u cos θ)2 +u2 sin2 θ

+

−

−

Ωu sin θ
(ΩR+u cos θ)2

2

Ωu sin θ
(ΩR+u cos θ)2

(ΩR + u cos θ)
1
ΩR+u cos θ

(ΩR + u cos θ)2 + u2 sin2 θ

−

u2 sin2 θ
(ΩR+u cos θ)3

Ωu sin θ
(ΩR+u cos θ)2
3Ωu2 sin2 θ
(ΩR+u cos θ)4

(ΩR + u cos θ)2 + u2 sin2 θ

−

Ω
(ΩR+u cos θ)2

u3 sin3 θ
u sin θ
2
2
2
+ ΩR+u
cos θ + β + ρair c3 c (ΩR + u cos θ) + u sin θ
3(ΩR+u cos θ)3
Ωu3 sin3 θ −
Ωu sin θ
u3 sin3 θ
u sin θ
−
+ ΩR+u
cos θ + β
(ΩR+u cos θ)4
(ΩR+u cos θ)2
3(ΩR+u cos θ)3
Ωu3 sin3 θ −
Ωu sin θ
u3 sin3 θ
ρair c3 cu2 sin2 θ
−
+ u sin θ +β
4
2
(ΩR+u cos θ)
(ΩR+u cos θ)
3(ΩR+u cos θ)3 ΩR+u cos θ
+
(ΩR+u cos θ)2 +u2 sin2 θ

−

+ 0.5ρair cΩ (ΩR + u cos θ)2 + u2 sin2 θ 3c1
3 3
− u sin θ 3
3(ΩR+u cos θ)

+

u sin θ
ΩR+u cos θ

2

+β

+ c2

1
ΩR+u cos θ

1
ΩR+u cos θ

43

−

−

u2 sin2 θ
(ΩR+u cos θ)3

u2 sin2 θ
(ΩR+u cos θ)3

+

3
u sin θ +β +c −
u3 sin3 θ
u3 sin3 θ
+
+ u sin θ +β
0.5ρair cΩu sin θ c1 −
2
3
ΩR+u
cos
θ
3(ΩR+u cos θ)
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)2 +u2 sin2 θ

0.5ρair cΩ(ΩR+u cos θ)2 3c1 β−

f101 =

2
u sin θ
u3 sin3 θ
+
+c2
3(ΩR+u cos θ)3 ΩR+u cos θ

(ΩR+u cos θ)2 +u2 sin2 θ
3

3

u sin θ
u sin θ
+ ΩR+u
+ 0.5ρair cΩ (ΩR + u cos θ)2 + u2 sin2 θ 3c1 β −
cos θ
3(ΩR+u cos θ)3

+ 3.ρair c1 c(ΩR + u cos θ)
β−

u3 sin3 θ
3(ΩR+u cos θ)3

(ΩR + u cos θ)2 + u2 sin2 θ

Ωu3 sin3 θ
(ΩR+u cos θ)4

−

2

Ωu sin θ
(ΩR+u cos θ)2

u sin θ
+ ΩR+u
cos θ
u3 sin3 θ
+ u sin θ
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)2 +u2 sin2 θ

ρair c3 cΩu sin θ(ΩR+u cos θ) β−

+
+ ρair c3 cu sin θ

f011 =

(ΩR + u cos θ)2 + u2 sin2 θ

1
u2 sin2 θ
−
ΩR+u cos θ
(ΩR+u cos θ)3
3
3
u sin θ
+ u sin θ
ρair c3 cu2 sin2 θ β−
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)2 +u2 sin2 θ

(ΩR + u cos θ)2 + u2 sin2 θ

u sin θ
+ ΩR+u
cos θ

+

(ΩR + u cos θ)2 + u2 sin2 θ β −

+ ρair c3 cu sin θ

f300 = 61
−

Ωu sin θ
(ΩR+u cos θ)2

(ΩR+u cos θ)2 +u2 sin2 θ

u3 sin3 θ
3(ΩR+u cos θ)3

+ ρair c3 c

−

2
u3 sin3 θ
u sin θ
0.5ρair cu sin θ(ΩR+u cos θ) 3c1 β−
+
+c2
3(ΩR+u cos θ)3 ΩR+u cos θ

+ 3.ρair c1 c(ΩR + u cos θ)
β−

Ωu3 sin3 θ
(ΩR+u cos θ)4

(ΩR + u cos θ)2 + u2 sin2 θ

0.5ρair c 6c3 u sin θ

u3 sin3 θ
3(ΩR+u cos θ)3

u3 sin3 θ
3(ΩR+u cos θ)3
1
ΩR+u cos θ

Ωu3 sin3 θ
(ΩR+u cos θ)4

u sin θ
+ ΩR+u
cos θ + β

44

−

u sin θ
+ ΩR+u
cos θ

−

u2 sin2 θ
(ΩR+u cos θ)3

Ωu sin θ
(ΩR+u cos θ)2

+ c2






Ω2

+2
−
2
−

3/2
(ΩR+u cos θ)2 +u2 sin2 θ



(ΩR + u cos θ)2 + u2 sin2 θ

+c3 u sin θ
6

−

(ΩR+u cos θ)2 +u2 sin2 θ

Ω2 (ΩR+u cos θ)2

2Ω2 u sin θ − 4Ω2 u3 sin3 θ
(ΩR+u cos θ)3
(ΩR+u cos θ)5
20Ω3 u3 sin3 θ − 6Ω3 u sin θ
(ΩR+u cos θ)6
(ΩR+u cos θ)4
u3 sin3 θ
3(ΩR+u cos θ)3
Ωu3 sin3 θ
(ΩR+u cos θ)4
u3 sin3 θ
3(ΩR+u cos θ)3

Ωu3 sin3 θ
(ΩR+u cos θ)4

u sin θ
+ ΩR+u
cos θ + β

−

Ωu sin θ
(ΩR+u cos θ)2

3c3 Ωu(ΩR+u cos θ) sin θ
(ΩR+u cos θ)2 +u2 sin2 θ

+
2

Ωu sin θ
(ΩR+u cos θ)2

−

2Ω2 u sin θ
(ΩR+u cos θ)3

+2

2 3 3
− 4Ω u sin θ5
(ΩR+u cos θ)

u sin θ
+ ΩR+u
cos θ + β


+u sin θ 
c3 −


3Ω3 (ΩR+u cos θ)3
5/2
(ΩR+u cos θ)2 +u2 sin2 θ

u3 sin3 θ
3(ΩR+u cos θ)3

c1 −

+0.5ρair c
3

3Ω3 (ΩR+u cos θ)

−

3/2
(ΩR+u cos θ)2 +u2 sin2 θ
2

u sin θ
+ ΩR+u
cos θ + β

u3 sin3 θ
3(ΩR+u cos θ)3



+ c4

u sin θ
+ ΩR+u
cos θ + β

3

3

u sin θ
+c2 − u sin θ 3 + ΩR+u
cos θ + β
3(ΩR+u cos θ)





3Ω3 (ΩR+u cos θ)3

(ΩR + u cos θ) 

5/2
(ΩR+u cos θ)2 +u2 sin2 θ

3Ω3 (ΩR+u cos θ)

−

3/2
(ΩR+u cos θ)2 +u2 sin2 θ


+3Ω 
+3 3c1
+c2




Ω2
(ΩR+u cos θ)2 +u2 sin2 θ
Ωu3 sin3 θ
(ΩR+u cos θ)4

Ωu3 sin3 θ
(ΩR+u cos θ)4

−

−

−

Ω2 (ΩR+u cos θ)2
3/2
(ΩR+u cos θ)2 +u2 sin2 θ

Ωu sin θ
(ΩR+u cos θ)2

−

u3 sin3 θ
3(ΩR+u cos θ)3



u sin θ
+ ΩR+u
cos θ + β

2(ΩR+u cos θ)Ω2
(ΩR+u cos θ)2 +u2 sin2 θ

Ωu sin θ
(ΩR+u cos θ)2

2

+ (ΩR + u cos θ)


Ω2
(ΩR+u cos θ)2 +u2 sin2 θ

−

Ω2 (ΩR+u cos θ)2
3/2
(ΩR+u cos θ)2 +u2 sin2 θ

+(ΩR + u cos θ) (ΩR + u cos θ)2 + u2 sin2 θ c2
+c1



6

Ωu3 sin3 θ
(ΩR+u cos θ)4

−

Ωu sin θ
(ΩR+u cos θ)2

3

+ 18

45



20Ω3 u3 sin3 θ
(ΩR+u cos θ)6

2Ω2 u sin θ
(ΩR+u cos θ)3

−

6Ω3 u sin θ
(ΩR+u cos θ)4

2 3 3
− 4Ω u sin θ5
(ΩR+u cos θ)

−

u3 sin3 θ
3(ΩR+u cos θ)3
20Ω3 u3 sin3 θ
(ΩR+u cos θ)6

+3

−

6Ω3 u sin θ
(ΩR+u cos θ)4

Ω(ΩR+u cos θ)2
(ΩR+u cos θ)2 +u2 sin2 θ

+3

6 −

+3

2Ω2 u sin θ
(ΩR+u cos θ)3

f030 = 61

u3 sin3 θ
3(ΩR+u cos θ)3

−

4Ω2 u3 sin3 θ
(ΩR+u cos θ)5

u2 sin2 θ
(ΩR+u cos θ)2 +u2 sin2 θ

−

3

3 3
− u sin θ 3
3(ΩR+u cos θ)

−

u sin θ
ΩR+u cos θ

(ΩR+u cos θ)2 +u2 sin2 θ
2

−

2

+6c3 ΩR+u1 cos θ − u sin θ 3
(ΩR+u cos θ)
 
u 

+β

2

−

2

+ c4

(u sin θ


3u sin θ

3/2
2
2
2
(ΩR+u cos θ) +u sin θ


1

1
(ΩR+u cos θ)2 +u2 sin2 θ

2u sin θ
(ΩR+u cos θ)2 +u2 sin2 θ

−

2

2

u3 sin3 θ
4 −
+ u sin θ +β
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)3


+3 

Ωu sin θ
(ΩR+u cos θ)2

u3 sin3 θ
+ u sin θ +β
4u sin θ −
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)3

3

3u3 sin3 θ
5/2
(ΩR+u cos θ)2 +u2 sin2 θ

+

−





(ΩR + u cos θ)2 + u2 sin2 θ

+

u2 sin2 θ
(ΩR+u cos θ)3

−

u sin θ
+ c3 − u sin θ 3 + ΩR+u
cos θ + β
3(ΩR+u cos θ)




2

(ΩR + u cos θ)2 + u2 sin2 θ

0.5ρair c c3 u sin θ

1
ΩR+u cos θ

Ωu3 sin3 θ
(ΩR+u cos θ)4

u sin θ
+ ΩR+u
cos θ + β

2 2
1
12u sin θ
− u sin θ
ΩR+u cos θ (ΩR+u cos θ)3

−
(ΩR+u cos θ)3


2

u sin θ
+ ΩR+u
cos θ + β

(ΩR + u cos θ)2 + u2 sin2 θ

+Ω



+3c3


Ωu sin θ
(ΩR+u cos θ)2

(ΩR+u cos θ)

u3 sin3 θ
3(ΩR+u cos θ)3

+c1

−

−

2 3 3
− 4Ω u sin θ5

2Ω2 u sin θ
(ΩR+u cos θ)3

c2

Ωu3 sin3 θ
(ΩR+u cos θ)4

u sin θ
+ ΩR+u
cos θ + β

u2 sin2 θ
3/2
(ΩR+u cos θ)2 +u2 sin2 θ

−

u3 sin3 θ
3(ΩR+u cos θ)3



u sin θ
+ ΩR+u
cos θ + β


u2 sin2 θ
 sin θ+
3/2
(ΩR+u cos θ)2 +u2 sin2 θ

+ 0.5ρair c (3(ΩR + u cos θ)
46








1



(ΩR+u cos θ)2 +u2 sin2 θ
1
ΩR+u cos θ

3c1

u3 sin3 θ
3(ΩR+u cos θ)3

−

−

−

u2 sin2 θ

3/2
2
2
2
(ΩR+u cos θ) +u sin θ

u2 sin2 θ
(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ + β

2

1
ΩR+u cos θ

+ c2

−

+(ΩR
+ u cos θ)





3u3 sin3 θ
5/2
(ΩR+u cos θ)2 +u2 sin2 θ

3u sin θ

3/2
2
2
2
(ΩR+u cos θ) +u sin θ

−

3

c1 −

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ + β

+c2 −

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ + β

c1

6

1
ΩR+u cos θ

1
ΩR+u cos θ

−

−

−

u2 sin2 θ
(ΩR+u cos θ)3

u3 sin3 θ
3(ΩR+u cos θ)3

3

u2 sin2 θ
(ΩR+u cos θ)3

3u(ΩR+u cos θ) sin θ
(ΩR+u cos θ)2 +u2 sin2 θ

+

u2 sin2 θ
(ΩR+u cos θ)3

−

6 −

−

c1 6

(ΩR + u cos θ)2 + u2 sin2 θ

+ (ΩR + u cos θ)
36u sin θ
(ΩR+u cos θ)3

−

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ + β



2
u sin θ +β
u3 sin3 θ
+


2c2
3(ΩR+u cos θ)3 ΩR+u cos θ
−

 (ΩR+u cos θ)3 
(ΩR+u cos θ)3

1
ΩR+u cos θ

−

u2 sin2 θ
(ΩR+u cos θ)3

2

u sin θ
+ ΩR+u
cos θ + β





2
u3 sin3 θ
u sin θ +β
6u sin θ −
+


3(ΩR+u cos θ)3 ΩR+u cos θ
 − 2c2 u sin θ 
−
3
3

(ΩR+u cos θ)
(ΩR+u cos θ) 





f003 = 21 0.5ρair c −
−

u3 sin3 θ
3(ΩR+u cos θ)3

c3

u2

−

2

(ΩR+u cos θ)2 +u2 sin2 θ

u sin θ
2
+ ΩR+u
cos θ + β sin θ +

Ωu3 sin3 θ
(ΩR+u cos θ)4

u3 sin3 θ
3(ΩR+u cos θ)3

4c3 Ωu2 (ΩR+u cos θ)

−

Ωu sin θ
(ΩR+u cos θ)2

u sin θ
+ ΩR+u
cos θ + β

Ωu3 sin3 θ
(ΩR+u cos θ)4

−

Ωu sin θ
(ΩR+u cos θ)2

1
(ΩR+u cos θ)2 +u2 sin2 θ

2

+2

sin2 θ +

47

3/2

2Ω2 u sin θ
(ΩR+u cos θ)3

2 3 3
− 4Ω u sin θ5
(ΩR+u cos θ)

Ωu3 sin3 θ −
Ωu sin θ
sin θ
(ΩR+u cos θ)4 (ΩR+u cos θ)2
(ΩR+u cos θ)2 +u2 sin2 θ
Ω
u3 sin3 θ
3Ωu2 sin2 θ −
−
+ u sin θ +β sin θ
4c3 Ωu(ΩR+u cos θ)
(ΩR+u cos θ)4 (ΩR+u cos θ)2
3(ΩR+u cos θ)3 ΩR+u cos θ

2 2
1
4c3 Ωu(ΩR+u cos θ)
− u sin θ
ΩR+u cos θ (ΩR+u cos θ)3

+

(ΩR+u cos θ)2 +u2 sin2 θ
1
ΩR+u cos θ

+2c3 u


u2 sin2 θ
(ΩR+u cos θ)3

−

Ω2



(ΩR+u cos θ)2 +u2 sin2 θ

−

2Ω2 u sin θ
(ΩR+u cos θ)3
Ωu3 sin3 θ
(ΩR+u cos θ)4

3/2
(ΩR+u cos θ)2 +u2 sin2 θ

3

2 3 3
− 4Ω u sin θ5

+4

Ωu sin θ
(ΩR+u cos θ)2

+2

(ΩR+u cos θ)

−

3

u sin θ
− u sin θ 3 + ΩR+u
cos θ + β
3(ΩR+u cos θ)


+u 
c3

3Ω2 u(ΩR+u cos θ)2 sin θ
5/2
(ΩR+u cos θ)2 +u2 sin2 θ
3 3
− u sin θ 3
3(ΩR+u cos θ)

+

u sin θ
+ ΩR+u
cos θ + β


Ω2 (ΩR+u cos θ)2

−

1
ΩR+u cos θ

(ΩR + u cos θ)2 + u2 sin2 θ 2

+c3 u

u3 sin3 θ
3(ΩR+u cos θ)3

3Ωu2 sin2 θ
(ΩR+u cos θ)4
2Ω2
(ΩR+u cos θ)3

u2 sin2 θ
(ΩR+u cos θ)3

−
−

Ω
(ΩR+u cos θ)2
12Ω2 u2 sin2 θ
(ΩR+u cos θ)5

sin θ

Ω2 u sin θ

3/2
2
2
2
(ΩR+u cos θ) +u sin θ

−

u sin θ
ΩR+u cos θ

−

 sin θ

2

+β

+ c4

sin θ

Ωu sin θ
Ωu3 sin3 θ −
u3 sin3 θ
+ u sin θ +β
4c3 Ω(ΩR+u cos θ)
−
(ΩR+u cos θ)4 (ΩR+u cos θ)2
3(ΩR+u cos θ)3 ΩR+u cos θ
+
(ΩR+u cos θ)2 +u2 sin2 θ

+c3
2

(ΩR + u cos θ)2 + u2 sin2 θ
Ωu3 sin3 θ
(ΩR+u cos θ)4

−

Ωu sin θ
(ΩR+u cos θ)2

2

2 3

3

2Ω u sin θ − 4Ω u sin θ
+2
(ΩR+u cos θ)3
(ΩR+u cos θ)5


+

Ω2

u3 sin3 θ
3(ΩR+u cos θ)3

+c4 )) + 0.5ρair c 2

+

−

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ + β


Ω2 (ΩR+u cos θ)2

3/2
(ΩR+u cos θ)2 +u2 sin2 θ
2
u sin θ
ΩR+u cos θ + β

(ΩR+u cos θ)2 +u2 sin2 θ

c3 −

−

2

Ω(ΩR+u cos θ)2
(ΩR+u cos θ)2 +u2 sin2 θ

+Ω

48



(ΩR + u cos θ)2 + u2 sin2 θ

3Ωu2 sin2 θ
(ΩR+u cos θ)4

3c1

1
ΩR+u cos θ

6c1
−

Ωu3 sin3 θ
(ΩR+u cos θ)4

u2 sin2 θ
(ΩR+u cos θ)3

−

u3 sin3 θ
3(ΩR+u cos θ)3
3Ωu2 sin2 θ

−

3 3
− u sin θ 3
3(ΩR+u cos θ)

Ω
(ΩR+u cos θ)2

−

+

u sin θ
ΩR+u cos θ

+c2

Ωu3 sin3 θ
(ΩR+u cos θ)4
Ωu3 sin3 θ
(ΩR+u cos θ)4

+c2 −

+

u sin θ
+ ΩR+u
cos θ + β

Ωu sin θ
(ΩR+u cos θ)2 +u2 sin2 θ

3c1

+β

Ωu sin θ
(ΩR+u cos θ)2

Ω
−
+c2
(ΩR+u cos θ)4
(ΩR+u cos θ)2


+2 

2

−

−

u3 sin3 θ
3(ΩR+u cos θ)3


Ωu(ΩR+u cos θ)2 sin θ

−

3/2

(ΩR+u cos θ)2 +u2 sin2 θ

Ωu sin θ
(ΩR+u cos θ)2

−

Ωu sin θ
(ΩR+u cos θ)2

+

u sin θ
+ ΩR+u
cos θ + β


u3 sin3 θ
3(ΩR+u cos θ)3

c1



2

u sin θ
+ ΩR+u
cos θ + β

3 3
− u sin θ 3
3(ΩR+u cos θ)

+

u sin θ
ΩR+u cos θ

3

+β


3Ω2 u(ΩR+u cos θ)2 sin θ
5/2
(ΩR+u cos θ)2 +u2 sin2 θ

(ΩR + u cos θ) 



2Ω2 u(ΩR+u cos θ) sin θ

Ω2 u sin θ
−

3/2
3/2
2
2
2
2
2
2
(ΩR+u cos θ) +u sin θ
(ΩR+u cos θ) +u sin θ

−

+ 3c1
+c2


1
ΩR+u cos θ

1
ΩR+u cos θ

−

u2 sin2 θ
(ΩR+u cos θ)3

2

u sin θ
+ ΩR+u
cos θ + β

+ (ΩR + u cos θ)


(ΩR+u cos θ)2 +u2 sin2 θ

−

Ω2 (ΩR+u cos θ)2
3/2
(ΩR+u cos θ)2 +u2 sin2 θ

(ΩR + u cos θ)2 + u2 sin2 θ c2
+c1 6

−

3(ΩR+u cos θ)3

+3

2Ω2
(ΩR+u cos θ)3

+6

1
ΩR+u cos θ

−

2Ω2
(ΩR+u cos θ)3

u2 sin2 θ
(ΩR+u cos θ)3
2
Ωu sin θ
+ 12
(ΩR+u cos θ)2

1
ΩR+u cos θ

Ωu3 sin3 θ
(ΩR+u cos θ)4
u3 sin3 θ

−

u3 sin3 θ
3(ΩR+u cos θ)3

2(ΩR+u cos θ)Ω2
(ΩR+u cos θ)2 +u2 sin2 θ

u2 sin2 θ
(ΩR+u cos θ)3

−

Ω2



−

 + (ΩR + u cos θ)

2 2 2
− 12Ω u sin θ5
(ΩR+u cos θ)

−

u sin θ
+ ΩR+u
cos θ + β
2 2 2
− 12Ω u sin θ5

(ΩR+u cos θ)
u2 sin2 θ
(ΩR+u cos θ)3

3Ωu2 sin2 θ −
Ω
(ΩR+u cos θ)4
(ΩR+u cos θ)2
Ωu3 sin3 θ −
Ωu sin θ
(ΩR+u cos θ)4
(ΩR+u cos θ)2

u3 sin3 θ
u sin θ
+ ΩR+u
cos θ
3(ΩR+u cos θ)3
2Ω2 u sin θ − 4Ω2 u3 sin3 θ
(ΩR+u cos θ)3
(ΩR+u cos θ)5

−

49

2

+β

u sin θ
u3 sin3 θ
+ ΩR+u
cos θ
3(ΩR+u cos θ)3
1
(ΩR+u cos θ)2 +u2 sin2 θ

−
+

2Ω2 u sin θ
(ΩR+u cos θ)3

u(ΩR + u cos θ) sin θ c2
+c1 6

3 3
− u sin θ 3
3(ΩR+u cos θ)

2Ω2 u sin θ
(ΩR+u cos θ)3

+3

−



+


2



u sin θ
ΩR+u cos θ

4Ω2 u3 sin3 θ
(ΩR+u cos θ)5



f210 = 21 0.5ρair c −
Ωu3 sin3 θ
(ΩR+u cos θ)4

−

+β

2 3 3
− 4Ω u sin θ5
(ΩR+u cos θ)

Ωu3 sin3 θ
(ΩR+u cos θ)4

+β

3 3
− u sin θ 3
3(ΩR+u cos θ)

+

−

Ωu sin θ
(ΩR+u cos θ)2

u sin θ
ΩR+u cos θ

2

+β

4c3 Ωu2 (ΩR+u cos θ)
3/2
(ΩR+u cos θ)2 +u2 sin2 θ

Ωu sin θ
(ΩR+u cos θ)2

−

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
2
2
+ ΩR+u
cos θ + β sin θc3 u

2
Ωu sin θ
Ωu3 sin3 θ −
(ΩR+u cos θ)4 (ΩR+u cos θ)2
(ΩR+u cos θ)2 +u2 sin2 θ





2Ω2 u sin θ − 4Ω2 u3 sin3 θ
2
 (ΩR+u cos θ)3 (ΩR+u cos θ)5

u3 sin3 θ
+ u sin θ +β sin2 θ
−

3(ΩR+u cos θ)3 ΩR+u cos θ

(ΩR+u cos θ)2 +u2 sin2 θ

+2 

4c3 Ωu(ΩR+u cos θ)

+
4c3 Ωu(ΩR+u cos θ)

+
+2c3 u



2

1
ΩR+u cos θ

2 2
1
− u sin θ
ΩR+u cos θ (ΩR+u cos θ)3

Ωu3 sin3 θ −
Ωu sin θ
(ΩR+u cos θ)4 (ΩR+u cos θ)2
(ΩR+u cos θ)2 +u2 sin2 θ

3Ωu2 sin2 θ −
Ω
u3 sin3 θ
+ u sin θ +β sin θ
−
4
2
(ΩR+u cos θ)
(ΩR+u cos θ)
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)2 +u2 sin2 θ

−

u2 sin2 θ
(ΩR+u cos θ)3

Ω2
(ΩR+u cos θ)2 +u2 sin2 θ

+c3 u

sin θ

−

−

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ + β


Ω2 (ΩR+u cos θ)2
3/2
(ΩR+u cos θ)2 +u2 sin2 θ

(ΩR + u cos θ)2 + u2 sin2 θ 2

2Ω2 u sin θ − 4Ω2 u3 sin3 θ
(ΩR+u cos θ)3
(ΩR+u cos θ)5
3Ωu2 sin2 θ −
Ω
+4
(ΩR+u cos θ)4
(ΩR+u cos θ)2

1
ΩR+u cos θ

−

Ωu3 sin3 θ
(ΩR+u cos θ)4

50

 sin θ

u2 sin2 θ
(ΩR+u cos θ)3

−

Ωu sin θ
(ΩR+u cos θ)2

2Ω2
(ΩR+u cos θ)3

+2

2 2 2
− 12Ω u sin θ5

−

(ΩR+u cos θ)

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ + β



+u 

3Ω2 u(ΩR+u cos θ)2 sin θ
5/2
(ΩR+u cos θ)2 +u2 sin2 θ

c3 −

u3 sin3 θ
3(ΩR+u cos θ)3

sin θ

Ω2 u sin θ

3/2
(ΩR+u cos θ)2 +u2 sin2 θ

−

2

u sin θ
+ ΩR+u
cos θ + β

+ c4

sin θ

Ωu sin θ
u3 sin3 θ
Ωu3 sin3 θ −
−
+ u sin θ +β
4c3 Ω(ΩR+u cos θ)
4
2
(ΩR+u cos θ)
(ΩR+u cos θ)
3(ΩR+u cos θ)3 ΩR+u cos θ
+
(ΩR+u cos θ)2 +u2 sin2 θ

(ΩR + u cos θ)2

+c3

+ u2 sin2 θ

2

2 3

2

3

2Ω u sin θ − 4Ω u sin θ
+2
(ΩR+u cos θ)3
(ΩR+u cos θ)5


+

Ω2

u3 sin3 θ
3(ΩR+u cos θ)3

+

+6c1
−

−

1
ΩR+u cos θ

u2 sin2 θ
(ΩR+u cos θ)3

u3 sin3 θ
3(ΩR+u cos θ)3
3Ωu2 sin2 θ

u sin θ
+ ΩR+u
cos θ + β


Ω
(ΩR+u cos θ)2

−

(ΩR + u cos θ)2 + u2 sin2 θ

+Ω

u3 sin3 θ
3(ΩR+u cos θ)3

Ωu3 sin3 θ
(ΩR+u cos θ)4



−

u sin θ
+ ΩR+u
cos θ + β

3c1
+c2

u sin θ
+ ΩR+u
cos θ + β

Ωu sin θ
(ΩR+u cos θ)2 +u2 sin2 θ
Ωu3 sin3 θ
(ΩR+u cos θ)4
Ωu3 sin3 θ
(ΩR+u cos θ)4

+c2 −

2

Ωu sin θ
(ΩR+u cos θ)2

Ω
+c2
−
(ΩR+u cos θ)4
(ΩR+u cos θ)2


+2 

2

Ω2 (ΩR+u cos θ)2

−

3Ωu2 sin2 θ
(ΩR+u cos θ)4

−

u3 sin3 θ
3(ΩR+u cos θ)3

Ω(ΩR+u cos θ)2
(ΩR+u cos θ)2 +u2 sin2 θ

+c4 )) + 0.5ρair c 2
3c1

−

−

Ωu sin θ
(ΩR+u cos θ)2

3/2
(ΩR+u cos θ)2 +u2 sin2 θ
2
u sin θ
ΩR+u cos θ + β

(ΩR+u cos θ)2 +u2 sin2 θ

c3 −

Ωu3 sin3 θ
(ΩR+u cos θ)4

−

−

u3 sin3 θ
3(ΩR+u cos θ)3

−


Ωu(ΩR+u cos θ)2 sin θ

3/2
(ΩR+u cos θ)2 +u2 sin2 θ

Ωu sin θ
(ΩR+u cos θ)2

Ωu sin θ
(ΩR+u cos θ)2

−
+

u sin θ
+ ΩR+u
cos θ + β

u3 sin3 θ
3(ΩR+u cos θ)3

c1 −

u sin θ
+ ΩR+u
cos θ + β

u3 sin3 θ
3(ΩR+u cos θ)3

((ΩR + u cos θ)

51


2

u sin θ
+ ΩR+u
cos θ + β

3




3Ω2 u(ΩR+u cos θ)2 sin θ
5/2
(ΩR+u cos θ)2 +u2 sin2 θ



−

Ω2 u sin θ

3/2
2
2
2
(ΩR+u cos θ) +u sin θ


2Ω2 u(ΩR+u cos θ) sin θ

−

3/2
(ΩR+u cos θ)2 +u2 sin2 θ

+ 3c1
+c2


1
ΩR+u cos θ

1
ΩR+u cos θ

−

−

u2 sin2 θ
(ΩR+u cos θ)3

−

−

Ω2 (ΩR+u cos θ)2
3/2
(ΩR+u cos θ)2 +u2 sin2 θ
2Ω2
(ΩR+u cos θ)3

(ΩR + u cos θ)2 + u2 sin2 θ c2

2 2
+c1 6 ΩR+u1 cos θ − u sin θ 3
(ΩR+u cos θ)
2
Ωu sin θ
Ωu3 sin3 θ −
+ 12
(ΩR+u cos θ)4
(ΩR+u cos θ)2

+3

u sin θ
+ ΩR+u
cos θ + β

2

+ (ΩR + u cos θ)


(ΩR+u cos θ)2 +u2 sin2 θ

−

u3 sin3 θ
3(ΩR+u cos θ)3

2(ΩR+u cos θ)Ω2
(ΩR+u cos θ)2 +u2 sin2 θ

u2 sin2 θ
(ΩR+u cos θ)3

Ω2





u sin θ
u3 sin3 θ
+ ΩR+u
cos θ + β
3(ΩR+u cos θ)3
2Ω2
12Ω2 u2 sin2 θ
−
(ΩR+u cos θ)3
(ΩR+u cos θ)5

3 3
u sin θ
− u sin θ 3 + ΩR+u
cos θ
3(ΩR+u cos θ)
2Ω2 u sin θ − 4Ω2 u3 sin3 θ
(ΩR+u cos θ)3
(ΩR+u cos θ)5
3 3
u sin θ
− u sin θ 3 + ΩR+u
cos θ
3(ΩR+u cos θ)

2 2 2
− 12Ω u sin θ5
(ΩR+u cos θ)

3Ωu2 sin2 θ −
Ω
(ΩR+u cos θ)4
(ΩR+u cos θ)2
Ωu3 sin3 θ −
Ωu sin θ
(ΩR+u cos θ)4
(ΩR+u cos θ)2

2

+β

 + (ΩR + u cos θ)

+6

1
ΩR+u cos θ

−

u2 sin2 θ
(ΩR+u cos θ)3

+β

u(ΩR+u cos θ) sin θ
(ΩR+u cos θ)2 +u2 sin2 θ

+
c2

2Ω2 u sin θ
(ΩR+u cos θ)3

Ωu3 sin3 θ
(ΩR+u cos θ)4

−

−

u3 sin3 θ
3(ΩR+u cos θ)3

2 3 3
− 4Ω u sin θ5

(ΩR+u cos θ)
2
Ωu sin θ
(ΩR+u cos θ)2

+ c1 6 −

+3

u sin θ
+ ΩR+u
cos θ + β

u3 sin3 θ
3(ΩR+u cos θ)3

2Ω2 u sin θ
(ΩR+u cos θ)3

2

52

u sin θ
+ ΩR+u
cos θ + β

2 3 3
− 4Ω u sin θ5
(ΩR+u cos θ)

f201 = 21

2Ω2 (ΩR+u cos θ)
(ΩR+u cos θ)2 +u2 sin2 θ

0.5ρair c




Ω2
(ΩR+u cos θ)2 +u2 sin2 θ

β−

u3 sin3 θ
3(ΩR+u cos θ)3

+

+12c1

u3 sin3 θ
3(ΩR+u cos θ)3

−

+Ω

3

Ω2
(ΩR+u cos θ)2 +u2 sin2 θ
u3 sin3 θ
3(ΩR+u cos θ)3

+2c3 u sin θ

Ωu3 sin3 θ
(ΩR+u cos θ)4

+6

−

−

2

3/2
(ΩR+u cos θ)2 +u2 sin2 θ



u sin θ
+ ΩR+u
cos θ
2Ω2 u sin θ
(ΩR+u cos θ)3

(ΩR + u cos θ)2 + u2 sin2 θ

2 3

3

− 4Ω u sin θ5
 (ΩR+u cos θ)



1

0.5ρair c (ΩR + u cos θ)


(ΩR+u cos θ)2 +u2 sin2 θ

u2 sin2 θ

3/2
2
2
2
(ΩR+u cos θ) +u sin θ

3c1 β −

+

Ωu sin θ
(ΩR+u cos θ)2

(ΩR+u cos θ)

Ω2 (ΩR+u cos θ)2

Ωu3 sin3 θ −
Ωu sin θ
(ΩR+u cos θ)4 (ΩR+u cos θ)2
(ΩR+u cos θ)2 +u2 sin2 θ

f021 = 21

−

2 3 3
− 4Ω u sin θ5

+ 0.5ρair c (2c3 u sin θ


4c3 Ωu sin θ(ΩR+u cos θ)

+

2Ω2 u sin θ
(ΩR+u cos θ)3

(ΩR + u cos θ)2 + u2 sin2 θ

u sin θ
u sin θ
β−
+ ΩR+u
cos θ
3(ΩR+u cos θ)3


β−

 (3c1

Ωu sin θ
(ΩR+u cos θ)2

Ω(ΩR+u cos θ)2
(ΩR+u cos θ)2 +u2 sin2 θ



3/2

(ΩR + u cos θ)2 + u2 sin2 θ 6

u sin θ
+ ΩR+u
cos θ

Ωu3 sin3 θ
(ΩR+u cos θ)4

3

Ω2 (ΩR+u cos θ)2

−

(ΩR+u cos θ)2 +u2 sin2 θ
2
u sin θ
ΩR+u cos θ

+c2 ) + c1 (ΩR + u cos θ)
β−

+ (ΩR + u cos θ)


u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

2

+ c2

2 2
1
12c1 u sin θ(ΩR+u cos θ)
− u sin θ
ΩR+u cos θ (ΩR+u cos θ)3

β−

(ΩR+u cos θ)2 +u2 sin2 θ

53

u3 sin3 θ
+ u sin θ
3(ΩR+u cos θ)3 ΩR+u cos θ

(ΩR + u cos θ)2

+c1 (ΩR + u cos θ)
12u sin θ β−

+ u2 sin2 θ

1
ΩR+u cos θ

6

2

−



u3 sin3 θ
+ u sin θ
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)3




1

+0.5ρair c 2c3 u sin θ

(ΩR+u cos θ)2 +u2 sin2 θ


−

−

u2 sin2 θ
(ΩR+u cos θ)3



u2 sin2 θ
+
3/2
(ΩR+u cos θ)2 +u2 sin2 θ

u sin θ
u3 sin3 θ
+ ΩR+u
cos θ
3(ΩR+u cos θ)3
4c u2 sin2 θ (ΩR+u cos θ)2 +u2 sin2 θ
− 3
(ΩR+u cos θ)3

2u sin θ

(ΩR+u cos θ)2 +u2 sin2 θ

β−

u2 sin2 θ
(ΩR+u cos θ)2 +u2 sin2 θ

f120 = 12
3

−

1
ΩR+u cos θ

−

u2 sin2 θ
(ΩR+u cos θ)3

3

u sin θ
− u sin θ 3 + ΩR+u
cos θ + β
3(ΩR+u cos θ)




Ω(ΩR+u cos θ)
(ΩR+u cos θ)2 +u2 sin2 θ

+4c3

1
ΩR+u cos θ

Ωu3 sin3 θ
(ΩR+u cos θ)4

+4c3

u2 sin2 θ
(ΩR+u cos θ)3

(ΩR + u cos θ)2 + u2 sin2 θ

+

0.5ρair c 4c3

1
ΩR+u cos θ

+ 4c3

−

−

−

3/2
(ΩR+u cos θ)2 +u2 sin2 θ



u2 sin2 θ
(ΩR+u cos θ)3

Ωu sin θ
(ΩR+u cos θ)2

3Ωu2 sin2 θ
(ΩR+u cos θ)4



Ωu2 (ΩR+u cos θ) sin2 θ

−

Ω
(ΩR+u cos θ)2

u2 sin2 θ
(ΩR+u cos θ)2 +u2 sin2 θ

−

u3 sin3 θ
3(ΩR+u cos θ)3

+

(ΩR + u cos θ)2 + u2 sin2 θ

u sin θ
+ ΩR+u
cos θ + β

u2 sin2 θ
(ΩR+u cos θ)2 +u2 sin2 θ

+ (ΩR + u cos θ)2 + u2 sin2 θ + c3 u sin θ
4

3Ωu2 sin2 θ
(ΩR+u cos θ)4

−

Ω
(ΩR+u cos θ)2

(ΩR + u cos θ)2 + u2 sin2 θ

1
ΩR+u cos θ

54

−

u2 sin2 θ
(ΩR+u cos θ)3

4u sin θ

−

Ωu sin θ
Ωu3 sin3 θ −
(ΩR+u cos θ)4 (ΩR+u cos θ)2
(ΩR+u cos θ)3

12Ωu sin θ −

+

u3 sin3 θ
+ u sin θ +β
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)4





c Ωu(ΩR+u cos θ) sin θ

+ 3
(ΩR+u cos θ)2 +u2 sin2 θ

1
ΩR+u cos θ


2

−

3

u2 sin2 θ
(ΩR+u cos θ)3

2

−

u3 sin3 θ
+ u sin θ +β
4u sin θ −
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)3

3

u sin θ
+ c3 − u sin θ 3 + ΩR+u
cos θ + β
3(ΩR+u cos θ)

Ω(ΩR+u cos θ)

−

(ΩR+u cos θ)2 +u2 sin2 θ

3/2



2

+ c4







u sin θ 

3Ωu2 (ΩR+u cos θ) sin2 θ
5/2
(ΩR+u cos θ)2 +u2 sin2 θ



2Ωu(ΩR+u cos θ) sin θ

−

(ΩR+u cos θ)2 +u2 sin2 θ

Ωu3 sin3 θ −
Ωu sin θ
(ΩR+u cos θ)4
(ΩR+u cos θ)2
3 3
u sin θ
− u sin θ 3 + ΩR+u
cos θ + β
3(ΩR+u cos θ)

3/2



+2c3

(u



1



(ΩR+u cos θ)2 +u2 sin2 θ

−

u2 sin2 θ

3/2
(ΩR+u cos θ)2 +u2 sin2 θ


sin θ +

3c1

2u sin θ
(ΩR+u cos θ)2 +u2 sin2 θ
1
ΩR+u cos θ

−

+ 0.5ρair c −

u2 sin2 θ
(ΩR+u cos θ)3

1
u2 sin2 θ
ΩR+u cos θ − (ΩR+u cos θ)3
Ω
(ΩR+u cos θ)2 +u2 sin2 θ

+c2
+

c1 6

1
ΩR+u cos θ

−

u2 sin2 θ
(ΩR+u cos θ)3

2Ωu sin θ
3/2
(ΩR+u cos θ)2 +u2 sin2 θ

3 3
− u sin θ 3
3(ΩR+u cos θ)

+

u sin θ
ΩR+u cos θ

2

+β

(ΩR + u cos θ)2

2

−

u3 sin3 θ
3(ΩR+u cos θ)3



u sin θ
+ ΩR+u
cos θ + β


2
u3 sin3 θ
u sin θ +β
6u sin θ −
+


3(ΩR+u cos θ)3 ΩR+u cos θ
 − 2c2 u sin θ  (ΩR + u cos θ)2
−
3
3
 (ΩR+u cos θ) 
(ΩR+u cos θ)

+

2u sin θ
(ΩR+u cos θ)2 +u2 sin2 θ

3c1

3Ωu2 sin2 θ
(ΩR+u cos θ)4

−

Ω
(ΩR+u cos θ)2

3 3
− u sin θ 3
3(ΩR+u cos θ)

55

+

u sin θ
ΩR+u cos θ

2

+β

1
ΩR+u cos θ

+6c1

u3 sin3 θ
3(ΩR+u cos θ)3

−

−

Ωu3 sin3 θ
(ΩR+u cos θ)4

−

Ωu sin θ
(ΩR+u cos θ)2

u sin θ
+ ΩR+u
cos θ + β

3Ωu2 sin2 θ
(ΩR+u cos θ)4

+c2

u2 sin2 θ
(ΩR+u cos θ)3

−

Ω
(ΩR+u cos θ)2

1

(ΩR + u cos θ) +

(ΩR+u cos θ)2 +u2 sin2 θ


−

u2 sin2 θ

3/2
2
2
2
(ΩR+u cos θ) +u sin θ
3

3

u sin θ
− u sin θ 3 + ΩR+u
cos θ + β
3(ΩR+u cos θ)


+

3Ωu2 (ΩR+u cos θ) sin2 θ
5/2
(ΩR+u cos θ)2 +u2 sin2 θ

Ωu3 sin3 θ
(ΩR+u cos θ)4

3c1
2

−

Ωu sin θ
(ΩR+u cos θ)2

(ΩR + u cos θ)


Ω(ΩR+u cos θ)

−

3/2
(ΩR+u cos θ)2 +u2 sin2 θ
3

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ + β

+c2 −

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ + β
Ωu3 sin3 θ
(ΩR+u cos θ)4

+ c1 6

Ωu sin θ
(ΩR+u cos θ)2

Ωu3 sin3 θ
(ΩR+u cos θ)4

+ c2

c1 −

6c2 Ωu sin θ
(ΩR+u cos θ)4

−

−

3Ωu2 sin2 θ −
Ω
(ΩR+u cos θ)4
(ΩR+u cos θ)2
1
u2 sin2 θ
ΩR+u cos θ − (ΩR+u cos θ)3

+12



(ΩR + u cos θ) +
Ωu sin θ
(ΩR+u cos θ)2

−

u3 sin3 θ
3(ΩR+u cos θ)3

(ΩR + u cos θ)2 + u2 sin2 θ

1
ΩR+u cos θ

−

u2 sin2 θ
(ΩR+u cos θ)3

2

u sin θ
+ ΩR+u
cos θ + β

2
u sin θ +β
u3 sin3 θ
+
3(ΩR+u cos θ)3 ΩR+u cos θ
+
(ΩR+u cos θ)4
3
Ωu sin3 θ −
Ωu sin θ
u3 sin3 θ
12u sin θ
−
+ u sin θ +β
4
2
(ΩR+u cos θ)
(ΩR+u cos θ)
3(ΩR+u cos θ)3 ΩR+u cos θ
−
(ΩR+u cos θ)3
18Ωu sin θ −



 (ΩR + u cos θ)

2Ωu sin θ
(ΩR+u cos θ)2 +u2 sin2 θ

+





+Ω 

u2 sin2 θ

3/2
2
2
2
(ΩR+u cos θ) +u sin θ

c1

1
(ΩR+u cos θ)2 +u2 sin2 θ
3 3
− u sin θ 3
3(ΩR+u cos θ)

+c2 −

u3 sin3 θ
3(ΩR+u cos θ)3

1
ΩR+u cos θ

−

+

−

u sin θ
ΩR+u cos θ

3

+β

u sin θ
+ ΩR+u
cos θ + β

u2 sin2 θ
(ΩR+u cos θ)3

2

−

+Ω

u3 sin3 θ
3(ΩR+u cos θ)3

56

(ΩR + u cos θ)2 + u2 sin2 θ (c1 (6
u sin θ
+ ΩR+u
cos θ + β





2
u sin θ +β
u3 sin3 θ
+
6u sin θ −


3(ΩR+u cos θ)3 ΩR+u cos θ
 − 2c2 u sin θ 
−
3
3
 (ΩR+u cos θ) 
(ΩR+u cos θ)



f102 = 21 

u3 sin3 θ
+ u sin θ
3.ρair c1 cΩ(ΩR+u cos θ)2 β−
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)2 +u2 sin2 θ

(ΩR + u cos θ)2 + u2 sin2 θ β −
(ΩR + u cos θ)2 + u2 sin2 θ



f012 = 21 

u3 sin3 θ
3(ΩR+u cos θ)3

Ωu3 sin3 θ
(ΩR+u cos θ)4

−

u sin θ
+ ΩR+u
cos θ

Ωu sin θ
(ΩR+u cos θ)2

(ΩR + u cos θ)2 + u2 sin2 θ

1
ΩR+u cos θ

+
+

ρ c cΩu sin θ(ΩR+u cos θ)
+ air 3
(ΩR+u cos θ)2 +u2 sin2 θ

1
ΩR+u cos θ

−

u2 sin2 θ
(ΩR+u cos θ)3

+

ρair c3 cu2 sin2 θ
(ΩR+u cos θ)2 +u2 sin2 θ

(ΩR + u cos θ)2 + u2 sin2 θ

f111 =

−

+ 3.ρair c1 c(ΩR + u cos θ)

u3 sin3 θ
+ u sin θ
3.ρair c1 cu sin θ(ΩR+u cos θ) β−
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)2 +u2 sin2 θ

+3.ρair c1 c(ΩR + u cos θ)

+ρair c3 c

+ 3.ρair c1 cΩ

3.ρair c1 cΩ
(ΩR+u cos θ)2 +u2 sin2 θ

−

u2 sin2 θ
(ΩR+u cos θ)3

−

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
2
+ ΩR+u
cos θ + β (ΩR + u cos θ)

2
u3 sin3 θ
0.5ρair cΩu sin θ 3c1 −
+ u sin θ +β +c2 (ΩR+u cos θ)2
3
ΩR+u cos θ
3(ΩR+u cos θ)
3/2
(ΩR+u cos θ)2 +u2 sin2 θ
2 2
1
ρair c3 cΩu sin θ
− u sin θ
(ΩR+u cos θ)
ΩR+u cos θ (ΩR+u cos θ)3
(ΩR+u cos θ)2 +u2 sin2 θ
3.ρair c1 cu sin θ
(ΩR+u cos θ)2 +u2 sin2 θ
Ωu3 sin3 θ
(ΩR+u cos θ)4

−

Ωu sin θ
(ΩR+u cos θ)2

−

u3 sin3 θ
3(ΩR+u cos θ)3

57

u sin θ
+ ΩR+u
cos θ + β (ΩR + u cos θ)

ρair c3 cΩ −

+

−

u3 sin3 θ
+ u sin θ +β (ΩR+u cos θ)
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)2 +u2 sin2 θ
u3 sin3 θ
+ u sin θ +β (ΩR+u cos θ)
3(ΩR+u cos θ)3 ΩR+u cos θ
3/2
(ΩR+u cos θ)2 +u2 sin2 θ

ρair c3 cΩu2 sin2 θ −

1
ΩR+u cos θ

+ 0.5ρair c (ΩR + u cos θ)2 + u2 sin2 θ 6c1
Ωu sin θ
Ωu3 sin3 θ −
(ΩR+u cos θ)4
(ΩR+u cos θ)2
3Ωu2 sin2 θ −
Ω
+6c1
(ΩR+u cos θ)4
(ΩR+u cos θ)2

(ΩR + u cos θ) + ρair c3 cu sin θ

−

u3 sin3 θ
3(ΩR+u cos θ)3

−

u2 sin2 θ
(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ + β

(ΩR + u cos θ)2 + u2 sin2 θ

3Ωu2 sin2 θ
(ΩR+u cos θ)4

−

Ω
(ΩR+u cos θ)2

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

+ ρair c3 c
Ωu3 sin3 θ
(ΩR+u cos θ)4

(ΩR + u cos θ)2 + u2 sin2 θ

−

+

Ωu3 sin3 θ −
Ωu sin θ
(ΩR+u cos θ)4 (ΩR+u cos θ)2
(ΩR+u cos θ)2 +u2 sin2 θ

+ 3.ρair c1 cΩ

(ΩR + u cos θ)2 + u2 sin2 θ

Ωu sin θ
(ΩR+u cos θ)2

ρair c3 cu2 sin2 θ

−

+

u3 sin3 θ
3(ΩR+u cos θ)3

1
ΩR+u cos θ

−

u2 sin2 θ
(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ + β

2
u3 sin3 θ
u sin θ +β +c
0.5ρair cΩu sin θ 3c1 −
+
2
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)2 +u2 sin2 θ

M000 = 0.5ρair clc (ΩR + u cos θ)2 + u2 sin2 θ
c1 −

c2
2

+ c3 −

+(c2 + c4 ) β −

c4
6

β−

u3 sin3 θ
3(ΩR+u cos θ)3

u3 sin3 θ
3(ΩR+u cos θ)3

+

u sin θ
ΩR+u cos θ

u sin θ
+ ΩR+u
cos θ

M100 = ρair clc Ω(ΩR+u cos θ)

c
c
c1 − 22 + c3 − 64

u3 sin3 θ
u sin θ
+ ΩR+u
cos θ
3(ΩR+u cos θ)3
2
2
2
(ΩR + u cos θ) + u sin θ

+(c2 + c4 ) β −
+ 0.5ρair clc

3

58

β−

3

3 c1 −

c2
2

+(c2 + c4 )

+ c3 −

Ωu3 sin3 θ
(ΩR+u cos θ)4

c4
6

Ωu3 sin3 θ
(ΩR+u cos θ)4

M010 = ρair clc u sin θ
+(c2 + c4 ) β −

1
ΩR+u cos θ

c
c
c1 − 22 + c3 − 64

β−

u sin θ
+ ΩR+u
cos θ

1
ΩR+u cos θ

−

u3 sin3 θ
3(ΩR+u cos θ)3

−

u sin θ
+ ΩR+u
cos θ

u3 sin3 θ
3(ΩR+u cos θ)3
c2
c4
2 + c3 − 6

+(c2 + c4 ) β −

3

2Ω2 u sin θ
(ΩR+u cos θ)3

Ωu sin θ
(ΩR+u cos θ)2
2Ω2 u sin θ
(ΩR+u cos θ)3

−

β−

u2 sin2 θ
(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

u3 sin3 θ
3(ΩR+u cos θ)3

β−

u sin θ
+ ΩR+u
cos θ

2

c
c
3 c1 − 22 + c3 − 64

c
c
c1 − 22 + c3 − 64

u sin θ
+ ΩR+u
cos θ

β−

u3 sin3 θ
3(ΩR+u cos θ)3

β−

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

+

u sin θ
ΩR+u cos θ

2

+6

+ (c2 + c4 )

(ΩR+u cos θ)

u sin θ
+ ΩR+u
cos θ

3

+ lc (ΩR + u cos θ)2 + u2 sin2 θ

2 3 3
− 4Ω u sin θ5

u3 sin3 θ
3(ΩR+u cos θ)3

3

+ c2 + c4

Ωu3 sin3 θ
(ΩR+u cos θ)4

c
c
+4lc Ω(ΩR + u cos θ) 3 c1 − 22 + c3 − 64

β−

u3 sin3 θ
3(ΩR+u cos θ)3

2

4Ω2 u3 sin3 θ
(ΩR+u cos θ)5

2

+

2

u2 sin2 θ
(ΩR+u cos θ)3

M200 = 0.5 0.5ρair c 2lc Ω2

c1 −

β−

u sin θ
ΩR+u cos θ

+ 0.5ρair clc (ΩR + u cos θ)2 + u2 sin2 θ

M001 = 0.5ρair clc (ΩR + u cos θ)2 + u2 sin2 θ
β−

u3 sin3 θ
3(ΩR+u cos θ)3

Ωu sin θ
(ΩR+u cos θ)2

u3 sin3 θ
3(ΩR+u cos θ)3

c
c
3 c1 − 22 + c3 − 64

+(c2 + c4 )

−

−

Ωu sin θ
(ΩR+u cos θ)2

2

+ (c2 + c4 )

59

−

Ωu3 sin3 θ
(ΩR+u cos θ)4

Ωu sin θ
(ΩR+u cos θ)2

−

Ωu sin θ
(ΩR+u cos θ)2

Ωu3 sin3 θ −
(ΩR+u cos θ)4

M020 = 0.5 0.5ρair c 2lc

c1 −

c2
2

+ c3 −

c4
6

β−

u3 sin3 θ
3(ΩR+u cos θ)3

+

u sin θ
ΩR+u cos θ

3

+(c2 + c4 )
u3 sin3 θ
u sin θ
+ ΩR+u
cos θ
3(ΩR+u cos θ)3
1
u2 sin2 θ
ΩR+u cos θ − (ΩR+u cos θ)3

c
c
+ 4lc u sin θ 3 c1 − 22 + c3 − 64

β−

β−

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

2

+ (c2 + c4 )

1
ΩR+u cos θ

−

u2 sin2 θ
(ΩR+u cos θ)3

+lc (ΩR + u cos θ)2
+u2 sin2 θ
β−

c
c
c1 − 22 + c3 − 64

u3 sin3 θ
3(ΩR+u cos θ)3

1
ΩR+u cos θ

6

−

2

u2 sin2 θ
(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ





2
u3 sin3 θ
6u sin θ β−
+ u sin θ
 2u(c +c ) sin θ 
3
ΩR+u
cos
θ
3(ΩR+u cos θ)
2 4
−

−
 (ΩR+u cos θ)3 
(ΩR+u cos θ)3

c
c
M002 = 0.5 3.ρair clc c1 − 22 + c3 − 64

β−

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

(ΩR + u cos θ)2 + u2 sin2 θ

+ 0.

c
c
M110 = ρair clc u sin θ 3 c1 − 22 + c3 − 64

β−

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

2

+ (c2 + c4 )

+ 0.5ρair clc (ΩR + u cos θ)2 + u2 sin2 θ
3Ωu2 sin2 θ
(ΩR+u cos θ)4

−

Ωu3 sin3 θ
(ΩR+u cos θ)4

−

−

Ωu sin θ
(ΩR+u cos θ)2

Ωu sin θ
(ΩR+u cos θ)2

c
c
3 c1 − 22 + c3 − 64

Ω
(ΩR+u cos θ)2

u3 sin3 θ
u sin θ
+ ΩR+u
cos θ
3(ΩR+u cos θ)3
Ωu3 sin3 θ −
Ωu sin θ
(ΩR+u cos θ)4
(ΩR+u cos θ)2

β−

Ωu3 sin3 θ
(ΩR+u cos θ)4

2

c
c
+ 6 c1 − 22 + c3 − 64
u3 sin3 θ
3(ΩR+u cos θ)3

1
ΩR+u cos θ

−

u2 sin2 θ
(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

+ (c2 + c4 )

c
c
+ ρair clc Ω(ΩR + u cos θ) 3 c1 − 22 + c3 − 64

1
ΩR+u cos θ

β−

3Ωu2 sin2 θ
(ΩR+u cos θ)4

−

Ω
(ΩR+u cos θ)2

60

−

u2 sin2 θ
(ΩR+u cos θ)3
u3 sin3 θ
3(ΩR+u cos θ)3

β−

u sin θ
+ ΩR+u
cos θ

2

+ (c2 + c4 )

1
ΩR+u cos θ

u2 sin2 θ
(ΩR+u cos θ)3

−

c
c
M101 = ρair clc Ω(ΩR + u cos θ) 3 c1 − 22 + c3 − 64
u3 sin3 θ
3(ΩR+u cos θ)3

β−

+

u sin θ
ΩR+u cos θ

+c2 + c4 )
c
c
+ 3.ρair clc c1 − 22 + c3 − 64

2

(ΩR + u cos θ)2 + u2 sin2 θ

Ωu sin θ
Ωu3 sin3 θ −
(ΩR+u cos θ)4
(ΩR+u cos θ)2
u sin θ
u3 sin3 θ
+ ΩR+u
β−
cos θ
3(ΩR+u cos θ)3

M011 = ρair clc u sin θ 3 c1 −

c2
2

+ c3 −

c4
6

β−

u3 sin3 θ
3(ΩR+u cos θ)3

+

u sin θ
ΩR+u cos θ

2

+c2 + c4 )
c
c
+ 3.ρair clc c1 − 22 + c3 − 64

β−

u3 sin3 θ
3(ΩR+u cos θ)3

M300 = 16
β−

β−

c
c
0.5ρair c 6lc Ω2 3 c1 − 22 + c3 − 64

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

u3 sin3 θ
3(ΩR+u cos θ)3
u3 sin3 θ
3(ΩR+u cos θ)3

u3 sin3 θ
3(ΩR+u cos θ)3

2

+ (c2 + c4 )

c
c
c1 − 22 + c3 − 64
u sin θ
+ ΩR+u
cos θ

2

u sin θ
+ ΩR+u
cos θ

+lc (ΩR + u cos θ)2 + u2 sin2 θ
β−

−

u2 sin2 θ
(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

+6lc Ω(ΩR + u cos θ)
β−

1
ΩR+u cos θ

(ΩR + u cos θ)2 + u2 sin2 θ

+

u sin θ
ΩR+u cos θ

3

Ωu3 sin3 θ
(ΩR+u cos θ)4

Ωu3 sin3 θ
(ΩR+u cos θ)4
2Ω2 u sin θ
(ΩR+u cos θ)3

−

−

Ωu sin θ
(ΩR+u cos θ)2

Ωu sin θ
(ΩR+u cos θ)2

2 3 3
− 4Ω u sin θ5

(ΩR+u cos θ)
2
Ωu3 sin3 θ −
Ωu sin θ
+6
(ΩR+u cos θ)4
(ΩR+u cos θ)2
2Ω2 u sin θ − 4Ω2 u3 sin3 θ
+ (c2 + c4 )
(ΩR+u cos θ)3
(ΩR+u cos θ)5

c
c
c1 − 22 + c3 − 64
2

61

3

20Ω3 u3 sin3 θ
(ΩR+u cos θ)6

−

6Ω3 u sin θ
(ΩR+u cos θ)4

+18
β−

2Ω2 u sin θ − 4Ω2 u3 sin3 θ
(ΩR+u cos θ)3
(ΩR+u cos θ)5
u3 sin3 θ
u sin θ
+ ΩR+u
cos θ
3(ΩR+u cos θ)3

Ωu3 sin3 θ
(ΩR+u cos θ)4

+6

M030 = 61

−

+6

+ (c2 + c4 )


2 2
1
− u sin θ
ΩR+u cos θ (ΩR+u cos θ)3

1
ΩR+u cos θ

20Ω3 u3 sin3 θ
(ΩR+u cos θ)6

2

u sin θ
+ ΩR+u
cos θ


c4 
c2

+u2 sin2 θ 
 c1 − 2 + c3 − 6 −

−

+ (c2 + c4 )

3

Ωu sin θ
(ΩR+u cos θ)2

+lc (ΩR + 
u cos θ)2

36u sin θ

−

c
c
0.5ρair c 6lc 3 c1 − 22 + c3 − 64

u3 sin3 θ
3(ΩR+u cos θ)3

β−

Ωu3 sin3 θ
(ΩR+u cos θ)4

6 β−

β−

Ωu sin θ
(ΩR+u cos θ)2

1
ΩR+u cos θ

1
ΩR+u cos θ

−

−

−

6Ω3 u sin θ
(ΩR+u cos θ)4

u2 sin2 θ
(ΩR+u cos θ)3

u2 sin2 θ
(ΩR+u cos θ)3

2
u3 sin3 θ
u sin θ
+
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)3

u3 sin3 θ
+ u sin θ
3(ΩR+u cos θ)3 ΩR+u cos θ

(ΩR+u cos θ)3
3

−

2(c2 +c4 )
−
(ΩR+u cos θ)3

u2 sin2 θ
(ΩR+u cos θ)3

+ 6lc u sin θ

c1 −

c2
2

+ c3 −

c4
6

6

1
ΩR+u cos θ

−

u2 sin2 θ
(ΩR+u cos θ)3



β−

−

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

−

2
u sin θ
u3 sin3 θ
+
6u sin θ β−

3(ΩR+u cos θ)3 ΩR+u cos θ


(ΩR+u cos θ)3

2u(c2 +c4 ) sin θ
(ΩR+u cos θ)3

c
c
M003 = 61 3.ρair clc c1 − 22 + c3 − 64

(ΩR + u cos θ)2 + u2 sin2 θ + 0.

c
c
M210 = 0.25ρair c 2lc 3 c1 − 22 + c3 − 64

−

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ + β

2

+ (c2 + c4 )

62

1
ΩR+u cos θ
1
ΩR+u cos θ

−
−

u2 sin2 θ
(ΩR+u cos θ)3

u2 sin2 θ
(ΩR+u cos θ)3

Ω2

2

c
c
+4lc (ΩR + u cos θ) 3 c1 − 22 + c3 − 64

−

u3 sin3 θ
3(ΩR+u cos θ)3

1
ΩR+u cos θ

−

Ωu3 sin3 θ
(ΩR+u cos θ)4
3Ωu2 sin2 θ
(ΩR+u cos θ)4

(c2 + c4 )

−

Ωu sin θ
(ΩR+u cos θ)2

−

u2 sin2 θ
(ΩR+u cos θ)3

−

Ωu sin θ
(ΩR+u cos θ)2

3(ΩR+u cos θ)3

+3
+6

1
ΩR+u cos θ

u3 sin3 θ
3(ΩR+u cos θ)3

2

+ 12

u sin θ
+ ΩR+u
cos θ + β

+ (c2 + c4 )

(ΩR + u cos θ)2 + u2 sin2 θ
c
c
+ c1 − 22 + c3 − 64

3Ωu2 sin2 θ −
Ω
(ΩR+u cos θ)4
(ΩR+u cos θ)2
Ωu3 sin3 θ −
Ωu sin θ
(ΩR+u cos θ)4
(ΩR+u cos θ)2

2 2 2
3 3
u sin θ
− 12Ω u sin θ5
− u sin θ 3 + ΩR+u
cos θ
(ΩR+u cos θ)
3(ΩR+u cos θ)
u2 sin2 θ
2Ω2 u sin θ − 4Ω2 u3 sin3 θ
3
(ΩR+u cos θ)
(ΩR+u cos θ)3
(ΩR+u cos θ)5

−

u3 sin3 θ
3(ΩR+u cos θ)3

+3

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ + β

2

+β

u sin θ
+ ΩR+u
cos θ + β
2Ω2 u sin θ
(ΩR+u cos θ)3

+2lc u sin θ (c2 + c4 )
−

−

Ω
Ω + lc
(ΩR+u cos θ)2
2 2 2
2Ω2
− 12Ω u sin θ5
(ΩR+u cos θ)3
(ΩR+u cos θ)

2Ω2
(ΩR+u cos θ)3

−

Ω
(ΩR+u cos θ)2

c
c
+ 6 c1 − 22 + c3 − 64

−

Ωu3 sin3 θ
(ΩR+u cos θ)4
u3 sin3 θ

−

−

u2 sin2 θ
(ΩR+u cos θ)3

1
ΩR+u cos θ

6

2

u sin θ
+ ΩR+u
cos θ + β

3Ωu2 sin2 θ
(ΩR+u cos θ)4

2 3 3
− 4Ω u sin θ5
(ΩR+u cos θ)

Ωu3 sin3 θ
(ΩR+u cos θ)4

u sin θ
+ ΩR+u
cos θ + β

2Ω2 u sin θ
(ΩR+u cos θ)3

2 3 3
− 4Ω u sin θ5
(ΩR+u cos θ)

c
c
+ c1 − 22 + c3 − 64 (6

−

−

u3 sin3 θ
3(ΩR+u cos θ)3

c
c
M201 = 0.25ρair c 2lc Ω2 3 c1 − 22 + c3 − 64

Ωu sin θ
(ΩR+u cos θ)2

2
2

u sin θ
+ ΩR+u
cos θ + β

β−

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

+c2 + c4 )
c
c
+lc c1 − 22 + c3 − 64

β−

u3 sin3 θ
3(ΩR+u cos θ)3

(ΩR + u cos θ)2 + u2 sin2 θ

u sin θ
+ ΩR+u
cos θ

+6

Ωu3 sin3 θ
(ΩR+u cos θ)4

c
c
c1 − 22 + c3 − 64

63

6
−

2Ω2 u sin θ
(ΩR+u cos θ)3
Ωu sin θ
(ΩR+u cos θ)2

2 3 3
− 4Ω u sin θ5
(ΩR+u cos θ)

2

+ 24lc Ω

2

(ΩR + u cos θ)

Ωu3 sin3 θ
(ΩR+u cos θ)4

−

Ωu sin θ
(ΩR+u cos θ)2

β−

c
c
M021 = 0.25ρair c 2lc 3 c1 − 22 + c3 − 64

u3 sin3 θ
3(ΩR+u cos θ)3

β−

u sin θ
+ ΩR+u
cos θ

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

2

+c2 + c4 )
c
c
+24lc u sin θ c1 − 22 + c3 − 64

β−

u3 sin3 θ
3(ΩR+u cos θ)3

−

−

u2 sin2 θ
(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ

c
c
+lc c1 − 22 + c3 − 64
12u sin θ β−

1
ΩR+u cos θ

(ΩR + u cos θ)2 + u2 sin2 θ


u3 sin3 θ
+ u sin θ
3(ΩR+u cos θ)3 ΩR+u cos θ
(ΩR+u cos θ)3

6

1
ΩR+u cos θ

−

u2 sin2 θ
(ΩR+u cos θ)3




c
c
M120 = 0.25ρair c 4lc u sin θ 3 c1 − 22 + c3 − 64
3Ωu2 sin2 θ
(ΩR+u cos θ)4

−

u3 sin3 θ
3(ΩR+u cos θ)3

1
ΩR+u cos θ

−

Ωu3 sin3 θ
(ΩR+u cos θ)4
3Ωu2 sin2 θ
(ΩR+u cos θ)4
Ωu3 sin3 θ
(ΩR+u cos θ)4

−

−

Ω
(ΩR+u cos θ)2
u sin θ
+ ΩR+u
cos θ + β

u3 sin3 θ
3(ΩR+u cos θ)3

−

Ωu sin θ
(ΩR+u cos θ)2

−

−

Ω
(ΩR+u cos θ)2

c
c
+ 2lc 3 c1 − 22 + c3 − 64

−

Ωu sin θ
(ΩR+u cos θ)2
u sin θ
+ ΩR+u
cos θ + β

+lc (ΩR + u cos θ)2 + u2 sin2 θ

+12

c
c
+ 6 c1 − 22 + c3 − 64

u2 sin2 θ
(ΩR+u cos θ)3

u3 sin3 θ
3(ΩR+u cos θ)3

6

2

Ωu3 sin3 θ
(ΩR+u cos θ)4

−

3Ωu2 sin2 θ
(ΩR+u cos θ)4

−

2

Ωu3 sin3 θ
(ΩR+u cos θ)4

+ (c2 + c4 )

6(c2 +c4 )Ωu sin θ
(ΩR+u cos θ)4

Ωu sin θ
(ΩR+u cos θ)2

u sin θ
+ ΩR+u
cos θ + β

1
ΩR+u cos θ

Ω
(ΩR+u cos θ)2

64

−

+ (c2 + c4 )

Ωu sin θ
(ΩR+u cos θ)2

c
c
+ c1 − 22 + c3 − 64

−

u2 sin2 θ
(ΩR+u cos θ)3

2

2

1
u2 sin2 θ
ΩR+u cos θ − (ΩR+u cos θ)3
2
u3 sin3 θ
18Ωu sin θ −
+ u sin θ +β
3
ΩR+u
cos
θ
3(ΩR+u cos θ)
+
(ΩR+u cos θ)4
3
Ωu sin θ
u3 sin3 θ
Ωu sin3 θ −
−
+ u sin θ +β
12u sin θ
4
2
(ΩR+u cos θ)
(ΩR+u cos θ)
3(ΩR+u cos θ)3 ΩR+u cos θ
−
(ΩR+u cos θ)3

−

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ + β

+2lc Ω(ΩR + u cos θ)

c
c
c1 − 22 + c3 − 64

6

1
ΩR+u cos θ

−





u2 sin2 θ
(ΩR+u cos θ)3

2



−

−

u3 sin3 θ
3(ΩR+u cos θ)3

u sin θ
+ ΩR+u
cos θ + β

−

2
u3 sin3 θ
u sin θ +β
+
6u sin θ −

3(ΩR+u cos θ)3 ΩR+u cos θ

3

(ΩR+u cos θ)

2(c2 +c4 )u sin θ
(ΩR+u cos θ)3

c
c
M102 = 21 6.ρair clc Ω c1 − 22 + c3 − 64 (ΩR + u cos θ)
u3 sin3 θ
u sin θ
+ ΩR+u
cos θ
3(ΩR+u cos θ)3
c
c
+3.ρair clc c1 − 22 + c3 − 64 (ΩR + u cos θ)2
Ωu sin θ
Ωu3 sin3 θ −
(ΩR+u cos θ)4
(ΩR+u cos θ)2

β−

M012 = 12

u3 sin3 θ
u sin θ
+ ΩR+u
cos θ
3(ΩR+u cos θ)3
(ΩR + u cos θ)2 + u2 sin2 θ ΩR+u1 cos θ

c
c
6.ρair clc u sin θ c1 − 22 + c3 − 64

c
c
+3.ρair clc c1 − 22 + c3 − 64

−

+ u2 sin2 θ

β−

u2 sin2 θ
(ΩR+u cos θ)3

c
c
M111 = 6.ρair clc u sin θ c1 − 22 + c3 − 64
u3 sin3 θ
u sin θ
+ ΩR+u
cos θ
3(ΩR+u cos θ)3
1
u2 sin2 θ
β
ΩR+u cos θ − (ΩR+u cos θ)3

β−

Ωu3 sin3 θ
(ΩR+u cos θ)4

−

Ωu sin θ
(ΩR+u cos θ)2

c
c
+ 6.ρair clc Ω c1 − 22 + c3 − 64 (ΩR + u cos θ)

u3 sin3 θ
u sin θ
+ ΩR+u
cos θ
3(ΩR+u cos θ)3
c
c
+ 0.5ρair clc (ΩR + u cos θ)2 + u2 sin2 θ 6 c1 − 22 + c3 − 64

−

65

Ω
3Ωu2 sin2 θ −
(ΩR+u cos θ)4
(ΩR+u cos θ)2
u3 sin3 θ
u sin θ
β−
+ ΩR+u
cos θ
3(ΩR+u cos θ)3
Ωu3 sin3 θ −
Ωu sin θ
(ΩR+u cos θ)4
(ΩR+u cos θ)2

c
c
+ 6 c1 − 22 + c3 − 64

1
ΩR+u cos θ

−

u2 sin2 θ
(ΩR+u cos θ)3

Mathematica code for computing fijk ,
num[yd, z] = (yd - z Om + u Sin[th]);
den[y, zd] = (R + y) Om + zd + u Cos[th];
vsquared[y, yd, z, zd] = num[yd, z]^2 + den[y, zd]^2;
ratio[y, yd, z, zd] = num[yd, z]/den[y, zd];
mf[y, yd] = (ratio[y, yd, z, zd] /. {z -> 0} /. {zd -> 0});
Expand[Normal[Series[mf[y, yd], {y, 0, 3}, {yd, 0, 3}]]]
mfunc[y, yd] = (Om^2 y^2 yd)/(Om R + u Cos[th])^3 - (
Om y yd)/(Om R + u Cos[th])^2 + yd/(Om R + u Cos[th]) - (
Om^3 u y^3 Sin[th])/(Om R + u Cos[th])^4 + (
Om^2 u y^2 Sin[th])/(Om R + u Cos[th])^3 - (
Om u y Sin[th])/(Om R + u Cos[th])^2 + (u Sin[th])/(
Om R + u Cos[th]);
alp[y, yd, psi] =
bet + psi + Expand[((mfunc[y, yd]) - (1/3) (mfunc[y, yd])^3) + bet + psi];
Cl = c1*alp[y, yd, psi]*alp[y, yd, psi]*alp[y, yd, psi] +
c2*alp[y, yd, psi];
Cd = c3*alp[y, yd, psi]*alp[y, yd, psi] + c4;
z = 0;
zd = 0;
num[yd, 0] = (yd - 0 Om + u Sin[th]);
den[y, 0] = (R + y) Om + 0 + u Cos[th];
vsquared[y, yd, 0, 0] = num[yd, 0]^2 + den[y, 0]^2;
liftr[y, yd,
den[y, zd]
dragr[y, yd,
num[yd, z]

psi] =
Sqrt[vsquared[y, yd, 0, 0]] Cl*airdens*chord*0.5;
psi] =
Sqrt[vsquared[y, yd, 0, 0]] Cd*airdens*chord*0.5;

radial[y,yd,psi]= liftr[y,yd,psi]+dragr[y,yd,psi]

66

fhatexpansion[y,yd,psi] =
Normal[Series[radial[y,yd,psi], {y, 0, 3}, {yd, 0, 3}, {psi, 0, 3}]]
From fhatexpansion, the constant and linear coefficients are retained. In the example
below T ERM 0 = fˆ000 , T ERM 1 = fˆ100 , T ERM 2 = fˆ010 , T ERM 3 = fˆ001.
c1 = -2.9916;
c2 = 4.4287;
c3 = 1.185;
c4 = 0.0094;
chord = 0.533;
airdens = 1.292;
lc = 0.1302360914550899;
th = Om*t;
bet = 0;
TERM0 = (airdens*chord*0.5
(Om R + u Cos[th]) Sqrt[(Om R + u Cos[th])^2 +
u^2 Sin[th]^2] (c2 (bet + (u Sin[th])/(Om R + u Cos[th]) - (
u^3 Sin[th]^3)/(3 (Om R + u Cos[th])^3)) +
c1 (bet + (u Sin[th])/(Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
3 (Om R + u Cos[th])^3))^3) +
u Sin[th] Sqrt[(Om R + u Cos[th])^2 +
u^2 Sin[th]^2] (c4 +
c3 (bet + (u Sin[th])/(Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
3 (Om R + u Cos[th])^3))^2))
TERM1 = (airdens*chord*0.5
(Om R + u Cos[th]) Sqrt[(Om R + u Cos[th])^2 +
u^2 Sin[th]^2] (c2 (-((Om u Sin[th])/(Om R + u Cos[th])^2) + (
Om u^3 Sin[th]^3)/(Om R + u Cos[th])^4) +
3 c1 (-((Om u Sin[th])/(Om R + u Cos[th])^2) + (
Om u^3 Sin[th]^3)/(Om R + u Cos[th])^4) (bet + (u Sin[th])/(
Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
3 (Om R + u Cos[th])^3))^2) + (
Om (Om R +
u Cos[th])^2 (c2 (bet + (u Sin[th])/(Om R + u Cos[th]) - (
u^3 Sin[th]^3)/(3 (Om R + u Cos[th])^3)) +
c1 (bet + (u Sin[th])/(Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
3 (Om R + u Cos[th])^3))^3))/
Sqrt[(Om R + u Cos[th])^2 + u^2 Sin[th]^2] +
Om Sqrt[(Om R + u Cos[th])^2 +
u^2 Sin[th]^2] (c2 (bet + (u Sin[th])/(Om R + u Cos[th]) - (
u^3 Sin[th]^3)/(3 (Om R + u Cos[th])^3)) +
c1 (bet + (u Sin[th])/(Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
67

3 (Om R + u Cos[th])^3))^3) +
(2 c3 u Sin[th] Sqrt[(Om R + u Cos[th])^2 +
u^2 Sin[th]^2] (-((Om u Sin[th])/(Om R + u Cos[th])^2) + (
Om u^3 Sin[th]^3)/(Om R + u Cos[th])^4) (bet + (u Sin[th])/(
Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
3 (Om R + u Cos[th])^3)) + (
Om u (Om R + u Cos[th]) Sin[
th] (c4 +
c3 (bet + (u Sin[th])/(Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
3 (Om R + u Cos[th])^3))^2))/
Sqrt[(Om R + u Cos[th])^2 + u^2 Sin[th]^2]))
TERM2 = (airdens*chord*0.5
(Om R + u Cos[th]) Sqrt[(Om R + u Cos[th])^2 +
u^2 Sin[th]^2] (c2 (1/(Om R + u Cos[th]) - (
u^2 Sin[th]^2)/(Om R + u Cos[th])^3) +
3 c1 (1/(Om R + u Cos[th]) - (
u^2 Sin[th]^2)/(Om R + u Cos[th])^3) (bet + (u Sin[th])/(
Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
3 (Om R + u Cos[th])^3))^2) + (
u (Om R + u Cos[th]) Sin[
th] (c2 (bet + (u Sin[th])/(Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
3 (Om R + u Cos[th])^3)) +
c1 (bet + (u Sin[th])/(Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
3 (Om R + u Cos[th])^3))^3))/
Sqrt[(Om R + u Cos[th])^2 + u^2 Sin[th]^2] +
(2 c3 u Sin[th] Sqrt[(Om R + u Cos[th])^2 +
u^2 Sin[th]^2] (1/(Om R + u Cos[th]) - (
u^2 Sin[th]^2)/(Om R + u Cos[th])^3) (bet + (u Sin[th])/(
Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
3 (Om R + u Cos[th])^3)) + (
u^2 Sin[th]^2 (c4 +
c3 (bet + (u Sin[th])/(Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
3 (Om R + u Cos[th])^3))^2))/
Sqrt[(Om R + u Cos[th])^2 + u^2 Sin[th]^2] +
Sqrt[(Om R + u Cos[th])^2 +
u^2 Sin[th]^2] (c4 +
c3 (bet + (u Sin[th])/(Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
3 (Om R + u Cos[th])^3))^2)))
TERM3 = (airdens*chord*0.5
(Om R + u Cos[th]) Sqrt[(Om R + u Cos[th])^2 +
u^2 Sin[th]^2] (c2 +
3 c1 (bet + (u Sin[th])/(Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
3 (Om R + u Cos[th])^3))^2) +
68

(2 c3 u Sin[th] Sqrt[(Om R + u Cos[th])^2 +
u^2 Sin[th]^2] (bet + (u Sin[th])/(Om R + u Cos[th]) - (
u^3 Sin[th]^3)/(3 (Om R + u Cos[th])^3))))

Mathematica code for computing Mijk ,
num[yd, z] = (yd - z Om + u Sin[th]);
den[y, zd] = (R + y) Om + zd + u Cos[th];
vsquared[y, yd, z, zd] = num[yd, z]^2 + den[y, zd]^2;
ratio[y, yd, z, zd] = num[yd, z]/den[y, zd];
mf[y, yd] = (ratio[y, yd, z, zd] /. {z -> 0} /. {zd -> 0});
Expand[Normal[Series[mf[y, yd], {y, 0, 3}, {yd, 0, 3}]]]
mfunc[y, yd] = (Om^2 y^2 yd)/(Om R + u Cos[th])^3 - (
Om y yd)/(Om R + u Cos[th])^2 + yd/(Om R + u Cos[th]) - (
Om^3 u y^3 Sin[th])/(Om R + u Cos[th])^4 + (
Om^2 u y^2 Sin[th])/(Om R + u Cos[th])^3 - (
Om u y Sin[th])/(Om R + u Cos[th])^2 + (u Sin[th])/(
Om R + u Cos[th]);
alp[y, yd, psi] = ((mfunc[y, yd]) - (1/3) (mfunc[y, yd])^3) + bet +
psi;
Momentcof[y, yd, psi] =
alp[y, yd, psi]*(c4 + c2) +
alp[y, yd, psi]*alp[y, yd, psi]*
alp[y, yd, psi]*(c3 + c1 - (c2/2) - (c4/6));
z = 0;
zd = 0;
num[yd, 0] = (yd - 0 Om + u Sin[th]);
den[y, 0] = (R + y) Om + 0 + u Cos[th];
vsquared[y, yd, 0, 0] = num[yd, 0]^2 + den[y, 0]^2;
moment[y, yd, psi] =
Momentcof[y, yd, psi] vsquared[y, yd, 0, 0]*chord*lc*airdens*0.5;
Mhatexpansion[y,yd,psi] =
Normal[Series[moment[y,yd,psi], {y, 0, 3}, {yd, 0, 3}, {psi, 0, 3}]]

From Mhatexpansion, the constant and linear coefficients are retained. In the example
ˆ 000 , T ERM 1 = M
ˆ 100 , T ERM 2 = M
ˆ 010 , T ERM 3 = M
ˆ 001.
below T ERM 0 = M
c1 = -2.9916;
c2 = 4.4287;
c3 = 1.185;
69

c4 = 0.0094;
chord = 0.533;
airdens = 1292;
lc = 0.1302360914550899;
u = 4.72;
th = Om*t;
bet = 0;
R = 6;
TERM0 = 0.5‘ airdens chord lc ((Om R + u Cos[th])^2 +
u^2 Sin[th]^2) ((c2 + c4) (bet + (u Sin[th])/(
Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
3 (Om R + u Cos[th])^3)) + (c1 - c2/2 + c3 - c4/6) (bet + (
u Sin[th])/(Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
3 (Om R + u Cos[th])^3))^3)
TERM1 = 0.5‘ airdens chord lc ((Om R + u Cos[th])^2 +
u^2 Sin[th]^2) ((c2 +
c4) (-((Om u Sin[th])/(Om R + u Cos[th])^2) + (
Om u^3 Sin[th]^3)/(Om R + u Cos[th])^4) +
3 (c1 - c2/2 + c3 - c4/
6) (-((Om u Sin[th])/(Om R + u Cos[th])^2) + (
Om u^3 Sin[th]^3)/(Om R + u Cos[th])^4) (bet + (u Sin[th])/(
Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
3 (Om R + u Cos[th])^3))^2) +
1.‘ airdens chord lc Om (Om R +
u Cos[th]) ((c2 + c4) (bet + (u Sin[th])/(Om R + u Cos[th]) - (
u^3 Sin[th]^3)/(3 (Om R + u Cos[th])^3)) + (c1 - c2/2 + c3 c4/6) (bet + (u Sin[th])/(Om R + u Cos[th]) - (
u^3 Sin[th]^3)/(3 (Om R + u Cos[th])^3))^3)
TERM2 = 0.5‘ airdens chord lc ((Om R + u Cos[th])^2 +
u^2 Sin[th]^2) ((c2 + c4) (1/(Om R + u Cos[th]) - (
u^2 Sin[th]^2)/(Om R + u Cos[th])^3) +
3 (c1 - c2/2 + c3 - c4/6) (1/(Om R + u Cos[th]) - (
u^2 Sin[th]^2)/(Om R + u Cos[th])^3) (bet + (u Sin[th])/(
Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
3 (Om R + u Cos[th])^3))^2) +
1.‘ airdens chord lc u Sin[
th] ((c2 + c4) (bet + (u Sin[th])/(Om R + u Cos[th]) - (
u^3 Sin[th]^3)/(3 (Om R + u Cos[th])^3)) + (c1 - c2/2 + c3 c4/6) (bet + (u Sin[th])/(Om R + u Cos[th]) - (
u^3 Sin[th]^3)/(3 (Om R + u Cos[th])^3))^3)
TERM3 = 0.5‘ airdens chord lc ((Om R + u Cos[th])^2 +
70

u^2 Sin[th]^2) (c2 + c4 +
3 (c1 - c2/2 + c3 - c4/6) (bet + (u Sin[th])/(
Om R + u Cos[th]) - (u^3 Sin[th]^3)/(
3 (Om R + u Cos[th])^3))^2)
ˆ ijk are to be used in equation (2.37) ultimately to express
These values of fˆijk and M
equation (2.44) and finally (2.45) and (2.46).

71

BIBLIOGRAPHY

72

BIBLIOGRAPHY
[1] Annual Energy Outlook 2014 with projections to 2040, U.S. Energy Information Administration, April 2014, DOE-EIA0383(2014)
[2] Picture
taken
from
Johnson
System
Inc.
webpage,
http://www.johnsonsysteminc.com/wp-content/uploads/2011/10/Turbine1.jpg
[3] Picture taken from Shanghai StarCreation Group Co., Limited webpage,
http://www.shstarcreation.com/0/www.shstarcreation.com/upload/content/2012827456508
82.jpg.
[4] Eriksson, S., Bernhoff, H., and Leijon, M., 2008. ”Evaluation of different turbine concepts
for wind power”. Renewable and Sustainable Energy Reviews, 12, pp. 1419-1434.
[5] Barone, M., and Paquette, J., 2012. ”Vertical-axis wind turbines revisited: A Sandia perspective”. In 2012 Sandia Wind Turbine Blade Workshop.
http://energy.sandia.gov/wp/wp-content/gallery/uploads/3-B-1-Barone-SAND20125503C.pdf.
[6] Sutherland, H. J., Berg, D. E., Ashwill, T. D., 2012. ”A Retrospective of VAWT Technology” SAND2012-03014.
[7] Gorban, N. A., Gorlov, A. M., Silantyev, V. M.”Limits of the Turbine Efficiency for Free
Fluid Flow”, Journal of Energy Resources Technology, Vol. 123, pages 311 317.
[8] Chiras, D., 2010. Wind Energy Basics: A Green Energy Guide. New Society Publishers,
Gabriola Island, BC, Canada.
[9] Dabiri, J. O., 2011. ”Potential order-of-magnitude enhancement of wind farm power
density via counter-rotating vertical-axis wind turbine arrays”. Journal of Renewable
and Sustainable Energy, 3(4). article number 043104.
[10] Hau, E. ”Wind Turbines Fundamentals, Technologies, Application, Economics”,
Springer, Munich, 2012.
[11] Whittlesey, R. W., Liska, S., and Dabiri, J. O., 2010. ”Fish schooling as a basis for
vertical axis wind turbine farm design”. Bioinspiration and Biomimetics, 5(3). article
number 035005.
73

[12] Kinzel, M., Mulligan, Q., and Dabiri, J. O., 2012. ”Energy exchange in an array of
vertical-axis wind turbines”. Journal of Turbulence, 13(38), pp. 1-13.
[13] Ramakrishnan, V., and Feeny, B. F., 2011. ”In-plane nonlinear dynamics of wind turbine
blades”. In ASME International Design Engineering Technical Conferences, 23th Biennial
Conference on Vibration and Noise. paper number DETC2011-48219, on CD-ROM.
[14] Inoue, T., Ishida, Y., and Kiyohara, T., 2012. ”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, volume 134, issue 3, 13 pages.
[15] G. Acar, B.F. Feeny, 2014, ”A bend-bend-twist model for horizontal axis wind turbine
blades”,, 17th US National Congress on Theoretical and Applied Mechanics East Lansing.
[16] Kallesøe, B. S., 2007. ”Equations of motion for a rotor blade, including gravity and
pitch action”. Wind Energy, 10(3), February, pp. 209230.
[17] B.F. Feeny, 2013, ”Vertical-Axis wind turbine blade vibration modeling” Chemical, Bioengineering, Environmental and Transports Systems (CBET) proposal no. 1336817 and
Civil, Mechanical and Manufacturing Innovation (CMMI) proposal no. 1435126 (2014),
National Science Foundation.
[18] Worstell, M. H., ”Aerodynamic Performance of the 17-M-Diameter Darrieus Wind Turbine in the Three-Bladed Configuration: An Addendum”, SAND79-1753, August 1982.
[19] Caruntu, D., 2008. ”On nonlinear forced response of nonuniform blades”. In Proceeding of the 2008 ASME Dynamic Systems and Control Conference, no. DSCC2008-2157,
ASME.
[20] Meirovitch, L., 1997. Principles and Techniques in Vibrations. Prentice Hall, Upper
Saddle River.
[21] Jacobs, N. E, Ward, K. E., Pinkerton, R. M., ”The Characteristics of 78 Related Airfoil
Sections From Tests In The Variable-Density Wind Tunnel”, NACA Report 460.
[22] Worstell, M. H., ”Aerodynamic Performance of the 17 Meter Diameter Darrieus Wind
Turbine”, SAND78-1737, September 1978.
[23] Charles L. L., Cuyler W. B., Acquilla S. H.,Darrell W. S., 1996, ”Computer Program
To Obtain Ordinates for NACA Airfoils”, NASA Technical Memorandum 4741.

74

[24] Dowell, E. H., Curtiss Jr., H. C., Scanlan, R. H., Sisto F., A Modern Course in Aeroelasticity, 1989.
[25] Beer, F. P., Johnston, E. R., 1981. Mechanics of Materials. McGraw-Hill.

75