MODAL ANALYSIS OF VERTICAL-AXIS DARRIEUS WIND TURBINE BLADE WITH A
TROPOSKEIN SHAPE
By
Amr Fawzy Abdel hakeem Saleh

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

ABSTRACT
MODAL ANALYSIS OF VERTICAL-AXIS DARRIEUS WIND TURBINE BLADE WITH A
TOPOSKEIN SHAPE
By
Amr Fawzy Abdel hakeem Saleh
Darrieus wind turbines with troposkein shaped blades are an important type of vertical
axis wind turbines. When designing these turbines it is important to consider how
vibrations may affect blade failure. In order to avoid resonance, the blade natural
frequencies need to be determined. The goal of this research is to determine the blade free
vibration mode shapes and modal frequencies, neglecting the variation of the upstream
wind speed for simplicity. The blade modal vibration is studied numerically using ANSYS
software, and analytically using thin-beam theory and an assumed modes method. Firstly,
the analysis is performed on a 17 m diameter Sandia simplified troposkein shaped blade
with a NACA 0015 airfoil. The first ten modes and the corresponding natural frequencies
were calculated using ANSYS. The analysis was done for stationary and spinning blade at
54 rpm. The spin case compared very well to previous Sandia results. The same analysis
was then applied to an ideal troposkein shaped blade, and a slightly modified cross section.
Results show consistency across the blade shapes. The Sandia blade shape was studied
analytically using thin-beam theory and an assumed modes method. Kinetic and potential
energies were derived. Mass and stiffness matrices are formulated for the discretized
model, to which modal analysis is applied. The natural frequencies were calculated
analytically for same simplified Sandia blade shape used for finite element model under
stationary condition. The results of the two methods were compared for consistency.

Copyright by
AMR FAWZY ABDEL HAKEEM SALEH
2017

I would like to dedicate this work to my mother and my father’s soul, and my lovely wife

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.

v

TABLE OF CONTENTS
LIST OF TABLES………………………………………………………………………….................………….………. viii
LIST OF FIGURES…………………………………………………………………………………………………............. ix
1

Chapter 1: Introduction and literature review......................................................................................1
1.1 Background ...................................................................................................................................................1
1.2 Types of wind converters .......................................................................................................................2
1.2.1
Constructional design (rotor axis orientation) based classification ........................2
1.2.2
Classification based on ways of extracting the energy from wind (lift or drag
type)...................................................................................................................................................................... 4
1.3 Literature survey of the structural vibration analysis of vertical axis troposkein
shaped wind turbine blade ..................................................................................................................................6
1.4 Motivation ......................................................................................................................................................9
1.5 Objective ...................................................................................................................................................... 10
1.6 Thesis outline ............................................................................................................................................ 10

2 Chapter 2: Finite element model for calculating mode shapes and modal frequencies of
troposkein shaped vertical axis wind turbine blades .............................................................................. 11
2.1 Introduction ............................................................................................................................................... 11
2.2 3-D model for 17m diameter blade of straight-circular-straight shape ....................... 12
2.2.1
Blade cross section ........................................................................................................................ 12
2.2.2
Blade profile...................................................................................................................................... 13
2.3 FEM model for 17m diameter blade of straight-circular-straight shape spins with
54 rpm constant rotational velocity ............................................................................................................. 14
2.4 FEM model for 17m diameter stationary blade of straight-circular-straight shape
..................................................................................................................................................................... 23
2.5 FEM model for 17m diameter ideal troposkein shaped blade .......................................... 24
2.5.1
3-D model for ideal troposkein shape blade ..................................................................... 25
2.5.2
FEM modal analysis for 17m diameter blade with ideal troposkein shape ....... 27
3 Chapter 3: Analytical beam-based formulation of troposkein shaped blade modal
analysis ........................................................................................................................................................................... 35
3.1 Introduction ............................................................................................................................................... 35
3.2 Blade potential energy .......................................................................................................................... 35
3.2.1
Derivation of strain equations ................................................................................................. 39
3.3 Kinetic energy ........................................................................................................................................... 48
3.4 Discretization of energy expressions for modal analysis ..................................................... 53
3.4.1
Derivation of stiffness matrix expressions ........................................................................ 55
3.4.2
Derivation of mass matrix expressions ............................................................................... 58
3.4.3
Meridian angle (θ) ......................................................................................................................... 60
3.5 Numerical calculation of modal frequencies of 17 m diameter Sandia simplified
blade shape...……………………………………………………………………….……………………………........... 60
3.5.1
17 m diameter simplified Sandia blade cross section properties ........................... 61
vi

3.5.2
4

Numerical values for first ten modal frequencies .......................................................... 62

Chapter 4 .............................................................................................................................................................. 65
4.1 Conclusion and comments .................................................................................................................. 65
4.2 Contribution ............................................................................................................................................... 66
4.3 Future work................................................................................................................................................ 66

APPENDIX………………………………………………………………………………………………….……………..... 66
BIBLIOGRAPHY………………………………………………………………………………………………….………... 70

vii

LIST OF TABLES
Table 2.1 First 10 natural frequencies of 17 m diameter Sandia blade…………………………. 16
Table 2.2 Sandia (MSU) stationary blade first ten modal frequencies ……………………...…... 23
Table 2.3 cross sectional properties for both Sandia and Ideal troposkein blades under
discussion…………………………………………………………………………………………………………...…….. 25
Table.2.4 first 10 natural frequencies of 17 m diameter ideal troposkein shape blade…... 28
Table 2.5 Comparison of the lower modal frequencies as obtained by finite element
analysis. "Sandia" refers to the straight-circular-straight blade model. Spinning occurs at 54
rpm…………………………………………………………………………………………………………………..………. 33
Table 3.1 Modal frequencies for 17 m diameter Sandia stationary blade (coupling between
all modes considered……………………………………………………..……………………………….………..… 62

viii

LIST OF FIGURES
Figure 1.1 Global wind power cumulative capacity [1]..............................................................................1
Figure 1.2 Schematic of vertical axis wind turbine [2] ...............................................................................2
Figure 1.3 Schematic of horizontal wind turbine [2] ..................................................................................3
Figure 1.4 Schematic of Savonius wind turbine [3] .....................................................................................4
Figure 1.5 Schematic of Darrieus wind turbine [4] .....................................................................................5
Figure 1.6 Schematic of a perfectly flexible cable rotating about a vertical axis [7] ....................6
Figure 1.7 Darrieus vertical-axis wind turbine (DOE/Sandia 34-m) [6] ...........................................8
Figure 2.1 plot of airfoil points ............................................................................................................................ 13
Figure 2.2 NASA 0015 blade cross section of 17 m diameter Sandia blade. (a) The original
Sandia one. (b) The one rebuilt for this work .............................................................................................. 13
Figure 2.3 17 m diameter Sandia blade profile (all dimensions in m) ............................................. 14
Figure 2.4 Sandia simplified blade shape clamped at both ends and allowed to rotate ......... 15
Figure 2.5 17 m diameter Sandia blade meshing pattern ...................................................................... 16
Figure 2.6 first mode of Sandia blade at 54 rpm (flat-wise only), 2.0348 Hz ............................... 17
Figure 2.7 second mode of Sandia blade at 54 rpm (edge-wise only), 2.5675 Hz...................... 17
Figure 2.8 third mode of Sandia blade at 54 rpm (flat-wise only), 4 HZ......................................... 18
Figure 2.9 fourth mode of Sandia blade at 54 rpm (edge-wise and torsion), 5.936 Hz........... 18
Figure 2.10 fifth mode of Sandia blade at 54 rpm (flat-wise only), 6.1904 Hz ............................ 19
Figure 2.11 sixth mode of Sandia blade at 54 rpm (flat-wise only), 8.5349 Hz........................... 19
Figure 2.12 seventh mode of Sandia blade at 54 rpm (flat-wise only), 11.749 Hz .................... 20
Figure 2.13 eighth mode of Sandia blade at 54 rpm (edge-wise and torsion), 13.52 Hz ........ 20
Figure 2.14 ninth mode of Sandia blade at 54 rpm (flat-wise only), 15.335 Hz.......................... 21
ix

Figure 2.15 tenth mode of Sandia blade at 54 rpm (flat-wise only), 19.354 Hz.......................... 21
Figure 2.16 mode shapes and frequencies of a single blade of a spinning rotor by Sandia .. 22
Figure 2.17 NASA 0015 blade cross section for 17 m diameter ideal troposkein shape blade
............................................................................................................................................................................................. 24
Figure 2.18 Troposkein shape ............................................................................................................................. 26
Figure 2.19 17 m diameter troposkein shape with H= 12.238 m....................................................... 27
Figure 2.20 first mode for the troposkein blade at 54 rpm (flat-wise only), 1.87 Hz............... 28
Figure 2.21 second mode for the troposkein blade at 54 rpm (edge-wise only), 2.34 Hz ..... 29
Figure 2.22 third mode for the troposkein blade at 54 rpm (flat-wise only), 3.68 Hz............. 29
Figure 2.23 fourth mode for the troposkein blade at 54 rpm (flat-wise only), 5.82 Hz .......... 30
Figure 2.24 fifth mode for the troposkein blade at 54 rpm (edge-wise & torsion), 6.2 Hz.... 30
Figure 2.25 sixth mode for the troposkein blade at 54 rpm (flat-wise only), 8.4 Hz ................ 31
Figure 2.26 seventh mode for the troposkein blade at 54 rpm (flat-wise only), 11.36 Hz .... 31
Figure 2.27eighth mode for the troposkein blade at 54 rpm (edge-wise / torsion), 13.3 Hz
............................................................................................................................................................................................. 32
Figure 2.28 ninth mode for the troposkein blade at 54 rpm (flat-wise only), 14.82 Hz ......... 32
Figure 2.29 tenth mode for the troposkein blade at 54 rpm (flat-wise only), 18.74 Hz ......... 33
Figure 3.1 vertical-axis wind turbine and coordinate systems [8] .................................................... 37
Figure 3.2 coordinate systems of blade cross section [8] ...................................................................... 37
Figure 3.3 axis of blade before and after deformation, and coordinate systems [8] ................ 38

x

1 Chapter 1: Introduction and literature review
1.1 Background
The development of renewable energy represents an obvious need especially when we see
based on modern studies that fossil fuel sources like oil and gas reserve will be depleted in
near future. Renewable energy can be solar energy, wind energy, geothermal energy,etc.
Among the renewable energy alternatives, wind energy introduces itself as one of the most
prominent sources of renewable energy. Wind turbines are the most popular machines that
can convert wind power into mechanical power. Wind turbines convert the wind power
into electricity via rotation in the generator. Now, after large technological improvements,
we are able to develop wind turbines to be more feasible, reliable, and dependable source
of energy especially for electricity. Indeed, a significant share of electricity for several
countries comes from wind turbines. Figure 1.1 indicates the growth of global wind power
production over years from 1996 to 2015 [1].

Figure 1.1 Global wind power cumulative capacity [1]

1

1.2 Types of wind converters
Wind energy converters can be classified firstly in accordance with their constructional
design (rotor axis orientation) and, secondly, according to their ways of extracting the
energy from wind (lift or drag type) [2].
1.2.1 Constructional design (rotor axis orientation) based classification
1.2.1.1 Vertical axis wind turbines (VAWTs):
As shown in Figure.1.2. The axis of rotation of this type is perpendicular to the wind
direction [2]. In this type the generator is on the ground which makes it more accessible
and no yaw system required.

Figure 1.2 Schematic of vertical axis wind turbine [2]

The constant blade profile and cross sectional shape of an H-rotor type VAWT along the
blade length is one of its advantages because it makes every blade section be subjected to
the same wind speed. As a result blade twisting not required. Troposkein-shaped turbines
also tend to have constant cross-sectional shapes.

2

The constant blade shape makes its design and manufacturing easier and the total
production cost lower.
1.2.1.2 Horizontal axis wind turbines (HAWTs):
As shown in figure 1.3 [2] the axis of rotation of this type must be oriented parallel to the
wind in order to capture the wind power. HAWTs are considered to be the most efficient
turbines. The disadvantages of this type include:
a. Operation at high starting wind velocity.
b. Low starting torque.
c. Yaw mechanism required to turn the rotor toward the wind.
d. Power loss when the rotors are tracking the wind directions.
e. High center of gravity.

Figure 1.3 Schematic of horizontal wind turbine [2]

3

1.2.2 Classification based on ways of extracting the energy from wind (lift or drag
type)
The rotor’s aerodynamics is classified based on whether the wind energy converter
captures its power from the aerodynamic drag of the air flow over the rotor surfaces, or
whether it captures its power from the aerodynamic lift generated from the air flow against
the airfoil surfaces [2]. Based on that, there are so-called drag-type rotors and lift- type
rotors.
1.2.2.1.1 The drag type
This turbine takes less energy from the wind but provides a higher torque and is suitable
for mechanical applications as water pumping [2]. The most representative model of a
drag-type VAWT as shown in Figure 1.4 [3] is the Savonius.

Figure 1.4 Schematic of Savonius wind turbine [3]

1.2.2.2 The lift type
This turbine generates power mainly by the generated lift force on the blade cross section.
It can move quicker than the free wind speed. The most important application of this kind

4

is electricity generation [2]. The most representative model of a lift-type VAWT as shown in
Figure 1.4 [4] is the Darrieus wind turbine.

Figure 1.5 Schematic of Darrieus wind turbine [4]

In 1931 the U.S. Patent Office patented Darrieus wind turbine in the name of G.J.M. Darrieus
[5]. The Darrieus patent states that “each blade should have a streamline outline curved in
the form of skipping rope.” [5] This means, the shape of the Darrieus rotor blades can be
approximated to the shape of a perfectly flexible cable, of uniform density and cross
section, anchored from two fixed points and rotating about its long axis; under the effect of
centrifugal forces such a shape minimizes inherent bending stresses. This blade shape is
called troposkien (from the Greek roots: τρoτs, turning and σχolυloυ, rope). A pure
troposkien shape (gravity neglected) does not depend on angular velocity [6]. The equation
that defines a troposkein blade profile was developed by Sandia in April 1974 [7]. The
analysis considered a perfectly flexible cable rotating about a fixed axis at a constant
angular velocity without gravity as shown in figure 1.6. The Darrieus wind turbine with a
troposkien shaped blades has advantage of working under the effect of high centrifugal forces
without failure.
5

Figure 1.6 Schematic of a perfectly flexible cable rotating about a vertical axis [7]

1.3 Literature survey of the structural vibration analysis of vertical axis
troposkein shaped wind turbine blade
The performance of a Darrieus wind-turbine was studied firstly by R.S. Rangi and P. South
using wind tunnel measurements in the National Research Council of Canada in March
1971 [8]. In February 1974, R.S. Rangi and P. South and their team investigated Darrieus
wind parameters including spoilers and aero-brakes effect on turbine performance and
reliability, in addition to the effect of efficiency related parameters like the number of
blades and the rotor’s solidity. Compared to Sandia’s simplified troposkein shape, the
engineers at the National Research Council of Canada (NRC) in the early 1970’s,
independently developed catenary shape as an approximation to the troposkein curved
6

blade [9]. In December 1979 a NASA team developed second degree nonlinear aero-elastic
partial differential equations of motion for a slender, flexible, non-uniform, Darrieus
vertical axis wind turbine blade using Hamilton’s principle. The analysis of NASA team
considered a blade undergoing combined flat-wise bending, edgewise bending, torsion, and
extension. Flatwise and edgewise are defined as follows. The airfoil cross section has a
major (chord-length) and minor axis. Bending about the major axis is referred to as “flat
wise”, and bending about the minor axis is called “edge-wise.” Extension refers to axial
deformation, and torsion occurs about the elastic axis [can cite an elasticity book regarding
elastic axis].The blade aero-dynamic loading was developed using strip theory based on a
quasi-steady approximation of two-dimensional incompressible unsteady airfoil theory [10].
The derivation of the equations of motion was done based on the geometric nonlinear
theory of elasticity in which the elongations and shears (and hence strains) are negligible
compared to unity. In this research work The NASA model was used to derive equations of
the blade strain energy and kinetic energy. These equations are suitable to study the blade
free vibration and the blade dynamic response [10].
In June 1979, Sandia introduced a modal analysis for a 17-m diameter simplified
troposkein shape rotating at 54 rpm constant angular velocity using ANSYS finite element
software. This simplified troposkein blade consists of two straight increments at both ends
and circular part in the middle (described in details later) with a NASA 0015 airfoil cross
section. The first eight natural frequencies were calculated and the corresponding mode
shapes were produced [11].
Another wind tunnel performance analysis for the Darrieus wind turbine with NACA 0012
Blades was performed in Sandia by Bennie F. Blackwell, Robert E. Sheldahl, and Louis V.
7

Feltz [12], where several configurations of Darrieus wind turbine blade with NASA 0012
airfoil cross section and 2- meters diameter were tested in a low-speed wind tunnel to
measure the output torque, rotational speed as a function of rotor solidity, Reynolds
number, and the free stream wind velocity. A significant step in the development of larger
and more efficient commercial Darrieus VAWT’s was the installation and operation of 34-m
Sandia-DOE VAWT in 1987, rated at 625 kW. The Sandia 34-m turbine (Fig. 1.7) was the
first curved-blade Darrieus turbine rotor originally designed to incorporate step-tapered
blades using varying blade-section airfoils and a blade airfoil section specifically designed
for VAWTs [6].

Figure 1.7 Darrieus vertical-axis wind turbine (DOE/Sandia 34-m) [6]

In 2014 a structural dynamic design tool developed by Sandia for large scale VAWTs for
studying the effect of geometry configuration, blade material, and number of blades on the
aeroelastic stability of VAWTs. This tool can describe quantitatively the aeroelastic
8

instabilities in VAWT design [13]. This is part of a recent strong research effort led by
Sandia to study VAWTs, motivated by the prospect of off-shore applications. Their work
has addressed the modeling of aeroelastic loading and the vibration responses, ultimately
for guiding the design and manufacture of VAWT systems, and has involved the
development of finite-element-based OWENS simulation toolkit, modal analysis of bladetower systems, potential for resonances, and field tests on the Sandia 34-m VAWT test bed
[14-16]. Further interest in structural modeling of VAWTs is represented by a series of
presentations at a recent conference dedicated to off-shore VAWTs [17].

1.4 Motivation
The following facts indicate why we need deep and active research on Darrieus vertical axis
wind turbines:
1. Limited knowledge and experience about Darrieus VAWTs
2. Large scaled HAWTs suffer failures due cyclic gravitational and aerodynamic loadings.
Darrieus VAWTs offers a good solution for this problem because of the blade shape
which minimizes inherent bending stresses and rotor vertical position which minimizes
the gravitational effect, so by developing a large scale Darrieus VAWTs and improving
VAWT dynamic behavior can provide an alternative to problematic HAWTs without
losing advantage of high energy production.
3. Further investigation and more study of VAWTs dynamics and vibration can increase its
life and as a result save a lot of money
4. All research in the field of troposkein shape blade dynamics was basically for the
simplified and approximated shape. Now with the help of high performance computers

9

and advanced calculation techniques, it is time to investigate the ideal troposkein blade
dynamics and how we can optimize this shape to get efficient blade performance.

1.5 Objective
This research aims to formulate a robust model for a troposkein shaped VAWT blade for
studying the blade vibration and its mode shapes and natural frequencies, taking into
account the blade loading and geometric complexities, which help us determine the safe
margins for the blade operating loads and angular velocities.
1.6 Thesis outline
This research is presented in the next two chapters:
Chapter 2 focuses on the finite-element analysis of the blade vibration, while chapter 3
sketches an analytical beam-based model for vibration studies.
More specifically, chapter 2 consists of the following:
1. Modal analysis is done using ANSYS workbench software for Darrieus wind turbine
blade with NASA 0015 airfoil for both simplified and ideal troposkein shapes.
2. Model validation is performed by comparing the mode shapes and natural frequencies
of the first eight modes to those of the Sandia model [11].
3. The effect of spinning on the model vibration is then examined
Chapter 3 consists of the following:
4. Kinetic and potential energies equations are derived and mass and stiffness matrices
are formulated to be used in Lagrange’s equation to calculate the natural frequencies
and mode shapes analytically.

10

2 Chapter 2: Finite element model for calculating mode shapes
and modal frequencies of troposkein shaped vertical axis
wind turbine blades
2.1 Introduction
The wind turbine structure will oscillate due to cyclic aerodynamic loads during operation.
Because the turbine is a continuous system with many degrees of freedom it has many
different modes of oscillation and many corresponding natural frequencies. The most
critical component of the wind turbine is the blades. In this chapter a finite element modal
analysis has been done using ANSYS WORKBENCH finite element software.
In the first part of this chapter a model for a 17-m diameter straight-circular-straight shape
suggested by Sandia as an approximation to the ideal troposkein shape has been built using
inventor software. This wind turbine blade model has a 17-m diameter with 25.15 m height
and 30.48 m total blade length. The blade cross section is a NASA-0015 airfoil with chord
length 0.61 m, airfoil thickness of 0.091 m and wall thickness 6.35 mm except for the
reinforced leading edge which has a wall thickness of 8.89 mm.
The blade cross section has four uniform internal stiffeners with the same thickness of the
wall. The stiffeners are strategically located to provide enough structural rigidity.
A free vibration modal analysis for this model has been performed using ANSYS
WORKBENCH under constant rotational velocity of 54 rpm with blade ends clamped. To
validate this model the results of the analysis compared to the results of Sandia report [11].

11

After validating our model, the vibrations of a stationary blade are studied. A stationary
blade condition is useful as a benchmark for stationary laboratory experiments or
comparison between models. It also helps to understand the effect of spin.
In the second part of this chapter a model for the ideal troposkein shape has been built and
analyzed to get the mode shapes and modal frequencies of the ideal troposkein shape and
to compare the approximated blade shape of the straight-circular-straight segmented beam
to the ideal shape.

2.2 3-D model for 17m diameter blade of straight-circular-straight
shape
2.2.1 Blade cross section
The blade airfoil is based on NASA’s four digit airfoils equation [18] as follows
𝑡

𝑥

𝑥

𝑥

𝑥

𝑥

𝑦𝑡 = 0.2 (0.2969√ 𝑐 − 0.126 𝑐 − 0.3516(𝑐 )2 + 0.2843(𝑐 )3 − 0.1015(𝑐 )4 )

(2.1)

where
c is the chord length
t is the maximum airfoil thickness which equals 0.15c
x is the coordinates of the airfoil along the chord or the major axis
𝑦𝑡 is the coordinates of the airfoil along the minor axis

For a NASA 0015 airfoil with 0.61 m cord length, the blade cross section shape can be
calculated using equation (2.1) and its coordinates as given in table A.1 in appendix A.
Plotting y versus x from table A.1 gives the airfoil shape as in figure 2.1

12

Figure 2.1 plot of airfoil points

Taking points of airfoil from table.1in the appendix and generating the shape of the outer
and inner airfoils by Inventor software using the offset feature, then a four stiffening spars
was generated by Inventor with 6.5 mm thickness each and strategically located to provide
proper structural rigidity [11]. The final shape of the airfoil is shown in figure 2.2

(a)

(b)

Figure 2.2 NASA 0015 blade cross section of 17 m diameter Sandia blade. (a) The
original Sandia one. (b) The one rebuilt for this work

2.2.2 Blade profile
The Sandia blade consists of three segments: one circular at the middle with radius 10.43 m
and 16.754 m length and two straight segments of 6.858 m length each as indicated in
figure 2.3.
13

Figure 2.3 17 m diameter Sandia blade profile (all dimensions in m)

2.3 FEM model for 17m diameter blade of straight-circular-straight
shape spins with 54 rpm constant rotational velocity
A modal analysis for the 17 m diameter Sandia blade has been done using ANSYS
workbench by meshing the blade structure using the SOLID 187 element. The SOLID 187 is
a ten nodes higher order 3-D element, which makes it well suited to model irregular
meshes. A SOLID 187 element has three degrees of freedom at each node (translations in
the nodal x, y, and z directions) [19]. The blade is clamped at both ends except for rotation
about the y-axis as shown in figure 2.4 where the blade assumed to be rotating at constant
rotational velocity of 54 rpm.

14

Figure 2.4 Sandia simplified blade shape clamped at both ends and allowed to rotate

Figure 2.5 indicates the mesh pattern generated by ANSYS. By magnifying a part of this
mesh we can note the elements shapes. For this analysis of 17 m diameter simplified Sandia
blade the element size refined several times starting from automatic mesh and going finer
to see the effect of mesh refinement on the results. After refining the mesh in limits of
ANSYS academic version available for our lab, we found that there is no significant
difference in the results of different meshing cases. The modal results of simplified Sandia
blade for both spinning and stationary cases calculated with mesh element size of 0.07m
and blade total number of elements equals 72562 elements.
The first 10 modes and the corresponding natural frequencies of the Sandia blade rotating
at 54 rpm constant rotational velocity have been calculated, and the results for the natural
frequencies are indicated in table 2.1, and the modes are shown in figures 2.6 to 2.15.

15

Figure 2.5 17 m diameter Sandia blade meshing pattern

Table 2.1.first 10 natural frequencies of 17 m diameter Sandia blade

16

Figure 2.6 first mode of Sandia blade at 54 rpm (flat-wise only), 2.0348 Hz

Figure 2.7 second mode of Sandia blade at 54 rpm (edge-wise only), 2.5675 Hz

17

Figure 2.8 third mode of Sandia blade at 54 rpm (flat-wise only), 4 HZ

Figure 2.9 fourth mode of Sandia blade at 54 rpm (edge-wise and torsion), 5.94 Hz

18

Figure 2.10 fifth mode of Sandia blade at 54 rpm (flat-wise only), 6.1904 Hz

Figure 2.11 sixth mode of Sandia blade at 54 rpm (flat-wise only), 8.5349 Hz

19

Figure 2.12 seventh mode of Sandia blade at 54 rpm (flat-wise only), 11.749 Hz

Figure 2.13 eighth mode of Sandia blade at 54 rpm (edge-wise and torsion), 13.52 Hz

20

Figure 2.14 ninth mode of Sandia blade at 54 rpm (flat-wise only), 15.335 Hz

Figure 2.15 tenth mode of Sandia blade at 54 rpm (flat-wise only), 19.354 Hz

21

Comparing the above results to Sandia results [11] shown in figure 2.16; we find that the
results are approximately the same. This validates the finite element analysis based on the
ANSYS WORKBENCH and SOLID 187 element, which will then be applied to stationary and
ideal troposkein blades.

Figure 2.16 mode shapes and frequencies of a single blade of a spinning rotor by Sandia

We can make some observations from these results. Figures 2.6 and 2.10 show
antisymmetric flat-wise bending modes. The zero-point on these mode shapes is not
exactly centered. Indeed, observing animations shows that the zero point moves during
pure modal vibrations. Since the stiffness and mass matrices are symmetric, this zero point
shifting suggests that these flat-wise modes also involve extension. This will be
22

investigated later with both finite-elements and beam-based modeling. The Sandia results
also suggest this behavior. Also, the dominantly flat-wise (or possibly flat-wise/extension)
and edge/torsion mode shapes suggests that the flat-wise and extension deformation
coordinates may be nearly decoupled from the edge-wise and torsion deformation
coordinates. A flatwise mode similar to a first mode of a clamped-clamped beam does not
show up among the lower modes. Such a mode of flatwise deformation will necessarily
involve significant extension and compression. Since extension is in a stiff orientation of a
blade (axial deformation tends to be stiffer than bending), such a mode is likely to have a
high frequency.

2.4 FEM model for 17m diameter stationary blade of straight-circularstraight shape
We have validated our model by comparing the calculated mode shapes and corresponding
modal frequencies of the Sandia shape spinning blade to those which were calculated by
Sandia. Now we will consider the mode shapes and modal frequencies of the same blade
analyzed in section 2.3 but under stationary conditions, which means the blade rotational
velocity is zero. For brevity the mode shapes will not be plotted. For the same element size
and meshing pattern, the stationary case modal analysis results can be summarized as
shown in table 2.2 below.
Table 2.2 Sandia (MSU) stationary blade first ten modal frequencies
Mode No.

Freqency shape

1

1.2824

Flat-wise

2

2.4595

Edge-wise

3

2.8691

Flat-wise
23

Table 2.2 (cont’d)
4

4.8445

Flat-wise

5

5.6198

Edge/torsion

6

7.1157

Flat-wise

7

10.188

Flat-wise

8

13.177

Edge/torsion

9

13.693

Flat-wise

10

17.649

Flat-wise

2.5 FEM model for 17m diameter ideal troposkein shaped blade
In this part a modal analysis has been done for an ideal troposkein shape blade with 17 m
diameter and blade cross section as in figure 2.17. Here the cross section is a little bit
different than the simplified Sandia shape blade to overcome the limitations of the number
of elements in the ANSYS software academic version available in our lab at Michigan State
University.

Figure 2.17 NASA 0015 blade cross section for 17 m diameter ideal troposkein shape blade

24

Comparing the ideal troposkein blade cross section shown in 2.17 to the simplified Sandia
blade cross section shown in figure 2.2 we can summarize the differences as in table 2.3
below.
Table 2.3 cross sectional properties for both Sandia and Ideal troposkein blades under
discussion
Property

Sandia cross section

Airfoil / chord length

NASA 0015/0.61 m

Ideal troposkein cross
section
NASA 0015/ 0.61m

Cross section area
Flat wise moment of inertia
Edge wise moment of inertia
Total mass of 17 m diameter
blade
Blade material

0.98 ∗ 10−2 m2 = 15.129 in2
0.859 ∗ 10−5 m^4
2.814 ∗ 10−4 m^4

0.96613 ∗ 10−2 m2 = 14.975 in2
0.819 ∗ 10−5 m^4
2.816 ∗ 10−4 m^4

810.82 Kg

802 Kg

6063-T6 Aluminum alloy

6063-T6 Aluminum alloy

From table 2.3 we can see some differences in the blade total mass and edge wise moment
of inertia and we see later how these differences will contribute in the comparison of
modal analysis results for both blades

2.5.1 3-D model for ideal troposkein shape blade
2.5.1.1 Blade profile
The troposkein shape shown in figure 2.18 can be described by [6]
𝑧0
0.5𝐻

𝐹(𝑘 ;𝜙)

= 1 − 𝐹(𝑘𝑒 ;𝜋)
𝑒 2

(2.2)

25

Figure 2.18 Troposkein shape

where
H is the blade height
𝑥

𝜙 = 𝐴𝑟𝑐𝑆𝑖𝑛[ 𝑅0 ], where 𝜙 ranges from 0 to

𝜋
2

, R is the blade radius

𝐹(𝑘𝑒 ; 𝜙) is the incomplete elliptic integral of the first kind
𝜋

𝐹(𝑘𝑒 ; 2 ) is the complete elliptic integral of the first kind
2.5.1.2 Elliptic integral arguments
𝑘𝑒 is elliptic integral arguments [6] and it can be calculated from
2𝑘

1

𝛽 = 1−𝑘𝑒2 ∗ 𝐹(𝑘

𝜋
𝑒 ;2)

(2.3)

𝑅

In Eq. 2.3 𝛽 = 𝐻 .
Solving Eq. (2.3) using Matlab for R=8.5m and H=24.476 gives 𝑘𝑒 =0.457
2.5.1.3 Blade length
The blade length can be calculated from 𝑘𝑒 and the blade height 2H [6] as in Equation 2.4

26

𝐿
2𝐻

2

𝜋
2
𝜋
𝑒; 2)

𝐸(𝑘𝑒 ; )

= 1−𝐾2 ∗ 𝐹(𝑘

(2.4)

−1

Solving Eq. 2.4 using MATLAB for 2H=24.476 m and k=0.457 gives L=30.42 m which is
approximately the same length of Sandia 17 m diameter simplified shape discussed in
sections 2.2 and 2.3. Plotting Eq. (2.2) for x0 ranging from 0 to R on MATLAB, where R=8.5
and H=12.238 m, gives half of the blade shape, as in figure 2.19. The 3-D blade model can be
generated by combining the blade cross section figure 2.17 and blade profile figure 2.19.

Figure 2.19 17 m diameter troposkein shape with H= 12.238 m

2.5.2 FEM modal analysis for 17m diameter blade with ideal troposkein shape
A modal analysis for 17 m diameter ideal troposkein shape blade has been done using
ANSYS WORKBENCH by meshing the blade structure using SOLID 187 element and
clamping both ends of the blade except for rotation. The blade assumed to be rotating at
constant rotational velocity of 54 rpm. The mesh pattern used here is the same as the mesh
pattern used for simplified Sandia blade shape to compare the results of both blades
significantly.
27

The first 10 modes and the corresponding natural frequencies of the blade have been
calculated and the results for the natural frequencies indicated in table 2.4 and the modes
in figures 2.20 to 2.29. The stationary ideal troposkein shaped blade was also analyzed and
the modal frequencies are listed in table 2.5, but the mode shapes are not shown for brevity
Table.2.4.first 10 natural frequencies of 17 m diameter ideal troposkein shape blade

Figure 2.20 first mode for the troposkein blade at 54 rpm (flat-wise only), 1.87 Hz

28

Figure 2.21 second mode for the troposkein blade at 54 rpm (edge-wise only), 2.34 Hz

Figure 2.22 third mode for the troposkein blade at 54 rpm (flat-wise only), 3.68 Hz

29

Figure 2.23 fourth mode for the troposkein blade at 54 rpm (flat-wise only), 5.82 Hz

Figure 2.24 fifth mode for the troposkein blade at 54 rpm (edge-wise & torsion), 6.2 Hz

30

Figure 2.25 sixth mode for the troposkein blade at 54 rpm (flat-wise only), 8.4 Hz

Figure 2.26 seventh mode for the troposkein blade at 54 rpm (flat-wise only), 11.36 Hz

31

Figure 2.27eighth mode for the troposkein blade at 54 rpm (edge-wise / torsion), 13.3 Hz

Figure 2.28 ninth mode for the troposkein blade at 54 rpm (flat-wise only), 14.82 Hz

32

Figure 2.29 tenth mode for the troposkein blade at 54 rpm (flat-wise only), 18.74 Hz

Table 2.5. Comparison of the lower modal frequencies as obtained by finite element
analysis. "Sandia" refers to the straight-circular-straight blade model. Spinning occurs at 54
rpm. F (flat-wise), E (edge-wise), and ET (edge-wise/ torsion)
Mode #

Sandia(54

Sandia (MSU)

Troposkein

Sandia (MSU)

Troposkein

rpm)

54rpm

54 rpm

stationary

(stationary

Freq.

shape

Freq.

shape

Freq.

shape

Freq.

shape

Freq.

shape

1

1.99

F

2.04

F

1.87

F

1.28

F

1.28

F

2

2.54

E

2.57

E

2.34

E

2.46

E

2.3

E

3

3.69

F

4

F

3.68

F

2.87

F

2.81

F

4

5.94

ET

5.94

ET

5.82

F

4.85

F

4.78

F

5

6.06

F

6.2

F

6.2

ET

5.62

ET

5.6

ET

6

8.47

F

8.54

F

8.36

F

7.12

F

7.21

F

7

11.72

F

11.75

F

11.36

F

10.2

F

10.14

F

8

13.51

ET

13.52

ET

13.3

ET

13.2

ET

13.06

ET

9

15.34

F

14.82

F

13.7

F

13.55

F

10

19.354

F

18.737

F

17.65

F

17.43

F

33

Table 2.5 summarizes our results for comparing the FEA results of Sandia blade and ideal
troposkein blade for both stationary and rotating condition with constant rotational
velocity 54 rpm as follows.
1. For the spinning case of the simplified blade we find that our finite element model
results and Sandia finite element results are essentially the same.
2. Comparing the spinning blade (54 rpm) first ten modes for ideal troposkein shape and
Sandia simplified shape we found that the natural frequencies spinning of ideal
troposkein shape is little less than the natural frequencies of Sandia one. This difference
may be because the ideal troposkein shape is slightly different, and its cross section is
little bit different than the simplified Sandia shape blade to overcome the limitations of
the number of elements in the ANSYS software academic version available in our lab,
which makes the ideal troposkein shaped blade is little lighter than the Sandia one.

3. For stationary free vibration case the first ten mode shapes for Sandia and ideal
troposkein shapes are essentially the same, but the natural frequencies of ideal
troposkein shape are a little lower than for simplified Sandia shape. This difference in
modal frequencies may be because of the difference in profile or difference in cross
section to overcome the ANSYS limitation problem as we mentioned before.
4. The spinning blades have increased modal frequencies. The increase is more significant
in the flat-wise modes than the edge-wise/torsion modes. It is likely that the spin
stiffening induces tension in the blade, and has a more significant effect on flat-wise
bending than torsion. Since the flat-wise modal frequencies undergo larger changes
under spin, in some cases their frequencies become larger than a nearby edge/torsion
mode, and the modal ordering therefore changes slightly.
34

3 Chapter 3: Analytical beam-based formulation of troposkein
shaped blade modal analysis
3.1 Introduction
In this chapter we will study the vibration of Darrieus wind turbine blade with a
troposkein-shape. Because of the blade shape complexity we will use an assumed modes
method for discretizing the kinetic energy and potential energy for the case of free
vibration (no aerodynamic forces condition). The gravitational effect will be neglected and
Lagrange’s equation will be used to formulate the mass and stiffness matrices which can
then be used to calculate the natural frequency and mode shapes analytically.
In the current work we will use the clamped-clamped beam modes as our assumed modal
functions to formulate the displacement functions for the blade, and by substituting these
functions and their derivatives we can get descretized expressions for the energies, which
can be used to approximate the modes and frequencies.

3.2 Blade potential energy
The potential energy of the blade can be divided in to two parts:
1. Strain energy
2. Gravitational potential energy or gravitational work
Here we will consider the strain energy only. The gravitational potential energy will be
neglected in this initial study. The strain energy SE for an element of volume V can be
expressed as
1

S𝐸 = 2 𝑉𝜎𝜖

(3.1)

where
35

V is the volume
σ is the engineering stress
𝜖 is the engineering strain
To express the blade strain energy using equation (3.1), several orthogonal coordinate
systems will be employed. These coordinate systems are shown in figures 3.1, 3.2, and 3.3
and can be defined as follows
1. Inertial coordinate system (XI YI ZI - system). As shown in figure 3.1 the ZI-axis coincides
with the vertical axis of the rotor. The XI-axis aligned with the free stream wind velocity,
𝑉∞ . The YI-axis is normal to XI-ZI plane.
2. Rotating coordinate system (XR YR ZR - system). As shown in figure 3.1, this system is
obtained by rotating the XI YI ZI – system by an angle τ =Ω*t, where Ω is the turbine
rotational velocity.
3. Blade coordinate system B1 (XB1YB1ZB1-system) ; as shown in figure 3.1 this blade local
coordinate system has its origin coincident with the blade cross section shear center
−

−

−

and parallel to the XR YR ZR – system. The unit vectors of this system are 𝑒𝑥B1 , 𝑒𝑦B1 , 𝑒𝑧B1 in
XB1 ,YB1 ,ZB1 respectively
4. Blade coordinate system B2 (XB2YB2ZB2-system). As shown in figure 3.1, this system is
obtained by rotating the B1 system about the negative YB1-axis by an angle θ (θ is the
meridian angle, which we will discuss later in detail). The unit vectors of this system are
−
−
−
𝑒𝑥B2 , 𝑒𝑦B2 , 𝑒𝑧B2

in XB2 , YB2 , ZB2, respectively.

5. Blade principle axes system B3 (XB3YB3ZB3-system). As shown in figure 3.2, in this
coordinate system XB3 and YB3 are taken to be aligned with the minor and major

36

principle axes of the blade cross section, respectively. The unit vectors of this system
−

−

−

are 𝑒𝑥B3 , 𝑒𝑦B3 , 𝑒𝑧B3 in XB3 , YB3, ZB3, respectively.
6. The blade coordinate system B6 (XB6YB6ZB6-system) is shown in figure 3.3 and is
obtained by translating and rotating the B3-system as we will see later in this chapter.
7. The deformations of the blade elastic axis are denoted by u, v, and w, in the XB3 ,YB3 , ZB3
respectively, and ϕ represents the twisting deformation.

Figure 3.1 vertical-axis wind turbine and coordinate systems [8]

Figure 3.2 coordinate systems of blade cross section [8]

37

Figure 3.3 axis of blade before and after deformation, and coordinate systems [8]

The blade strain energy can be expressed based on eq. (3.1) and using the B3 coordinate
system as
𝑆
1

𝑌3

𝑋3

𝑆𝐸 = 2 ∫ ∫ ∫ (𝜎𝑧3𝑧3 𝛾𝑧3𝑧3 + 𝜎𝑧3𝑥3 𝛾𝑧3𝑥3 + 𝜎𝑧3𝑦3 𝛾𝑧3 𝑦3 ) ⅆ𝑥3 ⅆ𝑦3 ⅆ𝑧3
0

0

0

(3.2)

Here because of the blade slenderness assumption we neglected 𝛾𝑥3 𝑥3 , 𝛾𝑥3 𝑦3 , and 𝛾𝑦3 𝑦3 .
Assuming that the engineering strain components 𝛾𝑧3 𝑧3 , 𝛾𝑧3 𝑥3 , 𝛾𝑧3 𝑦3 are equal to the
corresponding components of Lagrangian strain 𝜖𝑧3𝑧3 , 𝜖𝑧3 𝑥3 , 𝜖𝑧3 𝑦3 , and using Hooke’s law we
can get
𝜎𝑧3 𝑧3 = 𝐸𝛾𝑧3 𝑧3 = 𝐸𝜖𝑧3 𝑧3
(3.3)

𝜎𝑧3 𝑥3 = 𝐺𝛾𝑧3 𝑥3 = 2𝐺𝜖𝑧3𝑥3
38

𝜎𝑧3 𝑦3 = 𝐺𝛾𝑧3 𝑦3 = 2𝐺𝜖𝑧3𝑦3
Also
𝛾𝑧3 𝑧3 = 𝜖𝑧3 𝑧3
(3.4)

𝛾𝑧3 𝑥3 = 2𝜖𝑧3 𝑥3
𝛾𝑧3 𝑦3 = 2𝜖𝑧3 𝑦3
Substituting equations (3.3) and (3.4) into (3.2), we get
𝑆

𝑌3

1

𝑋3

𝑆𝐸 = ∫ ∫ ∫ (𝐸𝜖𝑧23 𝑧3 + 4𝐺(𝜖𝑧23𝑥3 + 𝜖𝑧23 𝑦3 )) ⅆ𝑥3 ⅆ𝑦3 ⅆ𝑧3
2

0

0

0

(3.5)

From eq. (3.5) we can find that the blade strain energy is only a function of Lagrangian
strain components 𝜖𝑧3𝑧3 , 𝜖𝑧3 𝑥3 , 𝜖𝑧3 𝑦3 . So we need to calculate these strain components for
the troposkein shaped blade.
3.2.1 Derivation of strain equations
From continuum mechanics we know that the Lagrangian finite strain 𝜖 for a beam of
length L subjected to axial stress is given by [20]
1 𝑙2 −𝐿2

𝜖 = 2(

𝐿2

(3.6)

)

where
l is the extended length
From eq. (3.6) we can express the strain tensor for the blade as
−

−

1 ❘𝑑𝑟1 ❘2 −❘𝑑𝑟0 ❘2

𝜖𝑖𝑗 = 2 (

−

−
❘𝑑𝑟0 ❘2

(3.7)

)

–

where 𝑟0 and 𝑟1 are the position vectors of an arbitrary mass point in the cross section of
the blade before and after deformation respectively [10] with respect to the origin of the
inertial coordinate system (XI YI ZI ) as shown in figure 3.1
39

Rearranging eq. 3.7 we get
−

−

−

❘ⅆ𝑟1 ❘ 2 − ❘ⅆ𝑟0 ❘ 2 = 2❘ⅆ𝑟0 ❘ 2 𝜖𝑖𝑗

(3.8)

Using figures 3.1, 3.2, and 3.3, and assuming that the blade cross section deformation is
negligible, using a slender blade assumption we get
−

−
𝑟0

= 𝑅 + 𝑥3 𝑒𝑥𝐵3 + 𝑦3 𝑒𝑦𝐵3

−

−
𝑟1

= 𝑅1 + 𝑥3 𝑒𝑥𝐵6 + 𝑦3 𝑒𝑦𝐵6

−

−

−

−

(3.9)

where
−

−

1. R and R1 are the position vectors of the blade cross section shear center (P) before and
after deformation with respect to the inertial coordinate system (XI YI ZI ) which is
shown in figure 3.3.
2. 𝑥3 and 𝑦3 are the coordinates of any arbitrary point in the blade cross section with
respect to the coordinate system (𝑋𝐵3 𝑌𝐵3 ) shown in figure 3.2, where in this system 𝑋3
and 𝑌3 are taken to be coincident with blade cross section minor and major axes,
−

−

respectively, and exB3 and eyB3 are the unit vectors in the direction of 𝑋𝐵3 and 𝑌𝐵3
respectively.
−

−

3. 𝑒𝑥𝐵6 and 𝑒𝑦𝐵6 are position vectors in the directions of 𝑋𝐵6 and 𝑌𝐵6 ,respectively, where
𝑋𝐵6 and 𝑌𝐵6 are the minor and major axes of the blade cross section, respectively, after
deformation. We will indicate later how we get the B6- coordinate system from B3coordinate system.
−

–

Now from eq. (3.8) we need the differentials of 𝑟0 and 𝑟1 . Firstly,
−
ⅆ𝑟0

−

𝑑𝑅

= [ 𝑑𝑠 + 𝑥3

−

𝑑𝑒𝑥𝐵3
𝑑𝑠

−

+ 𝑦3

𝑑𝑒𝑦𝐵3
𝑑𝑠

−

−

]ⅆ𝑠 + 𝑒𝑥𝐵3 ⅆ𝑥3 + 𝑒𝑦𝐵3 ⅆ𝑦3

40

(3.10)

where s is the same as Z3, and represents a coordinate tangent to the blade elastic axis at
each cross section and perpendicular to the cross section minor and major axes. The
derivative of the unit vectors with respect to the axial blade coordinate s can be written as
−

𝑑𝑒𝑥𝐵3
𝑑𝑠
−

𝑑𝑒𝑦𝐵3
𝑑𝑠

−

−

−

−

= 𝜔𝑋𝐵3 𝑌𝐵3 𝑍𝐵3 × 𝑒𝑥𝐵3
(3.11)

= 𝜔𝑋𝐵3 𝑌𝐵3 𝑍𝐵3 × 𝑒𝑦𝐵3
−

The curvature vector of the undeformed blade ωXB3YB3 ZB3 is given by
−

−

−

(3.12)

𝜔𝑋B3 𝑌B3 𝑍B3 = −𝜃′𝑒𝑦B2 + 𝛾′𝑒𝑧B3
where

𝜃 is the blade meridian angle or the angle between the vertical and the tangent to the blade
elastic axis at any point.
𝛾 is the total section pitch angle which includes the built-in twist and section pitch change
due to control inputs. In this thesis we will consider a blade with zero pitch angle.
From figure 3.2 the coordinate system B2 is the coordinate system generated by rotating
the B1-coordinate system by an angle 𝜃 about the negative yB1 -axis to make the zB1 -axis
tangent to the blade elastic axis. The transformation matrix between the B2 and B3
coordinate systems is
−
𝑒𝑥𝐵2
−
(𝑒𝑦𝐵2 )
−
𝑒𝑧𝐵2

−

𝑐𝑜𝑠𝛾
= [ 𝑠𝑖𝑛𝛾
0

−𝑠𝑖𝑛𝛾
𝑐𝑜𝑠𝛾
0

𝑒𝑥𝐵3

0 −
0] (𝑒𝑦𝐵3 )
1 𝑒−

(3.13)

𝑧𝐵3

Eq. (3.13) gives
−

−

−

(3.14)

𝑒𝑦B2 = 𝑒𝑥B3 sinγ + 𝑒𝑦B3 cosγ
Substituting eq. (3.14) in eq. (3.12) gives
41

−

−

−

−

(3.15)

𝜔𝑋B3 𝑌B3 𝑍B3 = −𝜃′sinγ𝑒𝑥B3 − 𝜃′cosγ𝑒𝑦B3 + 𝛾′𝑒𝑧B3
Now eq. (3.11) can be written as
−

de xB3
ds

−

−

de yB3
ds

−

= 𝛾′eyB3 + θ′cosγezB3
−

−

(3.16)

= −𝛾′exB3 − θ′sinγezB3

Substituting equation (3.16) into equation (3.10) and letting
𝑘𝑥B3 = −𝜃′sinγ
𝑘𝑦B3 = −𝜃′cosγ
𝑘𝑧B3 = 𝛾′
leads to
−

−

−

ⅆ𝑟0 = (dx3 − 𝑦3 𝑘𝑧B3 ds)𝑒𝑥B3 + (𝑥3 𝑘𝑧B3 ds + dy3 )𝑒𝑦B3
−

(3.17)

+(1 + 𝑦3 𝑘𝑥B3 − 𝑥3 𝑘𝑦B3 )ds𝑒𝑧B3
−

−

With the same procedure as ⅆ𝑟0 we can express ⅆ𝑟1 as

−

−

−

ⅆ𝑟1 = (ⅆ𝑥3 − 𝑦3 𝑘𝑧𝐵6 ⅆ𝑠1 )𝑒𝑥𝐵6 + (ⅆ𝑦3 + 𝑥3 𝑘𝑧𝐵6 ⅆ𝑠1 )𝑒𝑦𝐵6
−

(3.18)

+(1 − 𝑥3 𝑘𝑦𝐵6 + 𝑦3 𝑘𝑥𝐵6 )ⅆ𝑠1 𝑒𝑧𝐵6

Now from eq. (3.8) the Lagrangian strain tensor can be defined as follows
−
−
ⅆ𝑟1 . ⅆ𝑟1

−

−
−
ⅆ𝑟0 . ⅆ𝑟0

𝜖𝑥3 𝑥3
𝜖𝑖𝑗 = (𝜖𝑦3 𝑥3
𝜖𝑧3 𝑥3

ⅆ𝑥3
= 2[ⅆ𝑥3 ⅆ𝑥3 ⅆ𝑧3 ][𝜖𝑖𝑗 ] (ⅆ𝑥3 )
ⅆ𝑧3

𝜖𝑥3 𝑦3
𝜖𝑦3 𝑦3
𝜖𝑧3𝑦3

𝜖𝑥3 𝑧3
𝜖𝑦3 𝑧3 )
𝜖𝑧3 𝑧3

Expanding the right hand side of eq. (3.19) gives
42

(3.19)

ⅆ𝑥3
2[ⅆ𝑥3 ⅆ𝑥3 ⅆ𝑧3 ][𝜖𝑖𝑗 ] (ⅆ𝑥3 ) = 2(𝜖𝑥3 𝑥3 ⅆ𝑥3 2 + 𝜖𝑥3𝑦3 ⅆ𝑥3 ⅆ𝑦3 + 𝜖𝑥3 𝑧3 ⅆ𝑥3 ⅆ𝑧3 + 𝜖𝑦3 𝑦3 ⅆ𝑦3 2 +
ⅆ𝑧3
𝜖𝑦3 𝑥3 ⅆ𝑥3 ⅆ𝑦3 + 𝜖𝑦3 𝑧3 ⅆ𝑦3 ⅆ𝑧3 + 𝜖𝑧3𝑧3 ⅆ𝑧3 2 + 𝜖𝑧3 𝑥3 ⅆ𝑥3 ⅆ𝑧3 + 𝜖𝑧3𝑦3 ⅆ𝑧3 ⅆ𝑦3 )
Using the fact that the strain tensor is symmetric and neglecting any strain component in
the plane of the blade cross section (x3-y3) because of the blade slenderness, the above
equation reduces to
ⅆ𝑥3
2[ⅆ𝑥3 ⅆ𝑥3 ⅆ𝑧3 ][𝜖𝑖𝑗 ] (ⅆ𝑥3 ) = 4𝜖𝑧3 𝑥3 ⅆ𝑥3 ⅆ𝑧3 + 4𝜖𝑧3 𝑦3 ⅆ𝑦3 ⅆ𝑧3 + 2𝜖𝑧3𝑧3 ⅆ𝑧3 2
ⅆ𝑧3

(3.20)

By simplifying the left hand side of eq. (3.19)
−

−

−

−

2

2

ⅆ𝑟1 . ⅆ𝑟1 − ⅆ𝑟0 . ⅆ𝑟0 = {(ⅆ𝑥3 − 𝑦3 𝑘𝑧𝐵6 ⅆ𝑠1 ) + (ⅆ𝑦3 + 𝑥3 𝑘𝑧𝐵6 ⅆ𝑠1 ) + (1 − 𝑥3 𝑘𝑦𝐵6 +
2

2

2

𝑦3 𝑘𝑥𝐵6 ) ⅆ𝑠12 } − {(ⅆ𝑥3 − 𝑦3 𝑘𝑧𝐵3 ⅆ𝑠) + (ⅆ𝑦3 + 𝑥3 𝑘𝑧𝐵3 ⅆ𝑠) + (1 − 𝑥3 𝑘𝑦𝐵3 + 𝑦3 𝑘𝑥𝐵3 )2 ⅆ𝑠 2 }
(3.21)
The relation between ds1 and ds can be expressed as [10]
ⅆ𝑠1 = (1 + 2 𝜖𝑒 )0.5 ⅆ𝑠

(3.22)

where 𝜖𝑒 is the extensional component of the Green’s strain tensor
Substituting eq. (3.22) in to (3.21) and expanding leads to
−

−

−

−

ⅆ𝑟1 . ⅆ𝑟1 − ⅆ𝑟0 . ⅆ𝑟0 = −ⅆ𝑠 2 + ⅆ𝑠 2 (1 + 2𝜖𝑒 ) + 2ⅆ𝑠 2 𝑘𝑦𝐵3 𝑥3 − 2ⅆ𝑠 2 (1 + 2𝜖𝑒 )𝑘𝑦𝐵6 𝑥3 −
2ⅆ𝑠ⅆ𝑦3 𝑘𝑧𝐵3 𝑥3 + 2ⅆ𝑠(1 + 2𝜖𝑒 )ⅆ𝑦3 𝑘𝑧𝐵6 𝑥3 − ⅆ𝑠 2 𝑘𝑦2𝐵3 𝑥32 + ⅆ𝑠 2 (1 + 2𝜖𝑒 )𝑘𝑦2𝐵6 𝑥32 − ⅆ𝑠 2 𝑘𝑧2𝐵3 𝑥32 +
ⅆ𝑠 2 (1 + 2𝜖𝑒 )𝑘𝑧2𝐵6 𝑥32 − 2ⅆ𝑠 2 𝑘𝑥𝐵3 𝑦3 + 2ⅆ𝑠 2 (1 + 2𝜖𝑒 )𝑘𝑥𝐵6 𝑦3 + 2ⅆ𝑠ⅆ𝑥3 𝑘𝑧𝐵3 𝑦3 − 2ⅆ𝑠(1 +
2𝜖𝑒 )ⅆ𝑥3 𝑘𝑧𝐵6 𝑦3 + 2ⅆ𝑠 2 𝑘𝑥𝐵3 𝑘𝑦𝐵3 𝑥3 𝑦3 − 2ⅆ𝑠 2 (1 + 2𝜖𝑒 )𝑘𝑥𝐵6 𝑘𝑦𝐵6 𝑥3 𝑦3 − ⅆ𝑠 2 𝑘𝑥2𝐵3 𝑦32 + ⅆ𝑠 2 (1 +
2𝜖𝑒 )𝑘𝑥2𝐵6 𝑦32 − ⅆ𝑠 2 𝑘𝑧2𝐵3 𝑦32 + ⅆ𝑠 2 (1 + 2𝜖𝑒 )𝑘𝑧2𝐵6 𝑦32

43

Assuming that the extensional strain component 𝜖𝑒 is negligible compared to unity, we
have
−

−

−

−

ⅆ𝑟1 . ⅆ𝑟1 − ⅆ𝑟0 . ⅆ𝑟0 = 2ⅆ𝑠 2 𝜖𝑒 + 2ⅆ𝑠 2 𝑘𝑦𝐵3 𝑥3 − 2ⅆ𝑠 2 𝑘𝑦𝐵6 𝑥3 − 2ⅆ𝑠ⅆ𝑦3 𝑘𝑧𝐵3 𝑥3 +
2ⅆ𝑠ⅆ𝑦3 𝑘𝑧𝐵6 𝑥3 − ⅆ𝑠 2 𝑘𝑦2𝐵3 𝑥32 + ⅆ𝑠 2 𝑘𝑦2𝐵6 𝑥32 − ⅆ𝑠 2 𝑘𝑧2𝐵3 𝑥32 + ⅆ𝑠 2 𝑘𝑧2𝐵6 𝑥32 − 2ⅆ𝑠 2 𝑘𝑥𝐵3 𝑦3 +
2ⅆ𝑠 2 𝑘𝑥𝐵6 𝑦3 + 2ⅆ𝑠ⅆ𝑥3 𝑘𝑧𝐵3 𝑦3 − 2ⅆ𝑠ⅆ𝑥3 𝑘𝑧𝐵6 𝑦3 + 2ⅆ𝑠 2 𝑘𝑥𝐵3 𝑘𝑦𝐵3 𝑥3 𝑦3 − 2ⅆ𝑠 2 𝑘𝑥𝐵6 𝑘𝑦𝐵6 𝑥3 𝑦3 −
(3.23)

ⅆ𝑠 2 𝑘𝑥2𝐵3 𝑦32 + ⅆ𝑠 2 𝑘𝑥2𝐵6 𝑦32 − ⅆ𝑠 2 𝑘𝑧2𝐵3 𝑦32 + ⅆ𝑠 2 𝑘𝑧2𝐵6 𝑦32

Simplifying eq. (3.23) considering (dz3 = ds) along the blade elastic axis
−

−

−

−

ⅆ𝑟1 . ⅆ𝑟1 − ⅆ𝑟0 . ⅆ𝑟0 = 2(𝜖𝑒 + 𝑥3 (𝑘𝑦𝐵3 − 𝑘𝑦𝐵6 ) − 𝑦3 (𝑘𝑥𝐵3 − 𝑘𝑥𝐵6 ) + 𝑥3 𝑦3 (𝑘𝑥𝐵3 𝑘𝑦𝐵3 −
𝑘𝑥𝐵6 𝑘𝑦𝐵6 ) −

𝑦32
2

(𝑘𝑥2𝐵3 − 𝑘𝑥2𝐵6 ) +

𝑥32
2

(𝑘𝑦2𝐵6 − 𝑘𝑦2𝐵3 ) +

(𝑥32+𝑦32 )
2

(𝑘𝑧2𝐵6 − 𝑘𝑧2𝐵3 ))ⅆ𝑧3 2 − 2𝑥3 (𝑘𝑧𝐵3 −
(3.24)

𝑘𝑧𝐵6 )ⅆ𝑧3 ⅆ𝑦3 + 2𝑦3 (𝑘𝑧𝐵3 − 𝑘𝑧𝐵6 )ⅆ𝑧3 ⅆ𝑥3

Comparing the right hand side in eq. (3.20) to the left hand side in eq. (3.24) we find
ϵz3 z3 = ϵe + x3 (k yB3 − k yB6 ) − y3 (k xB3 − k xB6 ) + x3 y3 (k xB3 k yB3 − k xB6 k yB6 ) −
k 2xB6 ) +

x23
2

𝜖𝑧3 𝑥3 =

𝑦3

𝜖𝑧3 𝑦3 =

𝑥3

2

2

(k 2yB6 − k 2yB3 ) +

(x23 +y23 )
2

(k 2zB6 − k 2zB3 )

y23
2

(k 2xB3 −

(3.25)

(𝑘𝑧𝐵3 − 𝑘𝑧𝐵6 )

(3.26)

(𝑘𝑧𝐵3 − 𝑘𝑧𝐵6 )

(3.27)

For calculating the blade strain energy we need to get the strains in equations (3.25),
(3.26), and (3.27) in terms of the elastic deformations (u, v, and w) and the twisting
deformation ϕ, where u, v, and w are deformations in x3, y3, and z3 directions respectively.
As in Ref. [10] appendix A
1

(3.28)

𝜖𝑒 = 𝛼𝑧 + 2 (𝛼𝑥2 + 𝛼𝑦2 )
1

1

2

2

𝑘𝑥𝐵6 = 𝑘𝑥𝐵3 − 𝛼𝑥 𝑘𝑧𝐵3 − 𝛼𝑦 ′ − 𝑘𝑥𝐵3 𝛼𝑥2 − 𝑘𝑥𝐵3 𝜙 2 + 𝜙𝛼𝑥 ′ − 𝜙𝛼𝑦 𝑘𝑧𝐵3 + 𝜙𝑘𝑦𝐵3
44

1

𝑘𝑦𝐵6 = 𝑘𝑦𝐵3 − 𝛼𝑦 𝑘𝑧𝐵3 + 𝛼𝑥 ′ − 2 𝑘𝑦𝐵3 (𝛼𝑦2 + 𝜙 2 ) + 𝜙𝛼𝑦 ′ + 𝜙𝛼𝑥 𝑘𝑧𝐵3 − 𝜙𝑘𝑥𝐵3 − 𝑘𝑥𝐵3 𝛼𝑥 𝛼𝑦
(3.29)
1

𝑘𝑧𝐵6 = 𝛼𝑥 𝑘𝑥𝐵3 + 𝛼𝑦 𝑘𝑦𝐵3 + 𝑘𝑧𝐵3 − 2 𝑘𝑧𝐵3 (𝛼𝑦2 + 𝛼𝑥2 ) + 𝜙 ′ + 𝛼𝑥 ′ 𝛼𝑦
For case of zero section pitch angle [10]
𝑘𝑥B3 = 0
𝑘𝑦B3 = −𝜃 ′
(3.30)

𝑘𝑧B3 = 0
𝛼𝑥 = 𝑢′ − 𝑤𝜃 ′
𝛼𝑦 = 𝑣 ′
𝛼𝑧 = 𝑤 ′ + 𝑢𝜃 ′

Substituting equations (3.28), (3.29), and (3.30) in to equations (3.25), (3.26), (3.27) we get
𝜖𝑧3 𝑧3 =

(𝑣 ′)2
2

1

+ 𝑤 ′ + 𝑢𝜃 ′ + 2 (𝑢′ − 𝑤𝜃 ′ )2 + 𝑦3 (−𝜙𝜃 ′ + 𝜙(𝑢′′ − 𝑤′𝜃′ − 𝑤𝜃′′) − 𝑣 ′′ ) −

𝜃′

1

𝑥3 ( 2 (𝜙 2 + (𝑣 ′ )2 ) + (𝑢′′ − 𝑤′𝜃′ − 𝑤𝜃′′) + 𝜙𝑣 ′′ ) + 𝑥32 (− 2 𝜙 2 (𝜃 ′ )2 − 𝜃 ′ (𝑢′′ − 𝑤′𝜃′ − 𝑤𝜃′′) +
1
2

((𝑢′′ − 𝑤′𝜃′ − 𝑤𝜃′′))2 − 𝜙𝜃 ′ 𝑣 ′′ ) + 𝑥3 𝑦3 (−𝜙(𝜃 ′ )2 + 2𝜙𝜃 ′ (𝑢′′ − 𝑤′𝜃′ − 𝑤𝜃′′) − 𝜃 ′ 𝑣 ′′ +
1

1

(𝑢′′ − 𝑤′𝜃′ − 𝑤𝜃′′)𝑣 ′′ ) + 𝑦32 (2 𝜙 2 (𝜃 ′ )2 + 2 (𝑣 ′ )2 (𝜃 ′ )2 + 𝜙𝜃 ′ 𝑣 ′′ +

(𝑣 ′′ )2
2

)+

(𝑥32+𝑦32 )
2

((𝜙 ′ )2 −

(3.31)

𝑣′ 𝜃′ 𝜙′ )
1

𝜖𝑧3 𝑥3 = − 2 𝑦3 (−𝑣 ′ 𝜃 ′ + 𝜙 ′ + 𝑣 ′ (𝑢′′ − 𝑤′𝜃′ − 𝑤𝜃′′))
1

𝜖𝑧3 𝑦3 = 2 𝑥3 (−𝑣 ′ 𝜃 ′ + 𝜙 ′ + 𝑣 ′ (𝑢′′ − 𝑤′𝜃′ − 𝑤𝜃′′))

(3.32)
(3.33)

Now substituting equations (3.31), (3.32), and (3.33) in to equation (3.5), we can get the
strain energy as a function of u, v, w, and ϕ and their derivatives

45

𝑆
1

𝑌3

𝑋3

∫ ∫ {𝐸 [(

𝑆𝐸 = 2
∍
0

(𝑣 ′)

2

2

1

+ 𝑤 ′ + 𝑢𝜃 ′ + 2 (𝑢′ − 𝑤𝜃 ′ )2 + 𝑦3 (−𝜙𝜃 ′ + 𝜙(𝑢′′ − 𝑤 ′ 𝜃 ′ −

0

0

𝜃′

1

𝑤𝜃 ′′ ) − 𝑣 ′′ ) − 𝑥3 ( 2 (𝜙 2 + (𝑣 ′ )2 ) + (𝑢′′ − 𝑤 ′ 𝜃 ′ − 𝑤𝜃 ′′ ) + 𝜙𝑣 ′′ ) + 𝑥32 (− 2 𝜙 2 (𝜃 ′ )2 −
2

1

𝜃 ′ (𝑢′′ − 𝑤 ′ 𝜃 ′ − 𝑤𝜃 ′′ ) + ((𝑢′′ − 𝑤 ′ 𝜃 ′ − 𝑤𝜃 ′′ )) − 𝜙𝜃 ′ 𝑣 ′′ ) + 𝑥3 𝑦3 (−𝜙(𝜃 ′ )2 + 2𝜙𝜃 ′ (𝑢′′ −
2

1

1

𝑤 ′ 𝜃 ′ − 𝑤𝜃 ′′ ) − 𝜃 ′ 𝑣 ′′ + (𝑢′′ − 𝑤 ′ 𝜃 ′ − 𝑤𝜃 ′′ )𝑣 ′′ ) + 𝑦32 (2 𝜙 2 (𝜃 ′ )2 + 2 (𝑣 ′ )2 (𝜃 ′ )2 + 𝜙𝜃 ′ 𝑣 ′′ +
(𝑣 ′′ )
2

2

)+

(𝑥32+𝑦32 )
2

2

2

((𝜙

′ )2

′ ′

′)

1

− 𝑣 𝜃 𝜙 ) ] + 4𝐺 [(− 2 𝑦3 (−𝑣 ′ 𝜃 ′ + 𝜙 ′ + 𝑣 ′ (𝑢′′ − 𝑤 ′ 𝜃 ′ −
2

1

(3.34)

𝑤𝜃 ′′ ))) + (2 𝑥3 (−𝑣 ′ 𝜃 ′ + 𝜙 ′ + 𝑣 ′ (𝑢′′ − 𝑤 ′ 𝜃 ′ − 𝑤𝜃 ′′ ))) ]} ⅆ𝑥3 ⅆ𝑦3 ⅆ𝑧3

Because we seeking linear system of equations from Lagrange’s equation we will keep
quadratic terms only so we can simplify eq. (3.34) as follows:
𝑆
1

𝑌3

𝑋3

∫ ∫ {𝐸[(𝑤 ′ + 𝑢𝜃 ′ + 𝑦3 (−𝜙𝜃 ′ − 𝑣 ′′ ) − 𝑥3 (𝑢′′ − 𝑤′𝜃′ − 𝑤𝜃′′) + 𝑥32 (−𝜃 ′ (𝑢′′ −

𝑆𝐸 = 2
∍
0

0

0

2

1

1

𝑤′𝜃′ − 𝑤𝜃′′)) + 𝑥3 𝑦3 (−𝜙 (𝜃 ′ )2 − 𝜃 ′ 𝑣 ′′ ))2 ] + 4𝐺[(− 2 𝑦3 (−𝑣 ′ 𝜃 ′ + 𝜙 ′ )) + (2 𝑥3 (−𝑣 ′ 𝜃 ′ +
2

(3.35)

𝜙 ′ )) ]} ⅆ𝑥3 ⅆ𝑦3 ⅆ𝑧3
where

u, v, w, ϕ, and their derivatives represent u(s,t) , v(s,t) , w(s,t) , ϕ(s,t) and their derivatives.
Assuming that the blade cross section shown in figure 3.2 is symmetric about YB3-axis, the
cross sectional properties resulting from expanding eq. (3.35) are defined as follows
1. The blade cross section area A = ∍ ∍ dx3 dy3
46

2. AeA = ∍ ∍ y3 dx3 dy3 {eA is the chord-wise offset of area centeroid of cross section from
elastic axis (has a positive value when measured in front of elastic axis and vice versa)}
3. Edge-wise second moment of area of blade cross section Ix3 x3 = ∍ ∍ y3 2 dx3 dy3
4. Flat-wise second moment of area of blade cross section Iy3 y3 = ∍ ∍ x3 2 dx3 dy3
5. Polar second moment of area J = ∍ ∍ (y3 2 + x3 2 ) dx3 dy3
6. ∍ ∍ x3 dx3 dy3 = 0
7. ∍ ∍ x3 y3 dx3 dy3 = 0
8. ∍ ∍ x3 (x3 2 + y3 2 ) dx3 dy3 = 0
9. ∍ ∍ y3 x3 3 dx3 dy3 = 0
10. ∍ ∍ x3 y3 3 dx3 dy3 = 0
11. ∍ ∍ x3 y3 2 dx3 dy3 = 0
12. 𝑃1 = ∫ ∫ 𝑥3 4 d𝑥3 d𝑦3
13. 𝑃2 = ∫ ∫ 𝑥32 𝑦3 d𝑥3 d𝑦3
14. 𝑃3 = ∫ ∫ 𝑥32 𝑦32 d𝑥3 d𝑦3
Substituting the above blade cross sectional properties and expanding we get
1

𝑆

SE = 2 ∫ {𝐸[𝐴((𝑤 ′ )2 + 2𝑢𝑤 ′ 𝜃 ′ + 𝑢2 (𝜃 ′ )2 ) − 2𝐴𝑒𝐴 (𝜙𝑤 ′ 𝜃 ′ + 𝑢𝜙(𝜃 ′ )2 + 𝑤 ′ 𝑣 ′′ + 𝑢𝜃 ′ 𝑣 ′′ ) +
0

𝐼𝑥3 𝑥3 (𝜙 2 (𝜃 ′ )2 + 2𝜙𝜃 ′ 𝑣 ′′ + (𝑣 ′′ )2 ) + 𝑃3 (𝜙 2 (𝜃 ′ )4 + 2𝜙(𝜃 ′ )3 𝑣 ′′ + (𝜃 ′ )2 (𝑣 ′′ )2 ) +
4𝑃2 (−𝜙𝑤 ′ (𝜃 ′ )3 + 𝜙(𝜃 ′ )2 𝑢′′ − 𝑤 ′ (𝜃 ′ )2 𝑣 ′′ + 𝜃 ′ 𝑢′′ 𝑣 ′′ − 𝑤𝜙(𝜃 ′ )2 𝜃 ′′ − 𝑤𝜃 ′ 𝑣 ′′ 𝜃 ′′ ) +
𝐼𝑦3 𝑦3 (3(𝑤 ′ )2 (𝜃 ′ )2 + 2𝑢𝑤 ′ (𝜃 ′ )3 − 4𝑤 ′ 𝜃 ′ 𝑢′′ − 2𝑢(𝜃 ′ )2 𝑢′′ + (𝑢′′ )2 + 4𝑤𝑤 ′ 𝜃 ′ 𝜃 ′′ +
2𝑢𝑤(𝜃 ′ )2 𝜃 ′′ − 2𝑤𝑢′′ 𝜃 ′′ + 𝑤 2 (𝜃 ′′ )2 ) + 𝑃1 ((𝑤 ′ )2 (𝜃 ′ )4 − 2𝑤 ′ (𝜃 ′ )3 𝑢′′ + (𝜃 ′ )2 (𝑢′′ )2 +
2𝑤𝑤 ′ (𝜃 ′ )3 𝜃 ′′ − 2𝑤(𝜃 ′ )2 𝑢′′ 𝜃 ′′ + 𝑤 2 (𝜃 ′ )2 (𝜃 ′′ )2 )] + 𝐺𝐽[((𝑣 ′ )2 (𝜃 ′ )2 − 2𝑣 ′ 𝜃 ′ 𝜙 ′ + (𝜙 ′ )2 )]} d𝑠
(3.36)
47

3.3 Kinetic energy
−

The kinetic energy of the blade in terms of r 1 can be given by
𝑆

𝑌3

𝑋3

1

−

𝑑𝑟

𝑇 = 2 ∫ ∫ ∫ 𝜌( 𝑑𝑡1 .
0

0

0

−

𝑑𝑟1
𝑑𝑡

(3.37)

) ⅆ𝑥3 ⅆ𝑦3 ⅆ𝑧3 )

In eq. (3.36)
−

.

𝑑𝑟1

−

−

(3.38)

= 𝑟1 + 𝜔 × 𝑟1

𝑑𝑡

With reference to figure 3.3
−

𝜔 is the angular velocity of the B3-system.
−
𝑟1

is the position vector of an arbitrary mass point on the blade.
−

−

Note: 𝜔 and 𝑟1 are not shown in figure 3.3. Then
−
𝑟1

−

−

−

−

(3.39)

= 𝑅 + 𝛥𝑅 + 𝑥3 𝑒𝑥𝐵6 + 𝑦3 𝑒𝑦𝐵6

−

−

−

𝑅 = 𝑥0 𝑒𝑥𝐼 + 𝑧0 𝑒𝑧𝐼
For a nonrotating blade
−

−

−

𝑅 = 𝑥0 𝑒𝑥𝑅 + 𝑧0 𝑒𝑧𝑅
−

Because the R-system and B1-system are parallel we can express 𝑅 as
−

−

−

𝑅 = 𝑥0 𝑒𝑥𝐵1 + 𝑧0 𝑒𝑧𝐵1
The coordinate transformation between the B1 and B2 coordinate systems is
−

𝑒𝑥𝐵1

𝑐𝑜𝑠𝜃
−
(𝑒𝑦𝐵1 ) = [ 0
−
𝑠𝑖𝑛𝜃
𝑒
𝑧𝐵1

−

0
1
0

−𝑠𝑖𝑛𝜃 𝑒−𝑥𝐵2
0 ] (𝑒𝑦𝐵2 )
−
𝑐𝑜𝑠𝜃
𝑒

(3.40)

𝑧𝐵2

Using eq. (3.40)
48

−

−

−

𝑅 = (𝑥0 𝑐𝑜𝑠𝜃 + 𝑧0 𝑠𝑖𝑛𝜃)𝑒𝑥𝐵2 + (𝑧0 𝑐𝑜𝑠𝜃 − 𝑥0 𝑠𝑖𝑛𝜃)𝑒𝑧𝐵2
For zero section pitch angle B2 and B3 are the same so:
−

−

−

(3.41)

𝑅 = (𝑥0 𝑐𝑜𝑠𝜃 + 𝑧0 𝑠𝑖𝑛𝜃)𝑒𝑥𝐵3 + (𝑧0 𝑐𝑜𝑠𝜃 − 𝑥0 𝑠𝑖𝑛𝜃)𝑒𝑧𝐵3
With reference to figure 3.3
Note: u, v, w are not shown in figure 3.3
−

−

−

−

(3.42)

𝛥𝑅 = 𝑢𝑒𝑥𝐵3 + 𝑣𝑒𝑦𝐵3 + 𝑤𝑒𝑧𝐵3

The rotational transformation matrix [T] from B3 to B6 coordinates can be formulated by
Eulerian-type angles [8] β, θ, and Μ, which can be defined as follows:
1. β is the positive rotation angle about yB3 - axis which converts
𝑥𝐵3 𝑦𝐵3 𝑧𝐵3 to 𝑥𝐵4 𝑦𝐵4 𝑧𝐵4 (not shown)
2. ζ is the positive rotation angle about the negative 𝑥𝐵4 - axis which converts 𝑥𝐵4 𝑦𝐵4 𝑧𝐵4 to
𝑥𝐵5 𝑦𝐵5 𝑧𝐵5 (not shown)
3. σ is the positive rotation about zB5 - axis which converts 𝑥𝐵5 𝑦𝐵5 𝑧𝐵5 to 𝑥𝐵6 𝑦𝐵6 𝑧𝐵6 .
Applying these three rotations we get [T] as follows:
[𝑇] =
𝑐𝑜𝑠(𝛽)𝑐𝑜𝑠(𝜎) − 𝑠𝑖𝑛(𝜎)𝑠𝑖𝑛(𝜁)𝑠𝑖𝑛(𝛽) 𝑐𝑜𝑠(𝜁)𝑠𝑖𝑛(𝜎) −𝑐𝑜𝑠(𝜎)𝑠𝑖𝑛(𝛽) − 𝑠𝑖𝑛(𝜎)𝑠𝑖𝑛(𝜁)𝑐𝑜𝑠(𝛽)
[−𝑐𝑜𝑠(𝛽)𝑠𝑖𝑛(𝜎) − 𝑐𝑜𝑠(𝜎)𝑠𝑖𝑛(𝜁)𝑠𝑖𝑛(𝛽) 𝑐𝑜𝑠(𝜁)𝑐𝑜𝑠(𝜎) 𝑠𝑖𝑛(𝛽)𝑠𝑖𝑛(𝜎) − 𝑐𝑜𝑠(𝜎)𝑠𝑖𝑛(𝜁)𝑐𝑜𝑠(𝛽) ]
𝑐𝑜𝑠(𝜁)𝑠𝑖𝑛(𝛽)
𝑠𝑖𝑛(𝜁)
𝑐𝑜𝑠(𝛽)𝑐𝑜𝑠(𝜁)

(3.43)
Rewriting [T] matrix as a function of strain components u, v, w, and ϕ we get [10]
1−
[𝑇] =

𝜙2
2

1

− 2 (𝑢′ − 𝜃 ′ 𝑤)2

−𝜙 − (𝑢′ − 𝜃 ′ 𝑤)𝑣′
[

(𝑢′ − 𝜃′𝑤)

−(𝑢′ − 𝜃 ′ 𝑤) − 𝜙𝑣 ′

𝜙
1−

𝜙2
2

−

𝑣′2
2

𝜙(𝑢′ − 𝜃′𝑤) − 𝑣′
1

1 − 2 [(𝑢′ − 𝜃′𝑤)2 − 𝑣′2 ]]

𝑣′
49

(3.44)

The transformation from unit vectors in the B3- coordinate system to unit vectors in
B6- coordinate system can be written as
−

−

𝑒𝑥𝐵6

𝑒𝑥𝐵3

−

−

(3.45)

(𝑒𝑦𝐵6 ) = [𝑇] (𝑒𝑦𝐵3 )
−
𝑒𝑧𝐵6

−
𝑒𝑧𝐵3

Using eq. (3.45) with eq. (3.44)
−

𝑒𝑥𝐵6 = (1 −

𝜙2
2

−

𝛼𝑥 2 −
)𝑒𝑥𝐵3
2

−

−

−

+ 𝜙𝑒𝑦𝐵3 + (−𝛼𝑥 − 𝜙𝛼𝑦 )𝑒𝑧𝐵3
𝜙2

−

𝑒𝑦𝐵6 = (−𝜙 − 𝛼𝑥 𝛼𝑦 )𝑒𝑥𝐵3 + (1 −

2

−

𝛼𝑦 2 −
)𝑒𝑦𝐵3
2

−

(3.46)

+ (𝜙𝛼𝑥 − 𝛼𝑦 )𝑒𝑧𝐵3

Substituting equations (3.41), (3.42), and (3.46) in to equation 3.39 we get
−
𝑟1

𝜙2
2

= (𝑥0 𝑐𝑜𝑠𝜃 + 𝑧0 𝑠𝑖𝑛𝜃 + 𝑢 + 𝑥3 (1 −
−

𝛼𝑦 2
2

𝜙2

𝛼𝑥 2

−

2

2

−

) − 𝑦3 (𝜙 + 𝛼𝑥 𝛼𝑦 )) 𝑒𝑥𝐵3 + (𝑣 + 𝑥3 𝜙 + 𝑦3 (1 −

−

−

))𝑒𝑦𝐵3 + (𝑧0 𝑐𝑜𝑠𝜃 − 𝑥0 𝑠𝑖𝑛𝜃 + 𝑤 − 𝑥3 (𝛼𝑥 + 𝜙𝛼𝑦 ) + 𝑦3 (𝜙𝛼𝑥 − 𝛼𝑦 ))𝑒𝑧𝐵3

(3.47)

Substituting for Îąx and Îąy from eq. (3.30) in to eq. (3.47)
−
𝑟1

= [𝑥0 𝑐𝑜𝑠𝜃 + 𝑧0 𝑠𝑖𝑛𝜃 + 𝑢 + 𝑥3 (1 −
−

𝑣 ′ 𝑤𝜃 ′ ))] 𝑒𝑥𝐵3 + [𝑣 + 𝑥3 𝜙 + 𝑦3 (1 −

𝜙2

𝜙2
2

2

2

−

−

2

((𝑢 ′) −2𝑤𝑢 ′ 𝜃′ +𝑤 2(𝜃′ ) )

𝑣′

2
2

2

) − 𝑦3 (𝜙 + (𝑣 ′ 𝑢′ −

−

)] 𝑒𝑦𝐵3 + [𝑧0 𝑐𝑜𝑠𝜃 − 𝑥0 𝑠𝑖𝑛𝜃 + 𝑤 − 𝑥3 (𝑢′ −

−

(3.48)

𝑤𝜃 ′ + 𝜙𝑣 ′ ) + 𝑦3 (𝜙𝑢′ − 𝜙𝑤𝜃 ′ − 𝑣 ′ )]𝑒𝑧𝐵3
From eq. (3.48)
.

−
𝑟1

.

.

.

.

.

.

.

.

.

= [𝑢 + 𝑥3 (−𝜙𝜙 − 𝑢′ 𝑢′ − 𝜃 ′ (𝑤𝑢′ + 𝑤𝑢′ ) + 𝑤𝑤(𝜃 ′ )2 ) − 𝑦3 (𝜙 + 𝑣 ′ 𝑢′ + 𝑣 ′ 𝑢′ −
.

.

−

.

.

.

.

−

.

.

.

.

𝜃 ′ (𝑣 ′ 𝑤 + 𝑣 ′ 𝑤 ))] 𝑒𝑥𝐵3 + [𝑣 + 𝑥3 𝜙 + 𝑦3 (−𝜙𝜙 − 𝑣 ′ 𝑣 ′ )]𝑒𝑦𝐵3 + [𝑤 − 𝑥3 (𝑢′ − 𝑤 𝜃 ′ + 𝜙𝑣 ′ +
.

.

.

.

.

.

−

𝜙𝑣 ′ ) + 𝑦3 (𝜙𝑢′ + 𝜙𝑢′ − 𝜃 ′ (𝜙𝑤 + 𝜙𝑤) − 𝑣 ′ )]𝑒𝑧𝐵3

(3.49)

The angular velocity of the B3-coordinate system can be expressed as
50

−

−

−

−

(3.50)

𝜔 = 𝛺𝑒𝑧𝐼 = 𝛺𝑒𝑧𝑅 = 𝛺𝑒𝑧𝐵1
−

Substituting for ezB1 from equation (3.40) and for case of zero section pitch angle
−

−

−

(3.51)

𝜔 = 𝛺𝑠𝑖𝑛𝜃𝑒𝑥𝐵3 + 𝛺𝑐𝑜𝑠𝜃𝑒𝑧𝐵3
so
−

−

𝜔 × 𝑟1 = − [(𝑣 + 𝑥3 𝜙 + 𝑦3 (1 −

𝜙2
2

−

𝑣′

2

2

−

)) 𝛺𝑐𝑜𝑠𝜃] 𝑒𝑥𝐵3 − [(𝑧0 𝑐𝑜𝑠𝜃 − 𝑥0 𝑠𝑖𝑛𝜃 + 𝑤 −

𝑥3 (𝑢′ − 𝑤𝜃 ′ + 𝜙𝑣 ′ ) + 𝑦3 (𝜙𝑢′ − 𝜙𝑤𝜃 ′ − 𝑣 ′ ))𝛺𝑠𝑖𝑛𝜃 − (𝑥0 𝑐𝑜𝑠𝜃 + 𝑧0 𝑠𝑖𝑛𝜃 + 𝑢 +

𝑥3 (1 −

𝜙2
2

𝜙2

𝑦3 (1 −

2

2

−

−

2

((𝑢 ′) −2𝑤𝑢 ′𝜃′ +𝑤 2 (𝜃′ ) )
2

𝑣′

2

2

−

) − 𝑦3 (𝜙 + (𝑣 ′ 𝑢′ − 𝑣 ′ 𝑤𝜃 ′ ))) 𝛺𝑐𝑜𝑠𝜃] 𝑒𝑦𝐵3 + [(𝑣 + 𝑥3 𝜙 +

−

(3.52)

)) 𝛺𝑠𝑖𝑛𝜃]𝑒𝑧𝐵3

Substituting equations (3.49) and (3.52) in to equation 3.38 we find
−

𝑑𝑟1
𝑑𝑡

.

.

.

.

.

.

.

.

.

= [(𝑢 + 𝑥3 (−𝜙𝜙 − 𝑢′ 𝑢′ + 𝜃 ′ (𝑤 𝑢′ + 𝑤𝑢′ ) − 𝑤𝑤 (𝜃 ′ )2 ) − 𝑦3 (𝜙 + 𝑣 ′ 𝑢′ + 𝑣 ′ 𝑢′ −

.

.

𝜃 ′ (𝑣 ′ 𝑤 + 𝑣 ′ 𝑤 ))) − ((𝑣 + 𝑥3 𝜙 + 𝑦3 (1 −

.

𝜙2
2

−

𝑣′

2

2

−

.

.

)) 𝛺𝑐𝑜𝑠𝜃)] 𝑒𝑥𝐵3 + [(𝑣 + 𝑥3 𝜙 +

.

𝑦3 (−𝜙𝜙 − 𝑣 ′ 𝑣 ′ )) − ((𝑧0 𝑐𝑜𝑠𝜃 − 𝑥0 𝑠𝑖𝑛𝜃 + 𝑤 − 𝑥3 (𝑢′ − 𝑤𝜃 ′ + 𝜙𝑣 ′ ) + 𝑦3 (𝜙𝑢′ − 𝜙𝑤𝜃 ′ −

𝑣

′ ))

𝛺𝑠𝑖𝑛𝜃 + (𝑥0 𝑐𝑜𝑠𝜃 + 𝑧0 𝑠𝑖𝑛𝜃 + 𝑢 + 𝑥3 (1 −

𝜙2
2

51

2

−

2

((𝑢 ′ ) −2𝑤𝑢 ′𝜃′ +𝑤 2 (𝜃′ ) )
2

) − 𝑦3 (𝜙 +

−

.

.

.

.

.

.

.

(𝑣 ′ 𝑢′ − 𝑣 ′ 𝑤𝜃 ′ ))) 𝛺𝑐𝑜𝑠𝜃)] 𝑒𝑦𝐵3 + [(𝑤 − 𝑥3 (𝑢′ − 𝑤𝜃 ′ + 𝜙𝑣 ′ + 𝜙𝑣 ′ ) + 𝑦3 (𝜙𝑢′ + 𝜙𝑢′ −

.

.

.

𝜃 ′ (𝜙𝑤 + 𝜙𝑤 ) − 𝑣 ′ )) + ((𝑣 + 𝑥3 𝜙 + 𝑦3 (1 −

𝜙2
2

−

𝑣′

2

2

−

(3.53)

)) 𝛺𝑠𝑖𝑛𝜃)]𝑒𝑧𝐵3

Now we can get an expression for the kinetic energy as a function of u, v, w, and ϕ and the
inertial coordinates of the elastic axis by substituting eq. (3.53) into eq. (3. 37) to get
𝑆

𝑌3

𝑋3
.

.

1

.

∍
0
0 ∍
.

.

.

.

0
.

′ .

𝑣 ′ 𝑢′ + 𝑣 ′ 𝑢′ − 𝜃 ′ (𝑣 ′ 𝑤 + 𝑣 𝑤 ))) − ((𝑣 + 𝑥3 𝜙 + 𝑦3 (1 −

.

.

.

.

∫ {[(𝑢 + 𝑥3 (−𝜙𝜙 − 𝑢′ 𝑢′ − 𝜃 ′ (𝑤𝑢′ + 𝑤𝑢′ ) + 𝑤𝑤(𝜃 ′ )2 ) − 𝑦3 (𝜙 +

𝑇 = 2𝜌

.

𝜙2
2

−

𝑣′

2

2

2

)) 𝛺𝑐𝑜𝑠𝜃)] +

.

[(𝑣 + 𝑥3 𝜙 + 𝑦3 (−𝜙𝜙 − 𝑣 ′ 𝑣 ′ )) − ((𝑧0 𝑐𝑜𝑠𝜃 − 𝑥0 𝑠𝑖𝑛𝜃 + 𝑤 − 𝑥3 (𝑢′ − 𝑤𝜃 ′ + 𝜙𝑣 ′ ) +

𝑦3 (𝜙𝑢′ − 𝜙𝑤𝜃 ′ − 𝑣 ′ ))𝛺𝑠𝑖𝑛𝜃 − (𝑥0 𝑐𝑜𝑠𝜃 + 𝑧0 𝑠𝑖𝑛𝜃 + 𝑢 + 𝑥3 (1 −

𝜙2
2

−

2
2

2

((𝑢 ′ ) −2𝑤𝑢 ′𝜃′ +𝑤 2 (𝜃′ ) )
2

.

.

.

.

.

.

𝜙𝑣 ′ ) + 𝑦3 (𝜙𝑢′ + 𝜙𝑢′ − 𝜃 ′ (𝜙𝑤 + 𝜙𝑤) − 𝑣 ′ )) + ((𝑣 + 𝑥3 𝜙 + 𝑦3 (1 −

𝑣′

2

2

.

.

.

.

) − 𝑦3 (𝜙 + (𝑣 ′ 𝑢′ − 𝑣 ′ 𝑤𝜃 ′ ))) 𝛺𝑐𝑜𝑠𝜃)] + [(𝑤 − 𝑥3 (𝑢′ − 𝑤𝜃 ′ + 𝜙𝑣 ′ +

𝜙2
2

−

2

)) 𝛺𝑠𝑖𝑛𝜃)] } ⅆ𝑥3 ⅆ𝑦3 ⅆ𝑧3

(3.54)

Simplifying eq. (3.54) by keeping quadratic powers and neglecting any higher order terms
(linearize) and using the following blade cross sectional properties
52

1. The blade mass per unit length 𝑚 = ∫ ∫ 𝜌 d𝑥3 d𝑦3
2. me = ∫ ∫ ρy3 d𝑥3 d𝑦3 {e is the chord-wise offset of mass centeroid of cross section
from elastic axis (has a positive value when measured in front of elastic axis and vice
versa)}
3. Mass moment of inertia about Y3-axis 𝐼my3 = ∫ ∫ 𝜌𝑥32 d𝑥3 d𝑦3
4. Mass moment of inertia about X3-axis 𝐼mx3 = ∫ ∫ 𝜌𝑦32 d𝑥3 d𝑦3
5. ∫ ∫ 𝜌𝑥3 d𝑥3 d𝑦3 = 0
6. ∫ ∫ 𝜌𝑥3 𝑦3 d𝑥3 d𝑦3 = 0
Eq. 3.53 can be rewritten as
1

𝑆

.

.

. .

.

. .

.

. 2

.

. 2

. .

𝑇 = 2 ∫ {𝑚(𝑢2 + 𝑣 2 + 𝑤 2 ) + me(−2𝑢 𝜙 − 2𝑤𝑣 ′ ) + 𝐼𝑚𝑥3 (𝜙 2 + 𝑣 ′ ) + 𝐼𝑚𝑦3 (𝜙 2 + 𝑢′ − 2𝑤 𝑢′ 𝜃 ′ +
0

.

(3.55)

𝑤 2 (𝜃 ′ )2 )} d𝑠
.

∂

Note that in eq. (55); ( ) = ∂t

3.4 Discretization of energy expressions for modal analysis
In this section we will use modal functions of a uniform clamped-clamped beam as an
approximation for the assumed modal functions of the clamped-clamped Darrieus wind
turbine blade with troposkein shape.
For simplicity we will consider the Sandia blade to get expressions for mass and stiffness
matrices. The modal functions of a clamped-clamped beam [21] under the same
deformation conditions of the blade are:
1. flat-wise bending

53

𝑛

𝜓𝑢𝑗 (𝑠) = ∑𝑗=1{ 𝑐𝑜𝑠ℎ[

𝜆𝑗𝑠
𝑙

]−

𝜆𝑗 𝑠
]
𝑙

(𝑐𝑜𝑠ℎ[𝜆𝑗 ]−𝑐𝑜𝑠[𝜆𝑗 ])𝑠𝑖𝑛ℎ[
𝑠𝑖𝑛ℎ[𝜆𝑗 ]−𝑠𝑖𝑛[𝜆𝑗 ]

− 𝑐𝑜𝑠[

𝜆𝑗 𝑠
𝑙

]+

𝜆𝑗 𝑠
]
𝑙

(𝑐𝑜𝑠ℎ[𝜆𝑗 ]−𝑐𝑜𝑠[𝜆𝑗 ])𝑠𝑖𝑛[
𝑠𝑖𝑛ℎ[𝜆𝑗 ]−𝑠𝑖𝑛[𝜆𝑗 ]

}

(3.56)
2. Edge-wise bending
𝑛

𝜓𝑣𝑗 (𝑠) = ∑𝑗=1{ 𝑐𝑜𝑠ℎ[

𝜆𝑗𝑠
𝑙

]−

𝜆𝑗 𝑠
]
𝑙

(𝑐𝑜𝑠ℎ[𝜆𝑗 ]−𝑐𝑜𝑠[𝜆𝑗 ])𝑠𝑖𝑛ℎ[
𝑠𝑖𝑛ℎ[𝜆𝑗 ]−𝑠𝑖𝑛[𝜆𝑗 ]

− 𝑐𝑜𝑠[

𝜆𝑗 𝑠
𝑙

]+

𝜆𝑗 𝑠
]
𝑙

(𝑐𝑜𝑠ℎ[𝜆𝑗 ]−𝑐𝑜𝑠[𝜆𝑗 ])𝑠𝑖𝑛[
𝑠𝑖𝑛ℎ[𝜆𝑗 ]−𝑠𝑖𝑛[𝜆𝑗 ]

(3.57)
3. Extension
𝑛

𝜓𝑤𝑗 (𝑠) = ∑𝑗=1{ 𝑠𝑖𝑛[

𝑗𝜋𝑠
𝑙

]}

(3.58)

]}

(3.59)

4. Torsion
𝑛

𝜓𝜙𝑗 (𝑠) = ∑𝑗=1{ 𝑠𝑖𝑛[

𝑗𝜋𝑠
𝑙

We can express u, v, w, and ϕ as functions of the assumed modal functions ψ(s) and
assumed modal coordinates q (t) as follows
𝑢𝑗 (𝑠, 𝑡) = ∑𝑛𝑗=1{ 𝑞𝑢𝑗 (𝑡) 𝜓𝑢𝑗 (𝑠)}

(3.60)

𝑣𝑗 (𝑠, 𝑡) = ∑𝑛𝑗=1{ 𝑞𝑣𝑗 (𝑡) 𝜓𝑣𝑗 (𝑠)}

(3.61)

𝑤𝑗 (𝑠, 𝑡) = ∑𝑛𝑗=1{ 𝑞𝑤𝑗 (𝑡) 𝜓𝑤𝑗 (𝑠)}

(3.62)

𝜙𝑗 (𝑠, 𝑡) = ∑𝑛𝑗=1{ 𝑞𝜙𝑗 (𝑡) 𝜓𝜙𝑗 (𝑠)}

(3.63)

Substituting eqns. (3.56) - (3.59) in eqns. (3.60)-(3.63) we get
𝑢𝑗 (𝑠, 𝑡) = ∑𝑛𝑗=1{ 𝑞𝑢𝑗 (𝑡)(𝑐𝑜𝑠ℎ[
(𝑐𝑜𝑠ℎ[𝜆𝑗 ]−𝑐𝑜𝑠[𝜆𝑗 ])𝑠𝑖𝑛[
𝑠𝑖𝑛ℎ[𝜆𝑗 ]−𝑠𝑖𝑛[𝜆𝑗 ]

𝜆𝑗 𝑠
𝑙

]−

(𝑐𝑜𝑠ℎ[𝜆𝑗 ]−𝑐𝑜𝑠[𝜆𝑗 ])𝑠𝑖𝑛ℎ[
𝑠𝑖𝑛ℎ[𝜆𝑗 ]−𝑆𝑖𝑛[𝜆𝑗 ]

𝜆𝑗 𝑠
]
𝑙

} (3.64)

54

𝜆𝑗 𝑠
]
𝑙

− 𝑐𝑜𝑠[

𝜆𝑗 𝑠
𝑙

]+

}

𝑣𝑗 (𝑠, 𝑡) =

𝜆 𝑠
∑𝑛𝑗=1{ 𝑞𝑣𝑗 (𝑡)(𝑐𝑜𝑠ℎ[ 𝑗 ]
𝑙

−

(𝑐𝑜𝑠ℎ[𝜆𝑗 ]−𝑐𝑜𝑠[𝜆𝑗 ])𝑠𝑖𝑛ℎ[
𝑠𝑖𝑛ℎ[𝜆𝑗 ]−𝑠𝑖𝑛[𝜆𝑗 ]

𝜆𝑗 𝑠
]
𝑙

− 𝑐𝑜𝑠[

𝜆𝑗 𝑠
𝑙

]+

𝜆𝑗 𝑠
]
𝑙

(𝑐𝑜𝑠ℎ[𝜆𝑗 ]−𝑐𝑜𝑠[𝜆𝑗 ])𝑠𝑖𝑛[
𝑠𝑖𝑛ℎ[𝜆𝑗 ]−𝑠𝑖𝑛[𝜆𝑗 ]

}

(3.65)
𝑤𝑗 (𝑠, 𝑡) = ∑𝑛𝑗=1{ 𝑞𝑤𝑗 (𝑡)𝑠𝑖𝑛[
𝜙𝑗 (𝑠, 𝑡) = ∑𝑛𝑗=1{ 𝑞𝜙𝑗 (𝑡)𝑠𝑖𝑛[

𝑗𝜋𝑠
𝑙

𝑗𝜋𝑠
𝑙

]}

(3.65)

]}

(3.66)

3.4.1 Derivation of stiffness matrix expressions
Putting equations (3.60)-(3.63) into equation (3.36) we get
𝑆

𝑆𝐸 = ∑𝑛𝑗=1 ∑𝑛𝑘=1 0.5 ∫ {𝐴[𝐸𝑢𝑗 𝑢𝑘 (𝜃 ′ )2 + 𝐸𝑢𝑘 𝜃 ′ 𝑤𝑗 ′ + 𝐸𝑢𝑗 𝜃 ′ 𝑤𝑘 ′ + 𝐸𝑤𝑗 ′ 𝑤𝑘 ′ ] +
0

𝑃1 [𝐸(𝜃 ′ )4 𝑤𝑗 ′ 𝑤𝑘 ′ + 𝐸𝑤𝑘 (𝜃 ′ )3 𝑤𝑗 ′ 𝜃 ′′ + 𝐸𝑤𝑗 (𝜃 ′ )3 𝑤𝑘 ′ 𝜃 ′′ + 𝐸𝑤𝑗 𝑤𝑘 (𝜃 ′ )2 (𝜃 ′′ )2 − 𝐸(𝜃 ′ )3 𝑤𝑘 ′ 𝑢𝑗 ′′ −
𝐸𝑤𝑘 (𝜃 ′ )2 𝜃 ′′ 𝑢𝑗 ′′ − 𝐸(𝜃 ′ )3 𝑤𝑗 ′ 𝑢𝑘 ′′ − 𝐸𝑤𝑗 (𝜃 ′ )2 𝜃 ′′ 𝑢𝑘 ′′ + 𝐸(𝜃 ′ )2 𝑢𝑗 ′′ 𝑢𝑘 ′′ ] + 𝐴𝑒𝐴 [−𝐸𝑢𝑘 𝜙𝑗 (𝜃 ′ )2 −
𝐸𝑢𝑗 𝜙𝑘 (𝜃 ′ )2 − 𝐸𝜙𝑘 𝜃 ′ 𝑤𝑗 ′ − 𝐸𝜙𝑗 𝜃 ′ 𝑤𝑘 ′ − 𝐸𝑢𝑘 𝜃 ′ 𝑣𝑗 ′′ − 𝐸𝑤𝑘 ′ 𝑣𝑗 ′′ − 𝐸𝑢𝑗 𝜃 ′ 𝑣𝑘 ′′ − 𝐸𝑤𝑗 ′ 𝑣𝑘 ′′ ] +
𝐼𝑥3 𝑥3 [𝐸𝜙𝑗 𝜙𝑘 (𝜃 ′ )2 + 𝐺(𝜃 ′ )2 𝑣𝑗 ′ 𝑣𝑘 ′ − 𝐺𝜃 ′ 𝑣𝑘 ′ 𝜙𝑗 ′ − 𝐺𝜃 ′ 𝑣𝑗 ′ 𝜙𝑘 ′ + 𝐺𝜙𝑗 ′ 𝜙𝑘 ′ + 𝐸𝜙𝑘 𝜃 ′ 𝑣𝑗 ′′ +
𝐸𝜙𝑗 𝜃 ′ 𝑣𝑘 ′′ + 𝐸𝑣𝑗 ′′ 𝑣𝑘 ′′ ] + 𝑝3 [𝐸𝜙𝑗 𝜙𝑘 (𝜃 ′ )4 + 𝐸𝜙𝑘 (𝜃 ′ )3 𝑣𝑗 ′′ + 𝐸𝜙𝑗 (𝜃 ′ )3 𝑣𝑘 ′′ + 𝐸(𝜃 ′ )2 𝑣𝑗 ′′ 𝑣𝑘 ′′ ] +
𝐼𝑦3 𝑦3 [𝐺(𝜃 ′ )2 𝑣𝑗 ′ 𝑣𝑘 ′ + 𝐸𝑢𝑘 (𝜃 ′ )3 𝑤𝑗 ′ + 𝐸𝑢𝑗 (𝜃 ′ )3 𝑤𝑘 ′ + 3𝐸(𝜃 ′ )2 𝑤𝑗 ′ 𝑤𝑘 ′ − 𝐺𝜃 ′ 𝑣𝑘 ′ 𝜙𝑗 ′ −
𝐺𝜃 ′ 𝑣𝑗 ′ 𝜙𝑘 ′ + 𝐺𝜙𝑗 ′ 𝜙𝑘 ′ + 𝐸𝑢𝑘 𝑤𝑗 (𝜃 ′ )2 𝜃 ′′ + 𝐸𝑢𝑗 𝑤𝑘 (𝜃 ′ )2 𝜃 ′′ + 2𝐸𝑤𝑘 𝜃 ′ 𝑤𝑗 ′ 𝜃 ′′ + 2𝐸𝑤𝑗 𝜃 ′ 𝑤𝑘 ′ 𝜃 ′′ +
𝐸𝑤𝑗 𝑤𝑘 (𝜃 ′′ )2 − 𝐸𝑢𝑘 (𝜃 ′ )2 𝑢𝑗 ′′ − 2𝐸𝜃 ′ 𝑤𝑘 ′ 𝑢𝑗 ′′ − 𝐸𝑤𝑘 𝜃 ′′ 𝑢𝑗 ′′ − 𝐸𝑢𝑗 (𝜃 ′ )2 𝑢𝑘 ′′ − 2𝐸𝜃 ′ 𝑤𝑗 ′ 𝑢𝑘 ′′ −
𝐸𝑤𝑗 𝜃 ′′ 𝑢𝑘 ′′ + 𝐸𝑢𝑗 ′′ 𝑢𝑘 ′′ ] + 𝑃2 [−2𝐸𝜙𝑘 (𝜃 ′ )3 𝑤𝑗 ′ − 2𝐸𝜙𝑗 (𝜃 ′ )3 𝑤𝑘 ′ − 2𝐸𝑤𝑘 𝜙𝑗 (𝜃 ′ )2 𝜃 ′′ −
2𝐸𝑤𝑗 𝜙𝑘 (𝜃 ′ )2 𝜃 ′′ + 2𝐸𝜙𝑘 (𝜃 ′ )2 𝑢𝑗 ′′ + 2𝐸𝜙𝑗 (𝜃 ′ )2 𝑢𝑘 ′′ − 2𝐸(𝜃 ′ )2 𝑤𝑘 ′ 𝑣𝑗 ′′ − 2𝐸𝑤𝑘 𝜃 ′ 𝜃 ′′ 𝑣𝑗 ′′ +
2𝐸𝜃 ′ 𝑢𝑘 ′′ 𝑣𝑗 ′′ − 2𝐸(𝜃 ′ )2 𝑤𝑗 ′ 𝑣𝑘 ′′ − 2𝐸𝑤𝑗 𝜃 ′ 𝜃 ′′ 𝑣𝑘 ′′ + 2𝐸𝜃 ′ 𝑢𝑗 ′′ 𝑣𝑘 ′′ ]} ⅆ𝑠

(3.67)

Combining eqns. (3.60-3.63) with eq. (3.67) and simplifying using the blade cross sectional
properties, we can express the stiffness matrices equations as follows
1. Flat-wise only
55

𝑢𝑢
𝐾𝑗𝑘

𝑆

1

= 2 𝐸 ∫ {𝐴𝜃 ′ [𝑠]2 𝜓𝑢𝑗 [𝑠]𝜓𝑢𝑘 [𝑠] + 𝐼𝑦3 𝑦3 (−2𝜃 ′ [𝑠]2 𝜓𝑢𝑗 [𝑠]𝜓𝑢𝑘 ′′ [𝑠] + 𝜓𝑢𝑗 ′′ [𝑠]𝜓𝑢𝑘 ′′ [𝑠]) +
0

𝑃1 𝜃 ′ [𝑠]2 𝜓𝑢𝑗 ′′ [𝑠]𝜓𝑢𝑘 ′′ [𝑠]} ⅆ𝑠

(3.68)

2. Edge-wise only
S

1

1

′′
′′
′′
′′
′
′
′
2
′
2
K vv
jk = ∍ (2 E(Ix3 x3 ψvj [s]ψvk [s] + P3 θ [s] ψvj [s]ψvk [s]) + 2 G(J)θ [s] ψvj [s]ψvk [s]) ds
0

(3.69)
3. Extension only
𝑆

1

𝑤𝑤
𝐾𝑗𝑘
= 2 𝐸 ∫ (𝐴𝜓𝑤𝑗 ′ [𝑠]𝜓𝑤𝑘 ′ [𝑠] + 𝐼𝑦3 𝑦3 (𝜃 ′′ [𝑠]2 𝜓𝑤𝑗 [𝑠]𝜓𝑤𝑘 [𝑠] + 4𝜃 ′ [𝑠]𝜃 ′′ [𝑠]𝜓𝑤𝑗 [𝑠]𝜓𝑤𝑘 ′ [𝑠] +
0

3𝜃 ′ [𝑠]2 𝜓𝑤𝑗 ′ [𝑠]𝜓𝑤𝑘 ′ [𝑠]) + 𝑃1 (𝜃 ′ [𝑠]2 𝜃 ′′ [𝑠]2 𝜓𝑤𝑗 [𝑠]𝜓𝑤𝑘 [𝑠] + 2𝜃 ′ [𝑠]3 𝜃 ′′ [𝑠]𝜓𝑤𝑗 [𝑠]𝜓𝑤𝑘 ′ [𝑠] +
𝜃 ′ [𝑠]4 𝜓𝑤𝑗 ′ [𝑠]𝜓𝑤𝑘 ′ [𝑠])) ⅆ𝑠

(3.70)

4. Torsion only
𝑆

𝜙𝜙

1

1

𝐾𝑗𝑘 = ∫ (2 𝐸(𝐼𝑥3𝑥3 𝜃 ′ [𝑠]2 𝜓𝜙𝑗 [𝑠]𝜓𝜙𝑘 [𝑠] + 𝑃3 𝜃 ′ [𝑠]4 𝜓𝜙𝑗 [𝑠]𝜓𝜙𝑘 [𝑠]) + 2 𝐺𝐽𝜓𝜙𝑗 ′ [𝑠]𝜓𝜙𝑘 ′ [𝑠]) ⅆ𝑠
0

(3.71)
5. Flat wise-Edge wise
𝑆
𝑢𝑣
𝐾𝑗𝑘
=∍

1
2

0

𝐸1(−2𝐴𝑒𝐴 𝜃 ′ [𝑠]𝜓𝑢𝑗 [𝑠]𝜓𝑣𝑘 ′′ [𝑠] + 4𝑃2 𝜃 ′ [𝑠]𝜓𝑢𝑗 ′′ [𝑠]𝜓𝑣𝑘 ′′ [𝑠]) ⅆ𝑠 (3.72)

6. Flat wise – extension
𝑢𝑤
𝐾𝑗𝑘

1

𝑆

= 2 𝐸 ∫ {𝐴(2𝜃 ′ [𝑠]𝜓𝑢𝑗 [𝑠]𝜓𝑤𝑘 ′ [𝑠]) + 𝐼𝑦3 𝑦3 (2𝜃 ′ [𝑠]2 𝜃 ′′ [𝑠]𝜓𝑢𝑗 [𝑠]𝜓𝑤𝑘 [𝑠] +
0

2𝜃 ′ [𝑠]3 𝜓𝑢𝑗 [𝑠]𝜓𝑤𝑘 ′ [𝑠] − 2𝜃 ′′ [𝑠]𝜓𝑤𝑘 [𝑠]𝜓𝑢𝑗 ′′ [𝑠] − 4𝜃 ′ [𝑠]𝜓𝑤𝑘 ′ [𝑠]𝜓𝑢𝑗 ′′ [𝑠]) +

56

𝑃1 (−2𝜃 ′ [𝑠]2 𝜃 ′′ [𝑠]𝜓𝑤𝑘 [𝑠]𝜓𝑢𝑗 ′′ [𝑠] − 2𝜃 ′ [𝑠]3 𝜓𝑤𝑘 ′ [𝑠]𝜓𝑢𝑗 ′′ [𝑠])} ⅆ𝑠
(3.73)
7. Flat wise-torsion
𝑆

1

𝑢𝜙

𝐾𝑗𝑘 = 2 𝐸 ∫ {−2𝐴𝑒𝐴 𝜃 ′ [𝑠]2 𝜓𝑢𝑗 [𝑠]𝜓𝜙𝑘 [𝑠] + 4𝑃2 𝜃 ′ [𝑠]2 𝜓𝜙𝑘 [𝑠]𝜓𝑢𝑗 ′′ [𝑠]} ⅆ𝑠
0

8. Edge wise-extension
𝑣𝑤
𝐾𝑗𝑘

𝑆

1

= 2 𝐸 ∫ {−2𝐴𝑒𝐴 𝜓𝑤𝑘 ′ [𝑠]𝜓𝑣𝑗 ′′ [𝑠] + 𝑃2 (−4𝜃 ′ [𝑠]𝜃 ′′ [𝑠]𝜓𝑤𝑘 ′ [𝑠]𝜓𝑣𝑗 ′′ [𝑠] −
0

4𝜃 ′ [𝑠]2 𝜓𝑤𝑘 ′ [𝑠]𝜓𝑣𝑗 ′′ [𝑠])} ⅆ𝑠

(3.75)

9. Edge wise-torsion
𝑣𝜙

𝑆

𝐾𝑗𝑘 = ∫ {𝐸(𝐼𝑥3𝑥3 𝜃 ′ [𝑠]𝜓𝜙𝑘 [𝑠]𝜓𝑣𝑗 ′′ [𝑠] + 𝑃3 𝜃 ′ [𝑠]3 𝜓𝜙𝑘 [𝑠]𝜓𝑣𝑗 ′′ [𝑠]) −
0

𝐺𝐽𝜃 ′ [𝑠]𝜓𝑣𝑗 ′ [𝑠]𝜓𝑤𝑘 ′ [𝑠]} ⅆ𝑠 (3.76)
10. Extension-torsion
𝑤𝜙
𝐾𝑗𝑘

𝑆

= 𝐸 ∫ {−𝐴𝑒𝐴 𝜃 ′ [𝑠]𝜓𝜙𝑘 [𝑠]𝜓𝑤𝑗 ′ [𝑠] + 𝑃2 (−2𝜃 ′ [𝑠]2 𝜃 ′′ [𝑠]𝜓𝑤𝑗 [𝑠]𝜓𝜙𝑘 [𝑠] −
0

2𝜃 ′ [𝑠]3 𝜓𝜙𝑘 [𝑠]𝜓𝑤𝑗 ′ [𝑠])} ⅆ𝑠

(3.77)

Note that the stiffness matrix of the blade is symmetric so we have
𝑢𝑣
𝑣𝑢
𝐾𝑗𝑘
= 𝐾𝑗𝑘
𝑢𝑤
𝑤𝑢
𝐾𝑗𝑘
= 𝐾𝑗𝑘
𝑢𝜙

𝜙𝑢

𝐾𝑗𝑘 = 𝐾𝑗𝑘

(3.78)

𝑣𝑤
𝑤𝑣
𝐾𝑗𝑘
= 𝐾𝑗𝑘
𝑣𝜙

𝜙𝑣

𝐾𝑗𝑘 = 𝐾𝑗𝑘

57

(3.74)

𝑤𝜙

𝜙𝑤

𝐾𝑗𝑘 = 𝐾𝑗𝑘

3.4.2 Derivation of mass matrix expressions
In this section the mass matrix expressions were derived for the nonrotating blade
condition (for simplicity Ω = 0). For u, v, w, and ϕ as expressed in eqns. (3.60-3.63) we can
rewrite eq. (3.55) as
𝑛

𝑆

𝑛

1

𝑇 = 2𝜌∑

∑
𝑗=1

𝑋3

.

.

∫ ∫ ∫ {𝜓𝑢𝑗 [𝑠]𝜓𝑢𝑘 [𝑠]𝑞𝑢𝑗 [𝑡]𝑞𝑢𝑘 [𝑡] +
𝑘=1

.

𝑌3

0

0

0

.

.

.

.

.

𝜓𝑣𝑗 [𝑠]𝜓𝑣𝑘 [𝑠]𝑞𝑣𝑗 [𝑡]𝑞𝑣𝑘 [𝑡] + 𝜓𝑤𝑗 [𝑠]𝜓𝑤𝑘 [𝑠]𝑞𝑤𝑗 [𝑡]𝑞𝑤𝑘 [𝑡] + 𝑥32 𝜓𝜙𝑗 [𝑠]𝜓𝜙𝑘 [𝑠]𝑞𝜙𝑗 [𝑡]𝑞𝜙𝑘 [𝑡] +
.

.

.

.

.

.

𝑥32 𝜓𝑢𝑗 ′ [𝑠]𝜓𝑢𝑘 ′ [𝑠]𝑞𝑢𝑗 [𝑡]𝑞𝑢𝑘 [𝑡] − 2𝑦3 𝜓𝑢𝑗 [𝑠]𝜓𝜙𝑘 [𝑠]𝑞𝑢𝑗 [𝑡]𝑞𝜙𝑘 [𝑡] −
.

.

2𝑦3 𝜓𝑣𝑗 ′ [𝑠]𝜓𝑤𝑘 [𝑠]𝑞𝑣𝑗 [𝑡]𝑞𝑤𝑘 [𝑡] + 𝑦32 𝜓𝜙𝑗 [𝑠]𝜓𝜙𝑘 [𝑠]𝑞𝜙𝑗 [𝑡]𝑞𝜙𝑘 [𝑡] +
.

.

.

.

𝑦32 𝜓𝑣𝑗 ′ [𝑠]𝜓𝑣𝑘 ′ [𝑠]𝑞𝑣𝑗 [𝑡]𝑞𝑣𝑘 [𝑡] − 2𝑥32 𝜃 ′ [𝑠]𝜓𝑢𝑗 ′ [𝑠]𝜓𝑤𝑘 [𝑠]𝑞𝑢𝑗 [𝑡]𝑞𝑤𝑘 [𝑡] +
.

.

(3.79)

𝑥32 (𝜃 ′ [𝑠])2 𝜓𝑤𝑗 [𝑠]𝜓𝑤𝑘 [𝑠]𝑞𝑤𝑗 [𝑡]𝑞𝑤𝑘 [𝑡]} ⅆ𝑥3 ⅆ𝑦3 ⅆ𝑧3

From the kinetic energy expression shown in eq. (3.79) we can formulate the mass matrix
as follows
1. Flat-wise only
S

1

′
′
muu
jk = 2 ∍ {mψuj [s]ψuk [s] + Imy3 ψuj [s]ψuk [s]} ds

(3.80)

0

2. Edge-wise only
1

S

′
′
mvv
jk = 2 ∍ {mψvj [s]ψvk [s] + Imx3 ψvj [s]ψvk [s]} ds
0

3. Extension only

58

(3.81)

1

S

(3.82)

′
2
mww
jk = 2 ∍ {mψwj [s]ψwk [s] + Imy3 (θ [s]) ψwj [s]ψwk [s]} ds
0

4. Torsion only
𝑆

1

ϕϕ

(3.83)

𝑚jk = 2 (𝐼mx3 + 𝐼my3 ) ∫ {𝜓𝜙𝑗 [𝑠]𝜓𝜙𝑘 [𝑠]} d𝑠
0

5. Flat wise- edge wise
(3.84)

muv
jk = 0
6. Flat wise-extension
S
′
′
muw
jk = −Imy3 ∫ {θ [s]ψuj [s]ψwk [s]} ds

(3.85)

0

7. Flat wise-torsion
S

uϕ

(3.86)

mjk = −me ∫ { ψuj [s]ψϕk [s]} ds
0

8. Edge wise-extension
S

mvw
jk

= −me ∫ {ψvj ′ [s]ψwk [s]} ds

(3.87)

0

9. Edge wise-torsion
vϕ

(3.88)

𝑚jk = 0
10. Extension-torsion
wϕ

(3.89)

mjk = 0
Note that the mass matrix of the blade is symmetric so we have
vu
muv
jk = mjk
wu
muw
jk = mjk
uϕ

ϕu

mjk = mjk

59

(3.90)

wv
mvw
jk = mjk
vϕ

ϕv

wϕ

ϕw

mjk = mjk

mjk = mjk

3.4.3 Meridian angle (θ)
In the mass and stiffness expressions derived in section 3.4.2 we can see the parameter θ
and its derivatives with respect to z3 or s. This θ is known as the meridian angle and is
defined as the angle between the tangent to the blade at any point in the elastic axis and the
vertical axis at this point as shown in figure 2.17. For the simplified Sandia shape the
meridian angle can be calculated as:
S

θ = Rc

c

where
Sc: is the length of the circular segment of Sandia simplified shape blade
Rc: is the radius of the circular segment of Sandia simplified shape blade
θ′ =
θ′′ =

dθ
ds

1

=R

d2 θ
ds

c

=0

3.5 Numerical calculation of modal frequencies of 17 m diameter Sandia
simplified blade shape
In this section the analytical model developed in sections 3.2, 3.3, 3.4, was tested by using
the derived modal mass and stiffness matrices expressions to calculate the first ten modal
frequencies for 17 m diameter stationary simplified Sandia blade and comparing it with the
finite element modal analysis results.
60

3.5.1 17 m diameter simplified Sandia blade cross section properties
The blade cross section geometric and mass related parameters are calculated here for the
NASA 0015 blade cross section shown in figure 2.2 as follows
A = 0.98 ∗ 10−2 m2
Iy3 y3 = 0.859 ∗ 10−5 m4
AeA = −0.001434m3
Ix3 x3 = 2.814 ∗ 10−4 m4
P1 = 2.897 ∗ 10−9 m6
P2 = −6.8499 ∗ 10−7 m5
P3 = 1.816 ∗ 10−7 m6
Îť1 =4.730040744862704
Îť2 =7.853204624095838
Îť3 =10.995607838001671
Îť4 =14.137165491257464
Îť5 =17.27875965739948
Îť6 =20.42035224562606
Îť7 =23.561944902040455
E = 68.95 ∗ 109 N/m2
G = 25.86 ∗ 109 N/m2
m = 30.5kg/m
me = −3.97Kg
Imx3 = 1.5 kg. m
Imy3 = 0.0247 kg. m
61

Blade length S = 30.47m
Meridian angle θ= 46.36o for straight segments
Meridian angle θ is variable for circular segments
θ'= θ''=0 for straight segments
θ'= 0.096 rad/m for circular segment
θ''=0 for circular segment
3.5.2 Numerical values for first ten modal frequencies
Substituting the cross section parameters values calculated in section 3.5.1 in addition to
the blade length (S) and the meridian angle (θ) in the modal mass and stiffness submatrices equations and varying the summation parameters values (j and k) we got the best
convergence at ( j=k=7). When j=k=7, and all flat wise, edge wise, torsion, and extension
modes included, modal mass and stiffness gives the first 28 modal frequencies as shown in
table 3.1
Table 3.1 modal frequencies for 17 m diameter Sandia stationary blade (coupling between
all modes considered
modal frequencies (1 to 14)

modal frequencies (15 to 28)

1.5564

83.2925

3.00157

93.1968

4.32495

124.232

8.4747

134.915

10.254

162.029

11.7791

180.423

62

Table 3.1 (cont’d)
13.9191

218.931

15.5556

236.816

20.0909

263.703

23.5478

310.944

27.9191

327.779

49.2278

388.544

50.7068 i

462.103

55.4277

540.848

From table 3.1 we can see that the lowest seven modal frequencies are of a similar order of
magnitude as the first seven modal frequencies calculated using the FE model in chapter 2
with around 70% so the error here is around 30% which is not acceptable. Also the modal
frequency number 13 in red color is complex and this is not allowed because the mass and
stiffness matrices are supposed to be positive definite, and symmetric. The K matrix
contains some asymmetric terms. The terms are small compared to other terms in the K
matrix. They may be due to numerical errors, but in a continuing investigation we are
examining this further.
Because we followed a good derivation procedure to get expressions for kinetic and
potential energies and we checked that these equations are physically and dimensionally
correct, the sources of results inconsistency with our finite element model and Sandia finite
element model results may be because of some programming errors in the Mathematica
code for calculating modal frequencies and mode shapes, or because of unsuitable selection
63

of assumed modal functions. One of our future plans for this work is to replace these
assumed modal functions, which contain hyperbolic terms, with polynomial ones which are
better conditioned, and then seek better convergence for our analytical model.

64

4 Chapter 4
4.1 Conclusion and comments
As a result of this research work:
1. FEM model for troposkein-shaped Darrieus wind turbine blade free vibration has been
built for both ideal and simplified shapes. The blade vibration under flat-wise, edgewise, extensional, and torsional deformations was studied and the lowest ten modes
have been characterized for stationary and rotating blade conditions.
2. Energy equations for a troposkein blade under the same loading conditions mentioned
above were derived. We have made initial progress on the low-order model of the
curved blade. We have followed an analysis from a 1979 NASA paper by K and K, and
extended it by performing an assumed-mode discretization. However, the predicted
frequencies from the reduced-order model are not yet matching well. Possible
problems are listed below. Nonetheless, the work sets up the next student to review
and validate the analysis process, and then apply low-order modeling.
3. Possible problems in the low-order model are as follows
a) The use of clamped-clamped assumed modes. These modal functions involve
hyperbolic functions, which have extreme values when the arguments become large.
(However, the first mode won’t have this issue, and it is not producing satisfactory
estimations)
b) The use of clamped-clamped boundary conditions in the beam model may not
accurately represent the connections in the finite element analysis (However, this
should be a stiffening assumption, and lead to over-predicted frequencies)
65

c) The strain formulation followed in KK 1979 involves several rotating, deformation
dependent, coordinate transformations, and the algebra should be validated once
again. This is in progress.

4.2 Contribution
Because of the shortage in studying Darrieus wind turbine blade vibration, this work has
provided an initial contribution to research literature toward developing the study of
Darrieus wind turbine blade vibration.
The model developed here can be used for characterization of Darrieus wind turbine blade
modal vibration analytically and these results can be compared to the numerical (FEM)
results for validation.

4.3 Future work
Based on this work some follow-up studies will be useful. Some of the future work
includes.


More development in the analytical model to correct its erroneous results, which
makes it valid for analytical characterization of Darrieus wind turbine blade modal
vibration



Studying Darrieus wind turbine blade forced vibration under aerodynamic loadings



Modeling Darrieus wind turbine blade vibration considering higher order terms
using nonlinear techniques (for example perturbation methods)

66

APPENDIX

67

Table.1 blade airfoil coordinates
x
0
1
2
3
4

x/c
0
0.016393
0.032787
0.04918
0.065574

camber
line(Yc)
0
0
0
0
0

yt
0
1.640381
2.253687
2.691405
3.034709

θ
0
0
0
0
0

Xu
0
1
2
3
4

Yu
0
1.640381364
2.25368702
2.691405028
3.034709024

XL
0
1
2
3
4

YL
0
-1.640381364
-2.25368702
-2.691405028
-3.034709024

5
6
7
8
9
10

0.081967
0.098361
0.114754
0.131148
0.147541
0.163934

0
0
0
0
0
0

3.315231
3.54934
3.746877
3.914348
4.056357
4.176321

0
0
0
0
0
0

5
6
7
8
9
10

3.3152312
3.549340183
3.746876512
3.914348168
4.056356584
4.176321059

5
6
7
8
9
10

-3.3152312
-3.549340183
-3.746876512
-3.914348168
-4.056356584
-4.176321059

11
12
13
14
15
16
17
18
19
20
21
22
23
24

0.180328
0.196721
0.213115
0.229508
0.245902
0.262295
0.278689
0.295082
0.311475
0.327869
0.344262
0.360656
0.377049
0.393443

0
0
0
0
0
0
0
0
0
0
0
0
0
0

4.276882
4.360144
4.427822
4.481351
4.521944
4.550648
4.568374
4.575926
4.574017
4.563287
4.544312
4.517615
4.483671
4.442915

0
0
0
0
0
0
0
0
0
0
0
0
0
0

11
12
13
14
15
16
17
18
19
20
21
22
23
24

4.276882334
4.360143598
4.427822451
4.481351006
4.52194425
4.55064816
4.568374455
4.575926263
4.574017431
4.563287327
4.544312334
4.517614914
4.483670829
4.442914969

11
12
13
14
15
16
17
18
19
20
21
22
23
24

-4.276882334
-4.360143598
-4.427822451
-4.481351006
-4.52194425
-4.55064816
-4.568374455
-4.575926263
-4.574017431
-4.563287327
-4.544312334
-4.517614914
-4.483670829
-4.442914969

25
26

0.409836
0.42623

0
0

4.395746
4.342531

0
0

25
26

4.395746095
4.342530731

25
26

-4.395746095
-4.342530731

27
28
29
30
31
32

0.442623
0.459016
0.47541
0.491803
0.508197
0.52459

0
0
0
0
0
0

4.283606
4.219284
4.149852
4.075573
3.996694
3.913439

0
0
0
0
0
0

27
28
29
30
31
32

4.283606395
4.219284293
4.149851581
4.075573284
3.996693923
3.913438914

27
28
29
30
31
32

-4.283606395
-4.219284293
-4.149851581
-4.075573284
-3.996693923
-3.913438914

33
34
35
36

0.540984
0.557377
0.57377
0.590164

0
0
0
0

3.826016
3.734615
3.639412
3.540565

0
0
0
0

33
34
35
36

3.826015767
3.734615124
3.639411661
3.540564874

33
34
35
36

-3.826015767
-3.734615124
-3.639411661
-3.540564874

37
38

0.606557
0.622951

0
0

3.43822
3.332507

0
0

37
38

3.438219769
3.332507465

37
38

-3.438219769
-3.332507465

68

Table.1 (cont’d)
39
0.639344
40
0.655738

0
0

3.223546
3.111439

0
0

39
40

3.223545733
3.111439463

39
40

-3.223545733
-3.111439463

41
42
43
44
45
46

0.672131
0.688525
0.704918
0.721311
0.737705
0.754098

0
0
0
0
0
0

2.996281
2.878151
2.757118
2.633238
2.506559
2.377115

0
0
0
0
0
0

41
42
43
44
45
46

2.996281088
2.878150957
2.75711767
2.633238376
2.506559046
2.377114708

41
42
43
44
45
46

-2.996281088
-2.878150957
-2.75711767
-2.633238376
-2.506559046
-2.377114708

47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

0.770492
0.786885
0.803279
0.819672
0.836066
0.852459
0.868852
0.885246
0.901639
0.918033
0.934426
0.95082
0.967213
0.983607
1

0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

2.24493
2.110018
1.972382
1.832017
1.688904
1.543018
1.394322
1.242769
1.088305
0.930864
0.770372
0.606744
0.439888
0.269702
0.096075

0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

2.244929673
2.11001773
1.972382324
1.832016723
1.688904161
1.543017978
1.394321742
1.242769361
1.08830519
0.930864123
0.770371683
0.606744102
0.439888396
0.269702434
0.096075

47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

-2.244929673
-2.11001773
-1.972382324
-1.832016723
-1.688904161
-1.543017978
-1.394321742
-1.242769361
-1.08830519
-0.930864123
-0.770371683
-0.606744102
-0.439888396
-0.269702434
-0.096075

69

BIBLIOGRAPHY

70

BIBLIOGRAPHY
[1] Figure 1.1 taken from https://en.wikipedia.org/wiki/Wind_power web page.
[2] Wind Turbines: Fundamentals, Technologies, Application, Economics, by Erich Hau, 2 nd
edition, 2006, Springer, Germany.
[3] http://www.reuk.co.uk/wordpress/wind/savonius-wind-turbines.
[4] https://upload.wikimedia.org/wikipedia/commons/9/9e/Darrieus.jpg.
[5] G. J. M. Darrieus, “Turbine Having its Rotating Shaft Transverse to the Flow of the
Current, “United States Patent No. 1, 835, 018, December 8, 1931.
[6] Wind Turbine Design: With Emphasis on Darrieus Concept, By Ion Paraschivoiu,
Presses internationals Polytechnique, 2002, Canada.
[7] Blade shape for troposkien type of vertical-axis wind turbine B. F. Blackwell, G. E. Reis
Sandia National Laboratories; Albuquerque, NM 87185, April 1974.
[8] P. South and R. S. Rangi, Preliminary Tests of a High Speed - Vertical Axis Windmill
Model, National Research Council of Canada, LTR-LA-74, March 1971.
[9] P, South and R. S. Rangi, "The Performance and Economics of the VerticalAxis Wind
Turbine Developed at the National Research Council, Ottawa, Canada, " Agricultural
Engineer, February 1974, pp. 14-16. (See also American Society of Agricultural Engineers,
Paper No. PNW 73-303).
[10] “Aeroelastic equations of motion of Darrieus vertical axis wind turbine blade” by
Krishna Rao V. Kaza, and Raymond G. Kvaternik- DOE/NA SA11028-79125/NASA TM79295, December 1979.
[11] Alcoa Laboratories. Design and fabrication of a low cost Darrieus vertical axis wind
turbine system, phase 1. Technical Report, ALO-4272:1-236, June 22 1979.
[12] Bennie F. Blackwell, Robert E. Sheldahl , "Wind unnel performance data for the
Darrieus wind turnine with NACA 0012 blades" Aerothermodynamics Division, 1333 Louis
V. Feltz Mechanical Design Division, 1324 Sandia Laboratories Albuquerque, NM 87115,
May 1976.
[13] Owens, B.C., and Griffith, D.T., “Aeroelastic Stability Investigations of Large-scale
Vertical Axis Wind Turbines,” Journal of Physics Conference Series, Science of Making
Torque from Wind Conference, June 18–20, 2014, Lyngby, Denmark.

71

[14] Owens, B., Hurtado, J., Barone, M., and Paquette, J., “An Energy Preserving Time
Integration Method for Gyric Systems: Development of the Offshore Wind Energy
Simulation Toolkit” Proceedings of the European Wind Energy Association Conference &
Exhibition. Vienna, Austria, 2013.
[15] Owens, B.C., Hurtado, J.E., Paquette, J., Griffith, D.T., and Barone, M., “Aeroelastic
Modeling of Large Offshore Vertical-axis Wind Turbines: Development of the Offshore
Wind Energy Simulation Toolkit,” 54thAIAA/ASME/ASCE/AHS/ASC Structures, Structural
Dynamics, and Materials Conference, April 8–11, 2013, Boston, MA, USA, AIAA-2013-1552.
[16] Owens, B.C., Griffith, D.T., and Hurtado, J.E., “Modal Dynamics and Stability of Large
Multi-megawatt Deepwater Offshore Vertical-axis Wind Turbines: Initial Support Structure
and Rotor Design Impact Studies,” 32nd ASME Wind Energy Symposium, National Harbor,
MD, USA, January 2014.
[17] Simao Ferreira, C., Schmidt Paulson, U., Griffith, D. T., Borg, M., eds., 2016, Scientific
and Technical Challenges in Offshore Vertical Axis Wind Turbines, Euromech Colloquium
583, September 7-9, Delft, the Netherlands.
[18] Equation 2.2 taken from https://en.wikipedia.org/wiki/NACA_airfoil.
[19] www.sharcnet.ca/Software/Ansys/17.0/en-us/help/ans_elem/Hlp_E_SOLID187.
[20] G. Thomas Mase and George E. Mase, “Continuum Mechanics for Engineers”, 1999 by
CRC Press LLC, New York.
[21] Singiresu S. Rao, “Vibration of Continuous Systems”, 2007 by John Wiley & Sons, Inc,
New Jersey.

72