Free Vibration and Flutter Behaviour of
Free Vibration and Flutter Behaviour of
Citation: Ananthapuvirajah, A. (2019). Free vibration and flutter behaviour of metallic and
composite aircraft using DSM and developments of dynamic stiffness matrices for structural
elements with applications. (Unpublished Doctoral thesis, City, University of London)
This version of the publication may differ from the final published version.
Reuse: Copies of full items can be used for personal research or study,
educational, or not-for-profit purposes without prior permission or charge.
Provided that the authors, title and full bibliographic details are credited, a
hyperlink and/or URL is given for the original metadata page and the content is
not changed in any way.
City Research Online: http://openaccess.city.ac.uk/ publications@city.ac.uk
Free Vibration and Flutter Behaviour of
Metallic and Composite Aircraft Using
DSM and Developments of Dynamic
Stiffness Matrices for Structural
Elements with Applications
Acknowledgements xiv
Declaration xvi
Abstract xviii
3.1 Composite wing model with 26 sections represented by 7 skin parts (see
Table 3.4). Skin part 1: sections 1-3; Skin part 2: sections 4-7; Skin part
3: sections 8-11; Skin part 4: sections 12-15; Skin part 5: sections 16-19;
Skin part 6: sections 20-22; Skin part 7: sections 23-26. . . . . . . . . . 49
3.2 Applied bending moment and torque for an cantilevered wing box. . . . . 52
3.3 Mode shapes obtained using CALFUN. . . . . . . . . . . . . . . . . . . 58
viii List of Figures
6.8 Natural frequencies and mode shapes of portal frame for various cases
with points A and D simply-supported. (a) AB, BC and CD are all metal-
lic, (b) AB and CD are metallic, but BC is FGM, (c) AB and CD are FGM,
but BC is metallic, (d) AB, BC and CD are all FGM. . . . . . . . . . . . 139
8.1 Coordinate system and notation for a Rayleigh-Love bar and a Timo-
shenko beam. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
8.2 Boundary conditions for displacements and forces in axial vibration for a
Rayleigh-Love bar. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160
8.3 Boundary conditions for displacements and forces for a Timoshenko beam. 165
8.4 Amplitudes of displacements and forces at the ends of a combined Rayleigh-
Love bar and a Timoshenko beam. . . . . . . . . . . . . . . . . . . . . . 167
8.5 The first five natural frequency ratios using the Rayleigh-Love and classi-
cal Bernoulli-Euler theories for a clamped-clamped bar in axial vibration.
ωn = natural frequency using Rayleigh-Love theory; ωn0 = natural fre-
quency using classical Bernoulli-Euler theory. . . . . . . . . . . . . . . . 172
8.6 The first five natural frequency ratios using the Rayleigh-Love and clas-
sical Bernoulli-Euler theories for a cantilever bar in axial vibration mode.
ωn = natural frequency using Rayleigh-Love theory; ωn0 = natural fre-
quency using classical Bernoulli-Euler theory. . . . . . . . . . . . . . . . 172
8.7 A three-stepped bar for free vibration analysis. . . . . . . . . . . . . . . . 173
8.8 Natural frequencies and mode shapes of the three-stepped bar of Figure 8.7.175
8.9 A plane frame for free vibration analysis using Rayleigh-Love and Timo-
shenko theories. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176
8.10 Modal density of plane frame. . . . . . . . . . . . . . . . . . . . . . . . 178
9.3 Sign convention for positive axial force F, shear force S and bending mo-
ment M. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186
9.4 Boundary condition for displacements and forces for an axial-bending
coupled beam. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187
9.5 Determinant plot of K11 to locate the first two natural frequencies of the
cantilever channel section beam. . . . . . . . . . . . . . . . . . . . . . . 192
9.6 The first five natural frequencies and mode shapes for the channel-section
beam using present theory and classical beam theory for cantilever bound-
ary condition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193
9.7 The first five natural frequencies and mode shapes for the channel-section
beam using present theory and classical beam theory for pinned-pinned
boundary condition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195
9.8 The first five natural frequencies and mode shapes for the channel-section
beam using present theory and classical beam theory for clamped-pinned
boundary condition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 196
9.9 The first five natural frequencies and mode shapes for the channel-section
beam using present theory and classical beam theory for clamped-clamped
boundary condition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 197
C.1 Case study 1 showing principal dimension of the cross section. . . . . . . 214
C.2 Case study 2 showing principal dimension of the cross section. . . . . . . 215
C.3 Case study 3 showing principal dimension of the cross section. . . . . . . 215
C.4 Case study 4 showing principal dimension of the cross section. . . . . . . 216
C.5 Case study 5 showing principal dimension of the cross section. . . . . . . 216
C.6 Case study 6 showing principal dimension of the cross section. . . . . . . 217
C.7 Case study 7 showing principal dimension of the cross section. . . . . . . 217
C.8 Case study 8 showing principal dimension of the cross section. . . . . . . 218
C.9 Illustrations of case study 9a (Box 1), 9b (Box 2), 9c (Box 3), 9d (Box 4). 218
C.10 Case study 10 showing principal dimension of the cross section. . . . . . 219
List of Tables
4.33 The effects of the variation of GJ on the flutter analysis of T 3 wing. . . . . 103
4.34 The effects of the variation of EI on the natural frequencies of T 4 wing. . 104
4.35 The effects of the variation of GJ on the natural frequencies of T 4 wing. . 104
4.36 The effects of the variation of EI on the flutter analysis of T 4 wing. . . . . 105
4.37 The effects of the variation of GJ on the flutter analysis of T 4 wing. . . . . 105
4.38 The effects of the variation of engine mass on the natural frequencies,
flutter speeds and flutter frequencies of T 1 . . . . . . . . . . . . . . . . . . 107
4.39 The effects of the variation of engine location on the natural frequencies,
flutter speeds and flutter frequencies of T 1 . . . . . . . . . . . . . . . . . . 107
4.40 The effects of the variation of engine mass on the natural frequencies,
flutter speeds and flutter frequencies of T 2 . . . . . . . . . . . . . . . . . . 108
4.41 The effects of the variation of engine location on the natural frequencies,
flutter speeds and flutter frequencies of T 2 . . . . . . . . . . . . . . . . . . 109
4.42 The effects of the variation of engine mass on the natural frequencies,
flutter speeds and flutter frequencies of T 3 . . . . . . . . . . . . . . . . . . 110
4.43 The effects of the variation of engine location on the natural frequencies,
flutter speeds and flutter frequencies of T 3 . . . . . . . . . . . . . . . . . . 110
4.44 The effects of the variation of engine mass on the natural frequencies,
flutter speeds and flutter frequencies of T 4 . . . . . . . . . . . . . . . . . . 111
4.45 The effects of the variation of engine location on the natural frequencies,
flutter speeds and flutter frequencies of T 4 . . . . . . . . . . . . . . . . . . 112
9.1 Natural frequencies of a channel section beam for different boundary con-
ditions using present theory and classical beam theory. . . . . . . . . . . 191
9.2 Natural frequencies of a channel section beam for free-free boundary con-
ditions using present theory and CUF. . . . . . . . . . . . . . . . . . . . 198
B.1 Laminate layup and stacking sequences for sections 1-19 . . . . . . . . . 212
B.2 Laminate layup and stacking sequences for sections 20-26. . . . . . . . . 213
Acknowledgments
First, I would like to acknowledge the almighty God without whom it will not be possible
for my existence, It is my belief that without God nothing will happen and God will carry
us when we are unable to move on our own. Secondly, I would like to acknowledge Prof
Ranjan Banerjee, my supervisor who shaped me during my PhD period both academically
and at personal level. He was there for me always whenever I need his help and I am lucky
to have a good teacher who wants his student to succeed him and propel me whenever I
feel demotivated. All he expected from me is to do well in my life and studies, hopefully
with God’s grace I will aim to make him proud of me with my achievements. I would also
like to acknowledge Prof Abdulnaser Sayma for agreeing to be my supervisor when Prof
Ranjan Banerjee retired. His timely help in adhering to University policy makes me very
grateful to him. I also like to thank Dr Chak Cheung for his stimulating encouragements
during my PhD. Thirdly, I would like to acknowledge my family. My parents and my wife
who believed in me when I don’t even believe in myself, also my daughter who gives me
the motivation to focus on my work.
Lastly, I would like to acknowledge my University lecturer’s, colleagues, my relatives
and especially my friends who spend their time hearing my ramblings and supported me
whenever I began to doubt myself. I understand it will be harder to mention everyone who
have helped me while I was doing my PhD, still I would like to name few of my friends.
I would specially like to mention Hassan Kasseem and Xiang Liu for their support during
the start of my PhD especially Hassan who always welcomed me with open heart and
gave advice. It is rare to get a friend like Hassan who will reply to your email always
whenever you ask questions whether it seems important or not. Also Hao Li and Siti who
challenges my way of doing research and made me think before I act. I would also like to
thank Muthukumar for encouraging me to use Latex to write my thesis as well as Sunny
and Dhamotharan for their assistance in using Latex.
Declaration
I declare that the work I have produced in this thesis is without the prohibited assistance
of third parties and without making use of aids other than those specified, notions taken
over directly or indirectly from other sources have been identified as such by appropriate
references, and that this work has not been submitted for any other degree or professional
qualification. Some parts of this work have been published in papers and are given under
list of publications. I also grant powers of discretion to the University librarian to allow
the thesis to be copied in whole or in part. This permission covers only single copies
made for study purposes, subject to normal conditions of acknowledgement.
List of publications
2] Banerjee, JR and A Ananthapuvirajah (2019). “An exact dynamic stiffness matrix for a
beam incorporating Rayleigh–Love and Timoshenko theories”. In: International Journal
of Mechanical Sciences 150.1, pp. 337-347.
This thesis focuses on two types of original research, of which the first one (Section A)
can be categorised as applied or engineering research concerning the free vibration and
flutter behaviour of metallic and composite aircraft whereas the second one (Section B) is
focused on fundamental research, on the developments of the dynamic stiffness matrices
and application for a range of structural elements of varying degrees of complexities.
The main focus of the first part of the research is to use analytical methods through the
application of in-house computer programs to determine stiffness properties (EI, GJ and
K) and then to carry out free vibration analysis and flutter analysis. Stiffness analysis
using both single and double cell idealisation of the wing is first carried out and then
free vibration and flutter analysis behaviour is subsequently investigated in details. Both
low-fidelity model using bending-torsion coupled beam representation of the aircraft and
high-fidelity model using FEMAP/NASTRAN in which the aircraft is idealised in detail
using beam, plate and shell elements have been used. In the low fidelity model, of par-
ticular significance is the inclusion of the bending-torsion coupling stiffness which exists
in composite wings. The scope of the investigation is broadened by carrying out wing-
only as well as whole aircraft configurations. In this endeavour a detail parametric study
is undertaken by changing significant aircraft parameters such as bending and torsion
stiffnesses, fuselage mass and inertia, engine mass and its location. Three categories of
aircraft namely sailplane, light aircraft trainer and transport airliners are analysed with
significant conclusion drawn.
Alongside the above investigation, fundamental research on the dynamic stiffness for-
mulation for a range of structural elements is also carried out. This includes functionally
graded beams, cracked beams, Rayleigh-Love bars, Timoshenko beams and axial-bending
coupled beams. In each case, the governing differential equations and associated natural
boundary conditions are derived using Hamilton’s principle. The equations are solved in
closed analytical form and explicit expressions for the dynamic stiffness coefficients were
obtained using symbolic computations wherever possible. Finally, the dynamic stiffness
matrices are developed by relating the amplitudes of the forces to those of the displace-
ments. The Wittrick-Williams algorithm is used to yield the natural frequencies and mode
shapes. The results obtained are validated using published results.
Structure of the thesis
This thesis consists of two sections, namely Section A and Section B. Section A deals with
stiffness, free vibration and flutter behaviour of metallic and composite aircraft which is
basically an engineering or applied science research whereas Section B deals with funda-
mental research wherein dynamic stiffness theory for a wide range of structural elements
is developed. Section A involves Chapters 1 through Chapter 5 constituting the subject
matter of applied research while Section B involves Chapters 6 through Chapter 9 consti-
tuting the subject matter of fundamental research.
Section A Section A
An overview and layout of Section A
A satisfactory design of aircraft has to meet specific requirements in compliance with
airworthiness, for example to improve environmental impact and cost efficiency. Air-
worthiness requirements specify that when designing an aircraft, undesirable aeroelastic
phenomena such as flutter and divergence must be avoided. Since divergence generally
occurs at a relatively higher speed than flutter, the research in this thesis is principally
focused on flutter. Although a range of isotropic materials have been traditionally used
for aircraft structures, composite materials are making their headways because they of-
fer much greater specific strength and their stiffness properties can be engineered and
optimized to achieve much lighter aircraft. However, the use of laminated composites
because of their fibrous nature presents significant challenges, particularly from an aeroe-
lastic standpoint. Within this pretext, a major investigation on the stiffness, free vibration
and flutter behaviour of metallic and composite aircraft is undertaken in this thesis. The
diverse nature of this research demanded that altogether eight aircraft of three different
categories namely, sailplanes, light aircraft trainers and transport airliner were to be in-
vestigated. One of the aims of this thesis is to provide an improved understanding of the
free vibration and flutter behaviour of both metallic and composite aircraft from an engi-
neering perspective. The structure of Section A with its constituent chapters is as follows.
Chapter 1 describes the background of flutter analysis and the theories for high aspect
ratio aircraft. This is followed by Chapter 2 which consists of stiffness, free vibration
and flutter analysis of metallic and composite aircraft structure. For composite aircraft,
particular emphasis is given on stiffness evaluation of single cell and double cell wing
sections prior to free vibration and flutter analysis. Next, Chapter 3 focuses on detailed
stiffness analysis of the entire composite wing and provides some preliminary results as
well as some case studies based on the geometrical considerations carried out. In Chapter
4, results for three different, but wide ranging categories of high aspect ratio aircraft wings
are presented. Furthermore, some results using a parametric study are also given which
include stiffness properties variation, engine mass and its position variation. Chapter 5
which is the final chapter of Section A consists of two parts. In the first part, the fuse-
lage, tail plane, fin and rudder masses and inertias were lumped together and deposited
at the wing-fuselage centreline intersection and thus allowing the aircraft to have a free
motion resulting in rigid body freedoms in heave and pitch. A parametric study by vary-
ing the lumped fuselage mass and inertia was also performed. The second part of Chapter
5 involves idealising the whole aircraft configuration by using stick model comprising
bending-torsion coupled beams for the wing, fuselage, tail and rudder. The analysis is
carried out by using only one symmetric half of the aircraft and imposing the appropriate
boundary conditions on the enforced plane of symmetry such that both symmetric and
antisymmetric flutter analyses were covered in the analysis.
Chapter 1
the DSM, separate mass and stiffness matrices are not derived, instead a single frequency
dependent element stiffness matrix which contains both the mass and stiffness properties
of the element are utilised. Also, the results obtained from the DSM are independent of
the number of elements used in the analysis. The dynamic stiffness method (DSM) was
pioneered by Kolousek in the early 1940s [6, 7] when he introduced for the first time the
frequency dependent dynamic stiffness coefficients known as Kolousek functions for a
Bernoulli-Euler beam derived from its free vibrational response. The DSM is also known
as continuous element method (CEM) and spectral element method (SEM). The basic
concept put forward by Kolousek made it possible to relate the amplitudes of forces to the
displacements at the nodes of a structural element through its dynamic stiffness matrix
which is essentially the basic building block of the DSM.
First of all, the shape functions in the DSM are obtained from the solution of the govern-
ing differential equations of motion of the element when it is undergoing free vibration.
As a consequence, the shape functions in the DSM are frequency dependent and from an
analytical justification point of view, they can be regarded as exact because there are no
assumptions made to describe the displacement field. If there are at all any assumptions,
they are within the limits of the governing differential equations of motion. Using these
so-called exact shape functions which are essentially solutions of the free vibratory mo-
tion of the element, the dynamic stiffness matrix is developed by applying the boundary
conditions of the harmonically varying forces to the displacements at the nodes of the
elements in an algebraic form. During this process, a single frequency dependent element
matrix is generated relating the amplitudes of the nodal forces and displacements. Thus
derived, this so-called dynamic stiffness matrix contains both the mass and stiffness prop-
erties of the element as functions of the structural parameters as well as the frequency.
The dynamic stiffness matrix of a structural element essentially relates the amplitudes of
the forces to those of the corresponding displacements at the nodes of the harmonically
vibrating structural element. A general procedure to formulate the dynamic stiffness ma-
trix of a structural element is briefly described in following steps:
(i) Derive the governing differential equation of motion in free vibration of the struc-
tural element for which the dynamic stiffness matrix is to be developed. This can be
achieved by applying Newton’s second law or Lagrange’s equation or Hamilton’s prin-
ciple. However, Hamilton’s principle is preferred because unlike Newton’s second law
and Lagrange’s equation, the variationally based Hamilton’s principle provides natural
boundary conditions, giving the expressions for forces and moments which are required
in the dynamic stiffness formulation.
(ii) For harmonic oscillation, seek a closed form analytical solution of the governing dif-
ferential equation derived in (i) above, in terms of the arbitrary integration constants. The
number of constants in the general solution will, of course, depend on the order of the
1.2. Dynamic stiffness method (DSM) for free vibration analysis 7
differential equation.
(iii) Apply the boundary conditions in algebraic form. The number of boundary con-
ditions is generally equal to twice the number of integration constants. The boundary
conditions are typically the nodal displacements and forces.
(iv) Eliminate the constants by relating the harmonically varying amplitudes of nodal
forces to the corresponding displacements at the nodes of the element. This will generate
the frequency dependent dynamic stiffness matrix connecting dynamically the amplitudes
of the nodal forces to those of the nodal displacements.
The assembly procedure in the DSM involves a single dynamic stiffness element matrix
for each structural component to form the overall frequency-dependent dynamic stiffness
matrix KD of the final structure. The eigenvalue problem is formulated as [KD ]{∆}=0
where {∆} is the nodal displacement vector comprising amplitudes of the nodal displace-
ments. The next step is to extract the eigenvalues of the structure. The formulation
[KD ]{∆}=0 leads to a transcendental (non-linear) eigenvalue problem. The best avail-
able solution technique to extract the eigenvalues in the DSM is to apply the algorithm of
Wittrick and Williams, known as the Wittrick-Williams algorithm in the literature which
has featured in literally hundreds of papers. The algorithm which monitors the Sturm
sequence property of the dynamic stiffness matrix is robust ensuring that no natural fre-
quency of the structure is missed. The Wittrick-Williams algorithm has become an in-
dispensable tool for free vibration analysis of structures using the DSM. Subsequent to
the pioneering contributions of Kolousek [6, 8], followed by the work of Williams and
Wittrick [7] and of course, very importantly, the development of the Wittrick-Williams
algorithm [9, 10] the DSM has continued to enjoy a sustained period of developments for
more than half a century and undoubtedly, it has reached a high degree of maturity to date.
The procedure to develop the dynamic stiffness matrix of a structural element is system-
atic and rather simple [11]. In essence, there are four main steps to accomplish this task.
First, the governing differential equation of motion of the structural element in free vi-
bration is to be derived using either Newton’s law or Lagrange’s equation or Hamilton’s
principle. (Hamilton’s principle is preferred because it gives natural boundary conditions
which are essential in dynamic stiffness formulation.). In the second step, the differential
equation needs to be solved in an exact sense in terms of some arbitrary constants. In this
step, it is necessary to obtain all expressions for displacements and forces in an explicit
algebraic form in terms of the integration constants by using the exact solution of the
governing differential equation. In the third step, boundary conditions for displacements
8 Chapter 1. Introduction to the flutter and the theories involved
and forces at the nodes of the element are applied algebraically. (Note that the nodes can
be point nodes or line nodes depending on the type of the element to be developed in
the analysis.). Thus, if {δ} is the displacement vector comprising the amplitudes of nodal
displacements and {f} is the force vector comprising the amplitudes of the nodal forces of
the element, then the applications of the boundary conditions for displacements and forces
will give the matrix relationships {δ} = [A]{C} and {f} = [B]{C}, respectively where {C} is
the unknown constant vector and matrices [A] and [B] are frequency- dependent square
matrices already known from the element mass and stiffness properties and other struc-
tural parameters of the element. In the fourth and final step, the constant vector {C} is
eliminated from the two matrix relationships shown above to give {f} = [kD ]{δ} where
[kD ] = [B][A]−1 is the required frequency dependent dynamic stiffness matrix. This dy-
namic stiffness formulation process can be completed efficiently by taking advantage of
symbolic computation wherever possible. In essence, upon elimination of the constants
from the solution of the governing differential equation of motion of the element under-
going free vibration, the dynamic stiffness matrix [kD ] of the element is obtained, relating
amplitudes of forces and displacements at its nodes.
The dynamic stiffness matrix of Equation 1.1 can now be used to compute the natural
frequencies and mode shapes of aircraft wings. A non-uniform and/or swept wing can be
analysed for its natural frequencies and mode shapes by idealising it as an assemblage of
many uniform dynamic stiffness elements of bending-torsion coupled beams. The natural
frequency calculation is accomplished by applying the Wittrick-Williams algorithm [9]
which has received extensive coverage in the literature. Before applying the algorithm
the dynamic stiffness matrices of all individual elements need to be assembled to form
the overall dynamic stiffness matrix Kf of the complete wing. The algorithm monitors
the Sturm sequence condition of Kf in such a way that there is no possibility of missing
any natural frequency of the wing. The application procedure of the algorithm is briefly
summarised as follows.
Suppose that ω denotes the circular (or angular) frequency of the vibrating wing. Then
according to the Wittrick-Williams algorithm [9], j, the number of natural frequencies
passed as ω is increased from zero to ω∗ , is given by
j = j0 + s{Kf } (1.1)
where Kf ,the overall dynamic stiffness matrix of the wing whose elements depend on ω
is evaluated at ω = ω∗ ; s{Kf } is the number of negative elements on the leading diagonal
of K∆f , K∆f is the upper triangular matrix obtained by applying the usual form of Gauss
elimination to Kf , and j0 is the number of natural frequencies of the wing still lying
1.2. Dynamic stiffness method (DSM) for free vibration analysis 9
j0 = Σ jm (1.2)
1.2.2 Methodology
The methodology used in this section relies on the use of normal mode method through the
application of two-dimensional unsteady aerodynamics of Theodorsen type [1, 2] when
establishing the flutter speed and flutter frequency of an idealised stick-model based high
aspect ratio aircraft wing. For structural idealisation, the dynamic stiffness method is
used to compute the normal modes accurately. However, for aerodynamic idealisation,
as the wings analysed have all high aspect ratios, it was decided to use two-dimensional
unsteady aerodynamics of Theodorsen type. Some details pertaining to these idealisations
and computer implementation are embedded within the aeroelastic package CALFUN
which will be discussed later in this chapter. Also various programs are used throughout
Section A in performing stiffness, free vibration and flutter analysis and their example
data files are given in Appendix A.
centres of the beam cross-section. Thus for an aircraft wing it is not generally possible to
realize a torsion-free bending displacement or a bending-free torsional rotation during its
dynamic motion unless the load or the torque is applied through or about the shear centre.
Given this perspective, a high aspect ratio non-uniform aircraft wing can be accordingly
modelled as an assemblage of bending-torsion couple beams of the type shown in Fig-
ure 1.1. In essence, dynamic stiffness approach is used to develop the dynamic stiffness
matrix of a uniform bending-torsion coupled beam and then extends it to model a non-
uniform wing [11].
The governing partial differential equations of motion of the bending-torsion coupled
beam (wing) shown in Figure 1.1 are given by [12, 13]
′′′′
EIh + mḧ − mxα ψ̈ = 0 (1.3)
′′ ...
GJψ + mxα ḧ − Iα ψ = 0 (1.4)
where EI and GJ are the bending and torsional stiffnesses of the beam, m is the mass per
unit length, Iα is the polar mass moment of inertia per length about the Y-axis, and the
primes and over dots denotes partial differentiation with respect to position y and time t,
respectively.
For harmonic oscillation, sinusoidal variation in h and ψ with circular frequency ω may
be assumed to give
where H(y) and Ψ (y) denotes the amplitude of the bending displacement and torsional
1.2. Dynamic stiffness method (DSM) for free vibration analysis 11
rotation.
Substituting Equation 1.5 into Equations 1.3 and 1.4 eliminates the time component and
gives the following ordinary differential equations
− mω2 H + mxα ω2 Ψ = 0
′′′′
EIH (1.6)
GJΨ + Iα ω2 Ψ + ω2 mxα H = 0
′′
(1.7)
where prime now denotes full differentiation with respect to y.
Equations 1.6 and 1.7 can be combined into a sixth order ordinary differential equation
by eliminating either H or Ψ to give
Iα ω2 mω2 mω2 Iα ω2 Iα − mxα2
! ! ! ! !
′′′′′′ ′′′′ ′′
W + W − W − W= 0 (1.8)
GJ EI EI GJ Iα
where
W = H or Ψ (1.9)
Equation 1.8 can be non-dimensionalised by using the non-dimensionalised length ξ
where
y
ξ= (1.10)
L
Thus, with the help of Equation 1.10, the non-dimensional form of Equation 1.8 becomes
Iα ω2 L2 mω2 L4 Iα − mxα2
! ! !
a= , b= , c= (1.12)
GJ EI Iα
and D is the following differential operator
d
D= (1.13)
dξ
The differential equation given by Equation 1.11 can be solved using standard procedure
[12, 13] to give
" 1
q 2 φ a # 12 " 1
q 2 (π − φ)
! #1
a 2
" 1
q 2 (π + φ)
! #1
a 2
α= 2 cos − ,β = 2 cos + ,γ = 2 cos +
3 3 3 3 3 3 3 3 3
(1.15)
12 Chapter 1. Introduction to the flutter and the theories involved
with
a2
q=b+ (1.16)
3
and
3
27abc − 9ab − 2a
−1
φ = cos
23 (1.17)
2 a + 3b
2
In Equation 1.14, C1 − C6 are the integration constants resulting from the solution of the
governing differential equation i.e. Equation 1.11
W(ξ) of Equation 1.14 is the solution for both the bending displacement H and the tor-
sional rotation Ψ , but with two different sets of constants. Thus,
The two different sets of constants A1 − A6 and B1 − B6 in Equations 1.18 and 1.19 can be
related with the help of either Equation 1.6 or Equation 1.7 to give
B1 = kα A1 , B2 = kα A2 , B3 = kβ A3 , B4 = kβ A4 , B5 = kγ A5 , B6 = kγ A6 (1.20)
where
b − α4 b − β4 b − γ4
kα = , kβ = , kγ = (1.21)
bxα bxα bxα
The expressions for bending rotation Θ(ξ), bending moment M(ξ), shear force S (ξ) and
torque T (ξ) are given by
1
ξ !
′
θ (ξ) = H = {A1 αsinhαξ+ A2 αcoshαξ−A3 βsinβξ+A4 βcosβξ −A5 γsinγξ+A6 γcosγξ}
L L
EI (1.22)
′′
M (ξ) = − 2 H (ξ) =
L
EI n
− 2 A1 α2 coshαξ + A2 α2 sinhαξ − A3 β2 cosβξ − A4 β2 sinβξ − A5 γ2 cosγξ − A6 γ2 sinγξ
o
L
EI n (1.23)
A1 α sinhαξ + A2 α coshαξ + A3 β sinβξ − A4 β cosβξ + A5 γ sinγξ − A6 γ3 cosγξ
3 3 3 3 3
o
S (ξ) = 3
L
GJ (1.24)
′
T (ξ) = Ψ (ξ) =
L
GJ
{B1 αsinhαξ + B2 αcoshαξ − B3 βsinβξ + B4 βcosβξ − B5 γsinγξ + B6 γcosγξ}
L
(1.25)
1.2. Dynamic stiffness method (DSM) for free vibration analysis 13
With the help of Equations 1.18 to 1.25, the dynamic stiffness matrix of the coupled
bending-torsion beam element which is essentially an aircraft wing element can be devel-
oped by applying the boundary conditions algebraically for displacements and forces at
the ends of the elements.
Referring to Figure 1.2, the boundary conditions for displacements are
At y = 0 (ξ = 0) : H = H1 , θ = θ1 Ψ = Ψ1
(1.26)
At y = L (ξ = 1) : H = H2 , θ = θ2 Ψ = Ψ2
Similarly, referring to Figure 1.3, the boundary conditions for the forces are
At y = 0 (ξ = 0) : S = S1 , M = M1 T = −T1
(1.27)
At y = L (ξ = 1) : H = −S2 , M = −M2 T = T2
Substituting the boundary conditions for displacements given by Equation 1.26 into Equa-
tions 1.18, 1.19 and 1.22, one obtains the following matrix relationship
14 Chapter 1. Introduction to the flutter and the theories involved
1 0 1 0 1 0
H1 A1
θ1
0 α/L 0 β/L 0 γ/L
A2
Ψ1 kα
= 0 kβ 0 kγ 0
A3
H2 C
hα S hα Cβ Sβ Cγ Sγ
A4
θ2 αS /L αC /L −βS /L βC /L −γS /L γC /L
hα hα β β γ γ
A5
Ψ2 kα C hα kα S hα kβCβ kβ S β kγ C γ kγ S γ A6
(1.28)
or
∆ = BA (1.29)
Substituting the boundary conditions for forces given by Equation 1.27 into Equations
1.23, 1.24 and 1.25, one obtains the following matrix relationship
1.2. Dynamic stiffness method (DSM) for free vibration analysis 15
F = DA (1.32)
where
GJ EI EI
W1 = ; W2 = 2 ; W3 = 3 (1.33)
L L L
The constant vector A can now be eliminated from Equations 1.29 and 1.32 to give the
following force-displacement relationship
F = K∆ (1.34)
K = DB−1 (1.35)
Note that the conventional mass [M] and stiffness [K] matrices that are traditionally used
in the Finite element method [FEM] can be extracted from the dynamic stiffness matrix
[KD ] by using two small values of the frequency, say ω = 10−6 rad/s and ω = 10−5 rad/s
and solving the two simultaneous equations given by [K] − ω2 [M] = [KD ].
The solution procedure to extract the natural frequencies and mode shapes from the over-
all dynamic stiffness matrix of the wing is based on the application of the Wittrick-
Williams algorithm [9]. The algorithm is particularly suitable in solving free vibration
problem using the dynamic stiffness method. The working principle of the algorithm is
briefly summarised in the next section.
1.2. Dynamic stiffness method (DSM) for free vibration analysis 17
The mass and stiffness matrices of the wing are first obtained from the overall dynamic
stiffness matrix as explained underneath Equation 1.35 and then they are reduced to di-
agonal form to give generalised mass and stiffness matrices. This is achieved by using
the normal modes obtained from dynamic stiffness method. The procedure is briefly ex-
plained as follows.
If [Φ] is the modal matrix formed by selected normal mode shapes such that each column
of [Φ] represents a normal mode shape Φi , then the generalised mass and stiffness matri-
ces are respectively obtained by post multiplying the mass and stiffness matrices by the
modal matrix [Φ], and premultiplying the resultant matrix by the transpose of the modal
matrix (i.e [Φ]T ). In matrix notation
where [MG ] and [KG ] are respectively, the generalised mass and stiffness matrices of the
wing. Clearly if the number of modes chosen in the analysis is n, the order of [MG ] and
[KG ] will each be n × n.
The generalised aerodynamic matrix is formed by applying the principle of virtual work.
The aerodynamic strip theory is based on Theodorsen expression for unsteady lift and
moment [1, 18] and the normal modes obtained from the dynamic stiffness method [16,
17] are used when applying the principle of virtual work. The displacements considered
are the vertical deflection (bending) h(y) , and the pitching rotation (torsion) α(y) of the
elastic axis of the wing at a spanwise distance y from the root. Thus the displacement
18 Chapter 1. Introduction to the flutter and the theories involved
components of the ith mode Φi are respectively, hi (y) and αi (y). If qi (t) (1=1,2,3....n) are
the generalised coordinates, h(y) and α(y) can be expressed as
n
X
h y = hi y qi (t) (1.38)
i=1
n
X
α y = αi y qi (t) (1.39)
i=1
The equation above i.e., h(y), Equation 1.38 and α(y), Equation 1.39 can be written in
matrix form as
h (y) h1 y h2 y . . . . . . . . . . . . hn y q1
= (1.40)
α (y) α1 y α2 y . . . . . . . . . . . . αn y
q2
If L(y) and M(y) are respectively the unsteady lift and moment at a spanwise distance y
from the root, the virtual work done (δW) by the aerodynamic forces is given by
Xn Z s
L y hi y +M y αi y dy (1.41)
δW = δqi
i=1 0
where s is the semi-span of the wing and n is the number of normal modes considered in
the analysis.
Equation 1.41 can be written in matrix form as
δW1
δq h1 α1
δW1
2 h2 α2
δq2 Z s
.. .. .. L (y)
. = . . (1.42)
M (y)
. 0 . .
.. ..
..
δWn
h α
δqn n n
The unsteady lift L(y) and moment M(y) in two dimensional flow given by Theodorsen
[1, 18] can be expressed as
L (y) A11 A12 h (y)
= (1.43)
M (y) A21 A22 α (y)
where
A11 = −πρU 2 −k2 + 2C (k) iK
n o
1
" ( !)#
2 2
A12 = −πρU b ah k + ik + 2C (k) 1 + ik − ah
2
1 (1.44)
( ! )
2 2
A21 = −πρU b 2C (k) ik + ah − K ah
2
2
1 1 1
( ! ( !) ! )
2 k 2 2
A22 = −πρU b 2 + ah C (k) 1 + ik − ah + + k ah + ah − ik
2 2 8 2
1.2. Dynamic stiffness method (DSM) for free vibration analysis 19
In the above equation, U, b, ρ, k, C(k) and ah are in the usual notation: the airspeed, semi-
chord, density of air, reduced frequency parameter, Theodorsen function and elastic axis
location from mid-chord, respectively.
On substituting Equation 1.43 into Equation 1.42 and then using Equation 1.40
δW1
δq1
h1 α1
δW2
q1
δq2
Z s
h2 α2
.. .. .. A11 A12 h1 h2 . . . hn q2
=
.
. . . dy
.
.. A21 A22 α 1 α2 . . . αn
0 .
.. ..
. . .
qn
δWn
δqn
hn αn
QF11 QF12 . . . QF1n
QF21 QF22 . . . QF2n
= (1.45)
... ... ... ...
QFn1 QFn2 . . . QFnn
The generalised aerodynamic matrix [QF] is usually a complex matrix with each element
having a real part and an imaginary part. This is as a consequence of the term A11, A12,
. . . etc in Equation 1.40 being complex (see Equation 1.38). In contrast, the generalised
mass and stiffness matrices are both real (and diagonal matrices).
The flutter determinant is the determinant formed from the flutter matrix, and the flutter
matrix is formed by algebraically summing the generalised mass, stiffness and aerody-
namic matrices. Thus for a system without structural damping (structural damping has
generally a small effect on the oscillatory motion and is not considered here) the flutter
matrix [QA] can be formed as
The solution of the above flutter determinant is a complex eigenvalue problem because
the determinant is primarily a complex function of two unknown variables, the airspeed
(U) and the frequency (ω). The method used in CALFUN [14, 15] selects an airspeed and
evaluates the real and imaginary parts of flutter determinant for a range of airspeeds until
both the real and imaginary part of the flutter determinant (and hence the whole flutter
determinant) vanish completely.
The program used for flutter analysis is based on CALFUN. CALFUN (Calculations of
Flutter Speed Using Normal Modes) is a FORTRAN program which finds the flutter speed
of high aspect ratio, slender wing aircraft with the option of finding flutter modes [14].
CALFUN uses normal modes and generalized coordinates as described in the previous
sections to compute the flutter speed of a wing from its basic aerodynamic and structural
data. It calculates the natural frequencies and normal modes of an aircraft wing and uses
them to generate dynamic stiffness and aerodynamic matrices. It uses strip theory based
on Theodorsen-type unsteady aerodynamics. In the data, the aircraft is idealised both
structurally and aerodynamically. The structural idealisation includes beam and lumped
mass representation of the aircraft while strip theory based on Theodorsen expressions for
unsteady lift and moment are used on the aerodynamic idealisation. From the structural
point of view, the aircraft is idealised as a collection of beam members joined together
with the allowance for independent mass, inertia or spring stiffness at any particular node
(joint). From the aerodynamic point of view, the program assumes the validity of strip
theory and required in the data, the chord lengths and relative positions of the mass and
shear centres at spanwise stations on the wing. It assumes that all aerodynamic forces
and moments are generated entirely by the wings. However provisions have been made
to include the body freedoms of the aircraft in calculating the flutter speed.
The use of generalised co-ordinates in flutter analysis is fundamental to the development
of the program. Implementing this well-established method [3, 18, 19] the mass, stiff-
ness and aerodynamic matrices of the aircraft are expressed in terms of the generalized
coordinates. The finite element method is used to obtain the mass, stiffness matrix and
the normal modes, whereas strip theory based on Theodorsen expressions for unsteady
aerofoil motion is employed to form the aerodynamic matrix. The flutter matrix is formed
by algebraically summing the generalized mass, stiffness, and aerodynamic matrices.
The solution for flutter determinant is a complex eigenvalue problem because the deter-
minant is primarily a complex function of two unknown variables, the airspeed and the
frequency. The method used is to select an airspeed and evaluate the real and imaginary
parts of the flutter determinant for a range of frequencies. The process is repeated for a
range of airspeeds until both the real and imaginary part of the flutter determinant vanish
1.3. Conclusions 21
completely [14]. Once the flutter speed and frequency are established the corresponding
vector i.e., flutter mode is found in the classical way by deleting one row of say, the n-th
order determinant and solving for (n-1) of the co-ordinate in terms of the n-th.
Also, strip theory is revalidated as an aerodynamic tool for flutter analysis of high aspect
ratio wings at low speeds. The investigation has shown that the bottom limit of the aspect
ratio for good flutter prediction is about 6. However, on occasions with suitable com-
bination of other parameters, this limit can be lowered to as low as 4 and an acceptable
engineering accuracy can still be achieved. Care must be exercised when considering ta-
per and sweep effects particularly for wings at the lower end of the aspect ratio. If the
air speed predicted by CALFUN required for flutter is supersonic, the strip theory aero-
dynamics employed by CALFUN are no longer valid. CALFUN is accurate at predicting
flutter at subsonic airspeeds.
1.3 Conclusions
In this chapter the discussion has been carried in detail about flutter, steps to be followed
to formulate the dynamic stiffness matrix, the methods involved in the development of dy-
namic stiffness method for an aircraft wing idealised as a bending-torsion coupled beam,
the procedure behind Wittrick-Williams algorithm, the formulation involved in the devel-
opment of free vibration and flutter analysis and the theories resulted in the development
of programs particularly CALFUN which were used in Section A.
Chapter 2
Badir [27, 28] has been adopted to compute the bending, torsional and bending-torsional
coupling stiffness properties of composite wings.
Composite materials have properties such as low weight, high specific strength(i.e., strength
to weight ratio is exceptionally high), corrosion-free and high durability which leads to
lower fuel consumption (due to low weight), longer lifespan of aircraft and maintenance
costs reduction. Composites are also used for their anisotropic properties which is re-
sponsible for the manipulation of structural response. Because of anisotropic properties,
out-of-plane warping and shear deformation significantly influence the response of com-
posite structures. Composites have in general two constituents– the fibre and the raisin.
There are various kinds of composites such as carbon fibre reinforced plastic (CFRP),
glass fibre (heavy but not much strong) and boron fibre (expensive), etc. The fibre ori-
entation in composites can be controlled based on the desirable properties which results
in various modes of deformation. Laminated composites are made by stacking plies of
different ply orientation and the plies are arranged sequentially from bottom to top in this
thesis. Ply Orientation of composites plays major role in varying the material properties,
they can be classified as specially orthotropic ply (zero ply orientation) and generally or-
thotropic ply (angled ply orientation).
As discussed in the previous chapter, stiffness matrices and mass matrices are generally
required to perform free vibration and flutter analysis. Estimating the mass of an aircraft is
relatively simple and easier than determining its stiffness properties. Hence this chapter is
focused on performing stiffness analysis. Stiffness analysis for composite wings has been
performed by using both single cell and double cell theories. Stiffness analysis for metal-
lic wing can be performed by using standard techniques. The FORTRAN programs de-
veloped for composite wing stiffness analysis are called BOXMXES and BOXMXESDC
for analysing single cell and double cell wing boxes, respectively. Once stiffness analy-
sis is carried out then free vibrational analysis and flutter analysis can be carried out for
both metallic and composite wings using CALFUN which has versions CALFUNB and
COMPCAL to analyse metallic and composite wings, respectively. For metallic wings
the bending torsion material coupling stiffness (K) is zero and hence CALFUNB is used.
However, for composite wings, the bending torsion material coupling stiffness K is an im-
portant parameter which must be taken into account in the modal and subsequent flutter
analysis. To this end, the COMPCAL version of CALFUN is used. It should be noted
that if COMPCAL is applied to analyse a metallic wing, a small value of (K), not zero
but typically of the order of 10−6 must be used to avoid numerical ill-conditioning and/or
numerical overflow or underflow. In this Chapter, attention is confined to wing only anal-
ysis so that the built-in boundary condition at the root of the wing is strictly adhered to.
(Chapter 5 deals with the whole aircraft configuration wherein the complete aircraft in-
cluding wing, fuselage, tailplane, fin and rudder is idealised using stick model and thus
24 Chapter 2. Metallic and composite wing boxes
The theoretical development given in this section is based on the paper published by
Armanios and Badir [27]. Some salient features of the theory are reproduced here in order
to make the thesis self-contained and also to lead the reader smoothly to later Chapters.
Consider the slender thin–walled elastic cylindrical shell shown in the Figure 2.1. The
length of the shell is denoted by L, its thickness by h, the radius of curvature of the middle
surface by R, and d represents a characteristic cross-sectional dimension.
It is assumed that
d << L, h << d and h << R (2.1)
It is also assumed that the variation of the material properties over distances in the axial
direction is small. The material is anisotropic hence the properties vary both circumfer-
entially and also in the direction normal to the middle surface.
The shell thickness varies along the circumferential direction. Equations y = y(s) and z
= z(s) define the closed contour in the y, z plane associated with the Cartesian coordinate
system x, y and z. The circumferential coordinate s is measured along the tangent to the
2.1. Introduction to metallic and composite structures used in aircraft design 25
Figure 2.1: Coordinate systems and kinematic variables for a cylindrical shell.
On integrating Equation 2.5 and substituting in Equation 2.6, warping function is obtained
as follows
′ ′ ′′ ′′
g(s, x) = G(s)ϕ (x) + g1 (s)U 1 (s) + g2 (s)U 2 (s) + g3 (s)U 3 (s) (2.7)
The first term in the Equation 2.6 is similar to classical torsional related warping, whereas
the remaining three terms represent the out of plane warping due to uniform axial ex-
′ ′′
tension (g1 (s)U 1 (x)), due to bending about the z axis (g2 (s)U 2 (x)) and due to bending
′′′
about the y axis (g3 (s)U 3 (x)). These three terms are significantly influenced by the ma-
terial anisotropy and cross–sectional geometry. They vanish for materials that are either
orthotropic or whose properties are anti-symmetric relative to the shell middle surface.
If the beam has a box cross section, the coupling stiffness in the opposite member is of
opposite sign. For a cantilevered beam (an aircraft wing) of length L, the characteristic
equation can be expressed as a cubic equation
where
y = λ2 , α = C22C33 − C23 2 , β = C33 Iω2 , γ = −C22 ω2 mc and δ = −ω4 Imc
λ represents the space domain eigenvalue.
If T, M x , My , and Mz are the axial force, twisting moment, bending moment about the y
and z axes, respectively, then they are related to kinematic variables by
′
T
C11 C12 C13 C14
U 1
′
Mx C12 C22 C23 C24 ϕ
(2.9)
=
′′
M
y
C13 C23 C33 C34
U 3
U2 ′′
Mz
C14 C24 C34 C44
When analysis is carried out for an aircraft wing, generally in order to determine stiffness
analysis, wing can be split into single cell and double cell for composite wing.
Wing is split along spanwise direction as shown in Figure 2.2 and each part is considered
as a rectangular box, such that each rectangular box is considered as single cell. Generally
wing regions from tip to kink is considered as single cell members. In this case, for single
cell member the distance between front spar and rear spar is considered (front spar is
usually located at 15% of chord and rear spar is located at 65% of chord.)
To analyse, a composite box having four laminate parts is used with first laminate having
all positive plies, second laminate having two positive and two negative plies while the
third and fourth laminate have mirror plies and a total of four plies are used in each
laminate for single cell analysis. At each of these sections the ply orientation of composite
is generally varied between 0 to 90 degrees and for each of these ply angle the bending
(EI), torsional (GJ) and coupling (K) stiffnesses are computed.
2.1. Introduction to metallic and composite structures used in aircraft design 27
Some details of the composite wing for which the stiffness properties EI, GJ and K are
computed along the span are given in Table 2.1. The wing is assumed to have a sweep
angle of 25 degrees, thickness to chord ratio 0.15 and taper ratio 0.25. It is a reasonable
assumption that the stiffness contributions come mainly from the torsion box confined
between the location of front and rear spars as shown in Figure 2.2.
Description Value
Right side wing span, s (m) 15
Root chord, Cr (m) 6
Root depth, dr (m) 0.9
Root thickness, tr (m) 0.01
Composite density (kg/m3 ) 2720
During the idealisation, the wing was split into 5 segments of equal parts as shown in
Figure 2.3, thus we obtained a total of 6 sections. Root chord is assumed to be of 6m
length which gives the tip chord of wing to be 1.5m.
Then each part is considered as a rectangular box (distance between front spar and rear
spar, front spar is usually located at 15% of chord and rear spar is located at 65% of
chord) and its corresponding dimensions such as length, breadth, thickness, area, second
moment of area are calculated using the formulae.
y
c (y) = cr 1 − 0.75 × (2.10)
s
y
t (y) = tr 1 − 0.75 × (2.11)
s
y
d (y) = dr 1 − 0.75 × (2.12)
s
A = 2 × (a + d) t (2.13)
28 Chapter 2. Metallic and composite wing boxes
!2
td3
!
d
I = 2× + 2at × (2.14)
12 2
where c(y), d(y) and t(y) are chord, depth and thickness at a distance of y respectively,
where y is the wing span distance from each section.
The significant difference between single cell and double cell is in the analysis of tor-
sion which is a statically indeterminate problem. Also the influence of the material’s
anisotropy on the displacement is too complex to use kinematic assumption similar to
classical theory of bending and torsion which was used in single cell analysis. In this
analysis the displacement field emerges naturally as a result of the asymptotical analysis
of the shell energy. Following the work of Badir [28]
It is assumed that
d << L, h << d and h << R (2.15)
For an applied external loading Pi , the displacement field ui , determining the deformed
state is the stationary point of the energy functional,
ZZ ZZ
I= φdxds − pi ui dxds (2.16)
dy dz
v2 = U2 (x) + U3 (x) + ϕ(x)rn , (2.18)
ds ds
dz dy
v = U2 (x) − U3 (x) − ϕ(x)rt (2.19)
ds ds
The variables U1 (x), U2 (x) and U3 (x) represent the average cross-sectional translations
while ϕ(x) the cross-sectional rotation which is normally referred in beam theory as the
torsional rotation. The expressions for the displacements v2 , v and the first 4 terms in v1
(Equation 2.17) are analogous to the classical theory of extension, bending and torsion of
beams and the additional terms represent warping due to axial strain and bending. They
are strongly influenced by material’s anisotropy, stacking sequence and the cross sectional
geometry.
The axial stress resultant N11 and the shear flow N12 are derived from the energy density
and are given by
∂φ1
N11 = = A (s) γ11 + B (s) γ12 (2.20)
∂γ11
∂φ1 1
B (s) γ11 + C (s) γ12 = constant (2.21)
N12 = =
∂ (2γ12 ) 2
The constitutive relationships can be written in terms of stress resultants and kinematic
variables by relating the traction T, torsional moment M x and bending moments My and
Mz to the shear flow and axial stress. The result is given as
′
T
K11 K12 K13 K14
U 1
′
M x K12 K22 K23 K24 ϕ
(2.22)
=
′′
M
y
K13 K23 K33 K34
U 3
U2 ′′
Mz
K14 K24 K34 K44
For double cell analysis an aerofoil section is considered. Generally an aircraft wing
has two spars, first is located at around 15% of chord and the last spar at around 65%
of chord. Two cell box analyses are used in those sections. For the initial box distance
between 5% of chord and the front spar (15% of chord) is taken and for the second cell
distance between front spar and rear spar (65% of chord) are assumed as shown in the
Figure 2.4. For spanwise distribution of the stiffnesses EI, GJ and K, the methodology
used for single cell and depicted by the schematic diagram in Figure 2.3 is used for the
double cell analysis when averaging the properties between two subsequent sections.
To analyse, a composite box having 7 laminate parts is used. At each of these laminates
the ply orientation of composite is generally varied between 0 to 90 degrees and for each
30 Chapter 2. Metallic and composite wing boxes
of these ply angle, the bending (EI), torsional (GJ) and coupling (K) stiffnesses are found.
The composite wing laminate is made up of 8 plies, each ply is having a thickness of with
0.125 cm at the root. The wing thickness will vary from root to tip linearly similar to
single cell analysis and the double cell box has been split into 7 nodes as shown in Figure
2.5. Several cases have been considered, all having a ply orientation varying from 5 to 90
degrees with an increment of 5 degrees. Among them the following case is found to be
the most suitable one to produce desired aeroelastic effect based on its stiffness properties.
For the selected case, 3rd , 4th , 5th and 7th laminates have negative ply orientation while
the rest of the laminates have positive ply orientations. The elastic constraints used for
the composite material are that of carbon-fibre reinforced plastics (CFRP) and are given
in Table 2.2.
G12 = G13
E1 (N/m2 ) E2 (N/m2 ) Poisson’s Ratio (ν12 ) Density, ρ (kg/m3 )
= G23 (N/m2 )
0.14E12 0.95E10 0.28 0.58E10 1562
Initially wing is split into 6 sections which are of varying thickness similar to that of a
single cell analysis. Using an aerofoil, 6 aerofoil shapes at those 6 sections are deter-
mined. The coordinates for these aerofoil sections at each of the nodes are calculated
using AutoCAD. For these 6 sections, coordinates at 5%, 15% and 65% for lower aero-
foil surface and upper aerofoil surface are calculated by having origin from the nose of the
aerofoil. Then the origin is moved to the middle of 7th laminate (between nodes 1 and 4),
then based on the new centre coordinates the position of the other 6 nodes are calculated.
Thus the coordinates at various laminates are calculated for the wing span of 15 m, then
stiffness analysis is carried out.
2.1. Introduction to metallic and composite structures used in aircraft design 31
Since the thickness is small, thin walled assumptions are used in the analysis and the
following formulae is used to generate the results.
Area,
Area = t × {(a + b + c) + l1 + d + 2e} (2.23)
Mass,
mass = area × density = A × ρ (2.24)
Moment of inertia about x- axis,
!#2 !#2
a3 + b3 + c3
" "
b − c b − a
+ 2 e × c + (2.25)
I xx = t × + d × a +
12 2 2
For aluminium, the stiffness properties such as the bending (EI) and torsional (GJ) stiff-
nesses are calculated using the following formulae
E
G= (2.29)
2.6
4s2 4bav 2 dav 2 t
J = H ds = (2.30)
t
2 (bav + dav )
32 Chapter 2. Metallic and composite wing boxes
bouter douter 3 − binner dinner 3
I = Iouter − IInner = (2.31)
12
bav dav 3
Iav = (2.32)
12
Figure 2.6: A representation of sections obtained for vibrational analysis from an wing.
The flutter behaviour of composite wings can be very different from their metallic coun-
terparts. The analysis for composite wing is significantly influenced by the ply orientation
in the laminate. It should be also noted that sweep angle also plays an important role in
flutter analysis for both metallic and composite wings. The wing is idealized as a se-
ries of composite beams whose bending (EI), torsional (GJ) and coupling (K) stiffnesses
have already been established. Using these stiffness properties the natural frequencies
and mode shapes of the cantilever wing built-in at the root are computed by applying the
method of the dynamic stiffness matrix. First, the free vibrational analysis is carried out
for laminated composite wing of rectangular box cross sections with single cell represen-
tation having unidirectional ply orientation (θ) for each of the four laminates comprising
the wing box. However, for double cell representation which comprises seven laminates,
2.1. Introduction to metallic and composite structures used in aircraft design 33
the ply orientations for each of the laminates were different as shown in Table ??. When
obtaining the results, the ply orientations were varied from 5 to 90 degrees in steps of 5
degrees. For presentational purposes and also to avoid the clumsiness of the figures, it
was adjudged necessary to show the results for 10 degrees interval of ply orientations.
Figures 2.7 and 2.8 illustrates the bending (EI) and torsional (GJ) stiffnesses along wing
spanwise direction, respectively for various ply orientations when using the single cell
theory for the wing described in the previous section. Stiffness distribution for an alu-
minium wing of same dimensions is also shown in the Figures 2.7 and 2.8 so that a direct
comparison of stiffnesses are possible. From Figure 2.7, it can be ascertained that the
composite wing has higher bending stiffness than the aluminium one when the ply angle
orientation is below 20 degrees. This is to be expected because the effective Young’s mod-
ulus will be at its maximum at 0 degree or nearer to 0 degree ply orientations, particularly
for the wing box configuration under consideration. Clearly, the Figure 2.7 reveals that
for ply orientations above 25 degrees aluminium has better bending stiffness properties.
This is obvious because with increasing angle of ply orientation, the Young’s modulus
in the matrix direction which is generally small comes more and more into play in the
bending stiffness (EI) formulation. It should be noted that as ply orientation is increased
the bending stiffness decreases gradually as expected. From Figure 2.8, it can be seen
that aluminium wing always has higher torsional stiffness than the composite wing. This
is also to be expected because the torsion constant J, for a closed section such as the
rectangular box will have higher values for the aluminium than the composites. Figure
2.9 which is relevant only to composite wing as its shows the variation of the bending-
torsional material coupling parameter, K (which is non-existent in isotropic materials such
as aluminium) along the span. It is significant to note that the parameter K can be manip-
ulated to advantage in wings of composite construction which, of course, is impossible in
aluminium wings. It is observed in Figure 2.9 that the values of K around 15 and 25 de-
grees of ply orientation provide the maximum material coupling effect. Based on this two
angles of ply orientations, bending (EI), torsion (GJ) and the bending-torsion material
coupling (K) stiffnesses are plotted in Figure 2.10. In the next stage of the investigation,
the free vibration behaviour of the composite wing is studied in detail for different ply
orientations.
Table 2.3 shows the first five natural frequencies computed by varying the ply orientation
between 5 and 90 degrees in steps of 5 degrees. The bending, torsion and coupled modes
are designated by the letters B, T and C respectively. From Table 2.3, it can be observed
34 Chapter 2. Metallic and composite wing boxes
that the natural frequencies in general decrease with increasing ply angles. An interest-
ing and intriguing features of the results indicate the modal interchanges or flip-overs are
possible in composite wings of the same mass and the geometry of constructions. This
is certainly not possible for wings with metallic constructions. It can also be noted in
Table 2.3, that the fundamental node is always a bending one, whereas the second and
third mode can be changed from couple to bending or from torsion to coupled or from
coupled to bending depending on the ply orientation. However, the fifth mode is basically
a torsional mode for almost all ply orientation.
Figure 2.11 illustrates representative mode shapes corresponding to ply orientations of
2.1. Introduction to metallic and composite structures used in aircraft design 35
Figure 2.10: Stiffness properties at ply orientation of 15 and 25 degrees for single cell.
15 and 35 degrees with the solid black line showing bending displacements and the red
broken lines showing the torsional rotations. An inspection of the mode shapes in Figure
2.11 suggests that for 15 degrees ply angles the fundamental mode is bending, whereas
the second and third mode are bending-torsion coupled. For the same ply angle, the fourth
mode is a pure torsional mode whereas the fifth one is a coupled mode. Now turning at-
tention to the mode shapes corresponding to 35 degrees, somehow some similar as well as
different observations for the five modes can be made. For instance, the first three modes
follow more or less the same pattern, but the fourth mode for the 35 degrees case is a
36 Chapter 2. Metallic and composite wing boxes
coupled one unlike the 15 degrees case for which it was a pure torsional mode. The fifth
mode for the 35 degrees case is pure torsional mode whereas it was a coupled mode for
the 15 degrees case.
Using the modes described in Table 2.3 flutter analysis was carried out on the single cell
composite wing and the results are shown in Table 2.4, for the entire range of ply ori-
entations. The maximum achievable flutter speed appears to be at the higher end of ply
angles. The reason for this can be attributed to the fact that the fundamental bending fre-
quency goes down significantly at high ply angles whereas the reduction in fundamental
torsional frequency will not be that much and as a consequence, the separation between
the fundamental bending and torsional natural frequencies becomes much wider which
predictably increase the flutter speeds.
Table 2.3: Natural frequencies and characterisation of modes for various ply orientation
of a single cell composite wing. (B: Bending dominated mode; T: Torsion dominated
mode; C:Bending-torsion coupled mode.)
For an accurate representation of the stiffness properties of the composite wing, it is often
instructive to use double cell representation of the wing cross section (see Figure 2.4) as
opposed to a single cell one (see Figure 2.2). Therefore, the results for stiffnesses, nat-
ural frequencies, mode shapes, flutter speeds and flutter frequencies in this section are
computed based on double cell idealisation. Figures 2.12, 2.13 and 2.14 illustrates the
bending (EI), torsional (GJ) and bending-torsional coupling (K) stiffness distributions,
respectively along the span of the composite wing for different ply orientations using
double cell wing box idealisation. As can be seen in Figure 2.12, lower values of ply
orientations yield higher bending stiffness properties (EI) as was the case with single box
idealisation. The reason for this is of course due to the enormous strength of the fibres
as opposed to the matrix for lower values of ply orientations. Interestingly, Figure 2.13
shows that the torsional (GJ) stiffness is at its highest when the ply orientations is at 25
degrees whereas Figure 2.14 reveals that for both 15 and 25 degree ply orientations the
bending and torsional coupling (K) stiffness will reach its maximum value. The above
2.1. Introduction to metallic and composite structures used in aircraft design 39
results led to the presentation of Figure 2.15 which shows the interesting results for EI,
GJ and K corresponding to 15 and 25 degrees ply orientations in one graph.
Based on the stiffness properties illustrated in Figures 2.12 to 2.14, the natural frequen-
cies and mode shapes of a double cell composite wing was computed for further flutter
analysis.
Figure 2.15: Stiffness properties at ply orientation of 15 and 25 degrees for double cell.
Table 2.5 shows the first five natural frequencies computed by varying the ply orientation
between 5 and 90 degrees in steps of 5 degrees. The bending, torsion and coupled modes
are designated by the letters B, T and C respectively. It is observed that with increase in
ply orientations the natural frequencies decrease as was the case with single cell idealisa-
tion. However, unlike the single cell idealisation Table 2.5 shows that the first two modes
are predominantly bending modes, whereas the third mode can be either predominantly
coupled or bending. By contrast, the fourth mode is either torsional or coupled as can be
2.1. Introduction to metallic and composite structures used in aircraft design 41
seen. Interestingly, the fifth mode can be either coupled or torsional or bending and no
predictable pattern is evident.
Figure 2.16 shows representative mode shapes for ply orientation of 15 and 25 degrees for
the double cell composite wing with the solid black line showing bending displacements
and the red broken lines showing the torsional rotations. As mentioned before values ob-
tained are not similar to single cell because of the laminate layup chosen but the trends
in all the results follow closely with those of single cell member. For Figure 2.16, at ply
orientation of 15 degree it is observed that mode 1 is bending, modes 2 and 3 are strongly
coupled, mode 4 is torsion dominant, mode 5 is strongly coupled. Now on observing
mode shapes corresponding to the ply orientation of 35 degrees it is observed that mode 1
and 2 are bending dominant, mode 3 is strongly coupled, mode 4 is torsion dominant and
mode 5 is coupled. It is evident that as ply orientation changes bending and torsion will
be varied. Using the modes described in Table 2.5, flutter analysis was carried out on the
double cell composite wing and the results are shown in Table 2.6, for the entire range of
ply orientations. The maximum achievable flutter speed appears to be at the middle range
of ply angles.
42 Chapter 2. Metallic and composite wing boxes
Table 2.5: Natural frequencies and characterisation of modes for various ply orientation
of a double cell composite wing. (B: Bending dominated mode; T: Torsion dominated
mode; C:Bending-torsion coupled mode.)
In order to understand the analysis for free vibration and flutter of an metallic wing, an air-
craft wing made up of aluminium is used. Geometry of the metallic wing box is similar to
those of the wing box used for single and double cell analysis. The mode shapes obtained
from the analysis is similar in shape to the composite wing box analysis but the natural
frequency values are much higher for the metallic wing.It should be noted that standard
formulae described in section 2.1.1.3 are used to perform stiffness analysis. Figure 2.17
gives the mode shape for an aluminium wing with similar dimension used to obtained
single cell and double cell results. Table 2.7 shows the flutter speed and flutter frequency
obtained for metallic wing. It can be seen that the flutter speed of 510 m/s is obtained for
metallic wing.
2.2. Conclusions 45
Aluminium
Flutter speed (m/s) 510.0
Frequency (rad/s) 28.94
2.2 Conclusions
It can be concluded that by using laminated composite wings, aeroelastic properties can
be varied based on the user’s criteria. For metallic wing box, since there can be no con-
trol in the stiffness properties, its flutter speed cannot be varied while for composites it
can be changed using different ply orientation. Thus desired flutter results can be ob-
tained by varying ply orientations and also by using double cells improved results can
be achieved than by using single cell. Only geometric coupling is possible in metallic
structure while both geometric coupling and material coupling are possible in compos-
ites. Bending–torsional coupling, which is not possible in metallic wing, can be varied
for the composite wing to achieve desired aeroelastic effect. Also its effect on the flutter
46 Chapter 2. Metallic and composite wing boxes
analysis was discussed in detail in the next chapter. Since the analysis is carried on the
swept wing which gets more aeroelastic effect than that of rectangular wing, this result
is also applicable to unswept wings. Also, it should be noted that the thin-walled section
theory used in this chapter to perform single cell and double cell analysis which neglects
the effect of shear contribution and the assumptions about normal and shear stresses being
constant is also not applicable for 3D wing.
Chapter 3
Geometric
Value
parameter
Aspect ratio 10
Wing area (m2 ) 92.25
Wing span (m) 30.42
Root chord (m) 5.30
Tip chord (m) 1.34
Mean aerodynamic chord (m) 3.363
ples by utilising and suitably adopting wherever possible the stiffness data obtained from
BOXMXES and FEMAP/NASTRAN.
It should be noted that there were some marked variations in the stiffness properties from
the two types of analysis described above and this could be due to the complexities of
the wing design, particularly due to the presence of the manholes in the lower skin of the
wing. Also, whilst the practical 3D wing model in FEMAP/NASTRAN includes cut-outs,
sweep, taper ratio, flexible ribs, the 2D BOXMXES model uses only the cross section of
the wing in a 2D plane. This necessitated a further investigation of parametric nature
involving the wing box geometry for the stiffness analysis. Therefore, at the end of this
chapter, an investigation is undertaken to understand and provide an insight into the pa-
rameters affecting the wing box stiffnesses.
Figure 3.1: Composite wing model with 26 sections represented by 7 skin parts (see
Table 3.4). Skin part 1: sections 1-3; Skin part 2: sections 4-7; Skin part 3: sections 8-11;
Skin part 4: sections 12-15; Skin part 5: sections 16-19; Skin part 6: sections 20-22; Skin
part 7: sections 23-26.
First the composite wing model which has 27 ribs along semi-span is considered and it
was divided into 26 sections, numbered from tip to root, with each section chosen between
two successive ribs so that they lie across the centre line of each manhole, see Figure 3.1.
Note that, section 1 to section 19 (outboard wing) consist of only two spars, front spar and
rear spar. They are considered as single cell box sections for the analysis. On the other
50 Chapter 3. Wing analysis and parametric investigation
hand, section 20 to section 26 (inboard wing) consist of three spars, front spar, middle spar
and rear spar. They are considered as double cell box sections for the analysis. Section 1
represents the section near the wing tip whereas section 26 represents the section near the
wing root. The composite layup orientation and stacking sequence are given in Appendix
B. The thickness for each laminate for the top and bottom skin of the wing is taken as
1.83 × 10−4 m. It should be noted that wing model analysed has manhole openings on the
lower skin between two successive ribs which complicates the stiffness analysis.
For the given wing, BOXMXES program needs the coordinates of each part divided
around the circumference of the box section which is difficult to obtain directly from
FEMAP/NASTRAN and furthermore, the BOXMXES prefers the wing box section to be
roughly symmetrical about the horizontal axis. Therefore, as an acceptable alternative, the
dimensions needed for the wing box sections are determined from FEMAP/NASTRAN
model by taking proper account of the distances between nodes appearing along upper
skin, lower skin, front spar, rear spar (for sections 1-19) and middle spar (for sections
20-26). Upper skin and lower skin were further subdivided between two stringers. AU-
TOCAD is then used to plot these parts to recreate the cross-sections. By using AUTO-
CAD, coordinates of each part were identified along the circumference of the wing box
sections. Using the cross sectional coordinates, element properties, material properties,
laminate layup and stacking sequence, two input data file for BOXMXES program was
created, one for the single cell and the other is for the double cell. Then using BOXMXES
program bending (EI) and torsional (GJ) stiffnesses are calculated. To facilitate this work,
MATLAB is used to determine the preliminary estimates of EI and GJ values which are
subsequently refined and improved. It is to be noted that the BOXMXES program does
not take into account the contribution of the stringers. So by using the parallel axis theo-
rem, effect of the stringers was taken into account by lumping the appropriate areas. The
output from the BOXMXES program gives the stiffness properties, EI, GJ and the effec-
tive elastic constants E x , Ey and G xy . For the sections 1-26 the effective elastic constants
E x , Ey and G xy are given in Table 3.4. The elastic constants E x , Ey and G xy obtained
from BOXMXES program are then fed into FEMAP/NASTRAN to determine the stiff-
ness properties EI and GJ. The bending-torsion coupling (K) was not calculated using the
FEMAP/NASTRAN because of the difficulties encountered to adapt the package for this
purpose. These stiffness properties, EI and GJ are then used to create the data needed for
CALFUNB and COMPCAL programs to compute the free vibrational modes and flutter
speed.
3.3. Theoretical and numerical procedures 51
Skin part 1
Section 1,2,3 Ex Ey Gxy
Upper skin,
Lower Skin, 0.510E11 0.636E11 0.219E11
Rear Spar
Front spar 0.518E11 0.518E11 0.254E11
Skin part 2
Section 4,5,6,7 Ex Ey Gxy
Upper skin,
Front Spar, 0.576E11 0.576E11 0.222E11
Rear Spar
Lower Skin 0.569E11 0.569E11 0.213E11
Skin part 3
Section 8-11 Ex Ey Gxy
Upper skin,
0.620E11 0.529E11 0.220E11
Lower skin,
Front spar,
0.529E11 0.620E11 0.220E11
Rear Spar
Skin part 4
Section 12-15 Ex Ey Gxy
Upper skin, 0.554E11 0.629E11 0.212E11
Front spar, 0.572E11 0.572E11 0.217E11
Lower skin,
0.559E11 0.559E11 0.230E11
Rear Spar
Skin part 5
Section 16-19 Ex Ey Gxy
Upper skin, Front spar 0.529E11 0.594E11 0.230E11
Lower Skin,
0.591E11 0.526E11 0.227E11
Rear Spar
Skin part 6
Section 20-22 Ex Ey Gxy
Upper skin, Front Spar,
0.562E11 0.562E11 0.229E11
Lower skin
Rear Spar,
0.621E11 0.502E11 0.229E11
Mid wall
Skin part 7
Section 23-26 Ex Ey Gxy
Upper skin, Lower skin,
0.591E11 0.537E11 0.229E11
Rear spar, Mid wall
Front Spar 0.589E11 0.535E11 0.226E11
52 Chapter 3. Wing analysis and parametric investigation
Figure 3.2: Applied bending moment and torque for an cantilevered wing box.
To calculate the wing box bending stiffness, Euler-Bernoulli beam theory is used. Thus
the moment-curvature relationship is given by
′′
M (l) = EIk = EIw (y) (3.1)
where M is the bending moment, EI is the bending stiffness, k is the curvature and w is
the vertical displacement of the section and a ’prime’ denotes differentiation with respect
to the span wise displacement y.
For a bending moment M applied at the tip of a typical representative sectional length L,
it can be established that the bending stiffness EI follows the following relationships
ML
EI = (3.2)
θ (L)
3.4. Stiffness evaluation using finite element model 53
where θ(L) is the bending slope in radian at the tip of the section analysed.
In a similar manner, the torsional stiffness GJ is computed by applying a constant torque
T at the tip of the outboard section. Using elementary torsion theory the twist ψ(L) at the
tip of the section can be related to the torque T and torsional stiffness GJ as follows
TL
ψ (L) = (3.3)
GJ
Equation 3.3 can be rearranged to give the torsional stiffness GJ by the following equation
TL
GJ = (3.4)
ψ (L)
Table 3.5: Applied bending moment and bending rotation of each section.
The data given in Table 3.5 which were obtained through the application of FEMAP/NASTRAN
are further utilised to estimate the bending stiffness EI of each section by applying Equa-
tion 3.2. The results are shown in Table 3.6.
54 Chapter 3. Wing analysis and parametric investigation
Based on the underlying methodology described in section 3.4 the values of the applied
torque and the corresponding twist of each of the 26 sections are shown in Table 3.7 which
were computed by using FEMAP/NASTRAN.
The data given in Table 3.7 which were obtained through the application of FEMAP
/NASTRAN are further utilised to estimate the torsional stiffness GJ of each section by
applying Equation 3.4. The results are shown in Table 3.8.
3.4. Stiffness evaluation using finite element model 55
Table 3.7: Applied torque and twist of the wing box sections.
Free vibration and flutter analyses are carried out using both CALFUNB and COMPCAL
which are the derivatives of the original program CALFUN. These two programs are
completely independent in that CALFUNB is developed for metallic wings which does
not require the bending-torsion material coupling stiffness (K) in the input data whereas
COMPCAL is developed for composite wing requiring the stiffness K. However, in order
to verify the CALFUNB results a small values of the bending-torsion coupling parameter
K of the order 10−4 was used in COMPCAL as a degenerate case.
For ease of computation, the 26 wing box sections are reduced to 13 sections (see Table
3.9) by using the appropriate EI and GJ values obtained from FEMAP/NASTRAN. The
input data file for CALFUNB and COMPCAL were then created. The data preparation
included determining mass per unit length, polar mass moment of inertia per unit length,
distance between mass and element axis, projection of the element length along X and Y
axes, semi chord, etc. From the input file thus created, the in-house programs CALFUNB
and COMPCAL were activated and the mode shapes and flutter speed were computed.
The first six natural frequencies and mode shapes computed using CALFUNB and COM-
PCAL are shown in Figure 3.3 in which the solid black line shows bending displacement
and the broken red line shows torsional displacement. Predictably the program COMP-
CAL produced more or less the same natural frequencies and mode shapes as CALFUNB.
The modes generated by CALFUNB or COMPCAL needs some discussion. Referring to
Figure 3.3, mode 1 of the wing is basically the fundamental bending mode. The second
mode is also bending with very little torsion in it, whereas third mode is coupled in bend-
ing and torsion. The fourth and fifth modes are also to all intents and purposes couple
modes whereas the sixth mode can be regarded as a pure torsional mode. These 6 modes
are subsequently used in the flutter analysis and the results are shown in Table 3.10 which
shows the flutter speeds and flutter frequencies obtained using CALFUNB and COMP-
CAL for FEMAP/NASTRAN and BOXMXES.
3.4. Stiffness evaluation using finite element model 57
Although reasonably satisfactory results for flutter speeds and flutter frequencies are ob-
tained from the stiffnesses calculated using FEMAP/NASTRAN, it was felt to be neces-
sary to revisit the estimation of stiffness properties using the BOXMXES and FEMAP
/NASTRAN programs. Clearly, there were significant differences in the stiffness proper-
ties of the wing computed by the 2-D and 3-D idealisations using these two very different
programs. Table 3.11 shows the bending stiffness EI computed by the two programs with
percentage difference in the stiffnesses for the 26 sections varied from 10% to 35%. The
differences are unacceptably large. Similar observations are made for the torsional stiff-
ness (GJ) distributions for which the percentage difference was even greater for some
cross sections, see Table 3.12. The stiffnesses computed by BOXMXESC were larger
than the ones computed by FEMAP/NASTRAN and the reason may be attributed to the
fact that the former does not includes the effects of manhole and cut-outs. This prompted
a further investigation which includes:
(i) the presence of the manholes (cut-outs for access panels) appearing in the lower skin,
(ii) the rigidity of the ribs to retain the aerofoil shape,
(iii) the taper ratio of the wing planform and the taper ratio of the depth of the wing from
root to tip
(iv) the sweep angle of the leading edge.
60 Chapter 3. Wing analysis and parametric investigation
Table 3.11: Bending stiffness, EI of the wing boxes between sections 1-26.
EI from EI from
Percentage
Sections FEMAP/NASTRAN BOXMXE
difference
(Nm2 ) (Nm2 )
1 2.59E+06 3.31E+06 21.75
2 3.24E+06 4.03E+06 19.60
3 3.90E+06 4.79E+06 18.58
4 5.27E+06 6.35E+06 17.01
5 6.20E+06 7.49E+06 17.22
6 7.07E+06 8.77E+06 19.38
7 8.85E+06 1.07E+07 17.37
8 1.28E+07 1.50E+07 14.67
9 1.48E+07 1.71E+07 13.45
10 1.71E+07 1.95E+07 12.31
11 1.94E+07 2.21E+07 12.22
12 2.32E+07 2.61E+07 11.11
13 2.61E+07 2.93E+07 10.92
14 3.04E+07 3.39E+07 10.32
15 3.57E+07 3.98E+07 10.30
16 4.23E+07 4.71E+07 10.19
17 4.67E+07 5.13E+07 8.97
18 4.81E+07 5.43E+07 11.42
19 5.10E+07 7.05E+07 27.66
20 9.58E+07 1.10E+08 12.91
21 1.14E+08 1.53E+08 25.49
22 1.57E+08 2.13E+08 26.29
23 2.10E+08 3.05E+08 31.15
24 2.59E+08 3.95E+08 34.43
25 3.18E+08 4.94E+08 35.63
26 4.02E+08 6.33E+08 36.49
3.4. Stiffness evaluation using finite element model 61
Table 3.12: Torsional stiffness, GJ of the wing boxes between sections 1-26.
GJ from GJ from
Percentage
Sections FEMAP/NASTRAN BOXMXE
difference
(Nm2 ) (Nm2 )
1 8.87E+05 1.26E+06 29.60
2 1.52E+06 1.64E+06 7.32
3 2.12E+06 2.09E+06 -1.44
4 3.18E+06 3.11E+06 -2.25
5 4.04E+06 3.87E+06 -4.39
6 5.10E+06 4.74E+06 -7.59
7 6.60E+06 5.68E+06 -16.20
8 9.48E+06 7.98E+06 -18.80
9 1.13E+07 9.38E+06 -20.47
10 1.33E+07 1.10E+07 -20.91
11 1.56E+07 1.28E+07 -21.88
12 2.16E+07 1.79E+07 -20.67
13 2.49E+07 2.05E+07 -21.46
14 2.96E+07 2.34E+07 -26.50
15 3.41E+07 2.65E+07 -28.68
16 4.38E+07 3.48E+07 -25.86
17 5.57E+07 3.86E+07 -44.30
18 6.63E+07 4.44E+07 -49.32
19 7.68E+07 5.59E+07 -37.39
20 1.79E+08 1.57E+08 -14.01
21 2.04E+08 2.07E+08 1.45
22 2.75E+08 2.63E+08 -4.56
23 3.29E+08 3.88E+08 15.21
24 3.02E+08 4.84E+08 37.60
25 3.33E+08 5.86E+08 43.17
26 3.36E+08 7.21E+08 53.40
Table 3.13: Rib rigidity effect on the bending stiffness of box section 16.
that the BOXMXE results do not allow for the cross-sectional warping which essentially
means that the ribs have been assumed to be rigid undergoing no in-plane or out of plane
displacements. By contrast, the flexibility of the ribs are inherently accounted for in the
FEMAP/NASTRAN model. In order to bring parity in the sets of results an investigation
is carried out to examine the effect of the rib-rigidity on the bending and torsional stiff-
nesses for section 16 of the composite wing. This investigation is particularly relevant
to the analysis using FEMAP/NASTRAN, but not so for the BOXMXE program because
for the latter the question of rib rigidity does not arise as the theory used in BOXMXE
inherently assumes that the section does not warp implying that the rib is basically rigid.
As explained in the description of the wing analysed in section 3.2, the aircraft wing was
split into 26 sections numbered from tip to root. For the analysis, section 16 is considered
mainly because it has single cell and the parametric study using single cell box is preferred
for simplicity. Also section 16 is a critical section which lies closer to pylon in the kink
region. The results of this investigation are shown in Tables 3.13 and 3.14 for the bending
and torsional stiffnesses respectively. It is clear that the 3D results for the bending stiff-
ness EI from FEMAP/NASTRAN model using rigid rib are in close agreement with the
2D results with the effect of cut-out computed from BOXMXE. Apparently the flexibility
of the rib made a significant difference to the bending stiffness and the probable expla-
nation could be the fact that the length of the box element chosen was quite small and
comparable to its width, bearing in mind that the application of load to such a small ele-
ment is expected to have substantial localized effect. With regard to the torsional stiffness
GJ, Table 3.14 shows the result using FEMAP/NASTRAN and BOXMXE. Clearly the
comparative results show that the rigid rib assumption from the two programs is produc-
ing more differences in torsion stiffness GJ than the observed differences in the bending
stiffness EI. This suggests that a thorough but detailed investigation on the estimation of
the torsional stiffness GJ is needed to ascertain the suitability using FEMAP/NASTRAN
and BOXMXE model when establishing the torsional properties of the composite wing.
Therefore, an in-depth investigation to determine the effects of various parameters on the
torsional strength of wing sections is carried out in the subsequent section.
3.4. Stiffness evaluation using finite element model 63
Table 3.14: Rib rigidity effect on the torsional stiffness of box section 16.
The previous section is concerned with the effects of rib rigidity on the bending and tor-
sional stiffnesses of a typical composite wing section and as it turned out, the bending
stiffnesses can be computed reasonably accurately using the rigid rib assumptions, but
the estimation of torsional stiffnesses requires further investigations. Therefore, the pri-
mary focus is now confined to the examination of the effects of various parameters on
the torsional stiffness (GJ) of wing sections. One of the main problems in determining
the torsional stiffness using the 2D analysis through BOXMXES for a wing section with
cut-out stems from the fact that the relevant portion of the wing box cross section con-
taining the cut could not be removed because the section would then become an open
section which is not the real case to represent the actual wing as the cut-out does not ex-
tend through the neighbouring sections along the span wise direction of the wing. The
problem is further compounded by a number of other reasons which are explained next.
Altogether the following parameters are included in the investigation: (i) the presence of
the manholes (cut-outs for access panels) appearing in the lower skin, (ii) the rigidity of
the ribs to retain the aerofoil shape, (iii) the taper ratio of the wing planform and the taper
ratio of the depth of the wing from root to tip and (iv) the sweep angle of the leading edge.
In order to obtain an accurate measure of the torsional stiffness GJ and to make some en-
gineering judgement on the influence of each of the above parameters on GJ, the in-depth
investigation was undertaken by using FEMAP/NASTRAN together with the application
of the classical theory of aircraft structures. The results from the above two methods (one
numerical and the other theoretical) are compared and contrasted with some discussion.
It should be recognised that in the numerical method using FEMAP/NASTRAN, the tor-
sional stiffness GJ is obtained for a typical element of certain length by cantilevering the
element at its one end and applying a pure bending moment or a pure torque (pure, as
closely as possible) at the other end. This inevitably involves the length of the element
as a parameter in the data. By contrast, in the theoretical method the torsional stiffness is
simply a matter of cross sectional parameter without involving the length of the element
in the data. A part of the numerical investigation simulated the non-warping of the cross
section by assuming the ribs to be rigid. A number of case studies have been conducted
with primary emphasis on the torsional stiffness GJ as indicated and these are explained
64 Chapter 3. Wing analysis and parametric investigation
next. It is to be noted that the case studies from 1 through 8 correspond to wing sections
of metallic constructions (i.e. aluminium alloy) with no cut-outs (manholes) whereas the
case studies 9 and 10 are focused on the shape, size and material of composite wing of
similar dimensions of Section 16 with and without cut-outs (manholes). For comparative
purposes case study 9 included parallel investigations using aluminium material alongside
the original composite material of the wing whereas case study 10 is solely for composite
material. For each of the case studies, a general description showing the relevant dimen-
sions is illustrated in Figure 3.4 and the corresponding data for the dimensions of a to l
are shown in Table 3.15. Further details of each of the case studies are given in Appendix
C.
The definitions of taper ratio and its values for each of the relevant cases which have taper
are given in Table 3.16. The leading edge is considered to be on the right hand side of the
box whereas the trailing edge is on the left hand side (see Figure 3.4). Also the front of the
box is considered to be the tip and the rear to be the root in the analysis when imposing
the boundary conditions.
Figure 3.4: A general description of principal dimensions covering all 10 case studies.
3.4. Stiffness evaluation using finite element model 65
Table 3.16: Taper ratio for case studies referring to Figure 3.5 and Table 3.15.
In this case study, a rectangular metallic wing box of constant cross section and a given
element length is considered in the finite element analysis using FEMAP/NASTRAN to
establish the torsional stiffness GJ (see Figure C1 of Appendix C). The left hand end of
the element is cantilevered and the front face is loaded by two equal and opposite forces
in the front and the rear so as to produce a pure torque about the shear centre and thus
making sure that the loading caused no bending displacement of the cross section, but only
twisting deformation. A certain amount of trial and error was required to produce a pure
twist. The torsional stiffness GJ is then evaluated by using standard procedure, relating
the applied torque to the rate of twist of the element. In this case study, the in-plane
displacements of the front and rear faces are both restrained first and then unrestrained
later when computing the torsional stiffness. However, this restriction does not prevent
the torsional rotation. This essentially means that the ribs are rigid in the former and
flexible in the latter. For theoretical calculation of the torsion constraint J which requires
only the 2-D cross sectional details, the rear and the front face data are used to obtain GJ
and an average value was worked out at the midpoint so as to make the results somehow
comparable with the 3-D FEMAP/NASTRAN results. This strategy is consistently used
for the rest of the case studies.
Case study 2 is essentially the same as case study 1 except that a refined mesh in FEMAP/NASTRAN
is used to ensure that the convergence of results. Further refinement of the mesh did not
alter the results to any appreciable extent. Therefore, the mesh in case study 2 is consis-
tently used to ensure the accuracy of the results. Figure C2 of Appendix C corresponds
to case study 2 and they are merely a duplication of Figure C1 of Appendix C with the
understanding that case study 2 is simply based on refined mesh generation.
As the wing planform tapers, case study 3 principally focuses on the effect of taper ratio
on the torsional stiffness GJ of the wing. In this particular case, only the wing planform
taper is considered whereas the depth-wise taper in the span wise direction which is gen-
erally small is not taken into account by assuming the depth to be constant (see Figure
C3 of Appendix C). As in the case of case studies 1 and 2, the beam element of Figure
3.4. Stiffness evaluation using finite element model 67
C3 of Appendix C for case study 3 is still metallic and its rear end is cantilevered when
applying a pure torque on the front face about a point which is as close to the shear centre
as possible. (As before, this is achieved by trial and error method ensuring that the ap-
plication of the load produces only torsional rotation and not any bending displacement).
The wing planform tapers for both leading and trailing edges for this case as can be seen
in Figure C3 of Appendix C.
Case study 4 is similar but different from case study 3 in the sense that the wing planform
is not doubly tapered, but singly tapered linearly on the front which is essentially the
leading edge as shown in Figure C4 of Appendix C. The trailing edge is assumed to be
straight as shown. In this way, the case study 4 focuses on the effect of sweep on the
torsional stiffness (GJ) of the section. As in previous cases, the beam element of Figure
C4 of Appendix C for case study 4 is still metallic and its rear end is cantilevered when
applying a pure torque on the front face about a point which is as close to the shear centre
as possible. Case study 4 has an unsymmetrical cross section which is evident from the
Figure C4 of Appendix C.
Case study 5 principally focuses on the effect of depth-wise taper ratio on the torsional
stiffness GJ of the wing. In this particular case, only the depth-wise taper is considered
whereas the wing planform remains constant (see Figure C5 of Appendix C). As in the
case of previous case studies, the beam element for case study 5 is still metallic and its
rear end is cantilevered when applying a pure torque on the front face about a point which
is as close to the shear centre as possible. The depth-wise change in dimension is taken
into account for both leading and trailing edges for this case as can be seen in Figure C5
of Appendix C.
Case study 6 is concerned with the effect of taper ratio for both leading edge and trailing
edge on the torsional stiffness GJ of the wing. In this case both the wing planform as
well as the depth-wise taper is in the span wise direction (see Figure C6 of Appendix C).
As in the case of previous case studies, the beam element of Figure C6 of Appendix C is
based on aluminium material. As before the rear end of the element is cantilevered when
applying a pure torque on the front face about a point which is as close to the shear centre
as possible.
68 Chapter 3. Wing analysis and parametric investigation
Case study 7 is similar but different from case study 6 in the sense that the trailing edge
is assumed to be straight whereas the leading edge is linearly tapered as shown in Figure
C7 of Appendix C. Obviously the case study 7 focuses on the effect of sweep only. Here
again, the beam element of Figure C7 of Appendix C for the case study 7 is still metallic
and its rear end is cantilevered when applying a pure torque on the front face about a point
which is as close to the shear centre as possible. Case study 7 has an unsymmetrical cross
section which is evident from the Figure C7 of Appendix C.
Case study 8 in many ways is similar but different from case study 4 in the sense that
dimensions are different even though the general layout is same. The trailing edge is
assumed to be straight whereas the leading edge is linearly tapered as shown in Figure C8
of Appendix C. As in previous cases, the beam element of Figure C8 of Appendix C for
case study 8 is still metallic and its rear end is cantilevered when applying a pure torque
on the front face about a point which is as close to the shear centre as possible. Note that
the case study 8 has an unsymmetrical cross section.
Case study 9 is based on similar dimensions of composite wing at section 16, but without
the presence of stringers and has four subcases: (a) Subcase 1 – wing box with cut-out
(manhole) made up of aluminium, (b) Subcase 2 – wing box without cut-out (no manhole)
made up of aluminium, (c) Subcase 3- wing box with cut-out made up of composite and
(d) Subcase 4 – wing box without cut-out made up of composite (see Figure C9 of Ap-
pendix C). All of these 4 subcases of the wing box are created using FEMAP/NASTRAN.
However, some simplifying assumptions are made in that the wing box dimensions are
considered to be made up of corner nodes of the actual aircraft wing box so that the wing
box surfaces are straight lines rather than curved lines whereas in the actual composite
wing box both the upper and lower skin surfaces are curved. The dimension and location
of the cut-out (manhole) are identical with the composite wing box. The relevant data
for case study 9 are given in Figure C9 of Appendix C. This case study focuses on the
effect of manhole on the torsional stiffness GJ for the composite wing and the differences
it makes when using isotropic material as opposed to composites. Note that for case study
9, the rib is considered to be rigid for both metallic and composite subcases.
3.4. Stiffness evaluation using finite element model 69
Case study 10 is a portion of the actual composite wing which is actually Section 16 of
the wing which unlike case study 9 includes the stringers. In essence, case study 10 has
all the essential features covered in previous case studies (see Figure C10 of Appendix C).
In all of the previous case studies the skin thicknesses were assumed to be constant, but
for this case study, the upper skin, lower skin and the front spar have the same thickness
whereas the rear spar and the ribs have different thicknesses that are representative of
the original composite wing model. However, case study 10 has also two subcases: (a)
Subcase 1 - wing box is made of composites and the ribs are considered elastic and made
of composites, (b) Subcase 2 - wing box is made of composites, but the ribs are rigid,
modelled by artificially increasing the elastic constants to a larger value. In the above two
subcases, the presence of the manhole has been taken into account.
The results of the first 8 case studies described above are shown in Table 3.7. The results
of Table 3.17 shows that for case studies 1 and 2, the discrepancy between the results
from FEMAP/NASTRAN and classical theory is quite small, notably within engineering
accuracy, but for case studies 3 and 4, the difference in results is quite substantial. The
results using classical theory were obtained by multiplying the (effective) shear rigidity
2
with torsion constant J evaluated from J = H4Ads with A being the cell area and the inte-
H t
gration dst is carried out all around the cross section. This could be attributed to the fact
that the taper ratio for these two cases is significant which can alter the results sufficiently.
Also, it should be noted that case 3 has doubly tapered and its GJ values are lot higher
for case 4 which was singly tapered. Similar trends were observed for case studies 5 to 8
which had the added complexity of sweep angle in addition to the taper ratio. From case
studies 6 and 7, it can be observed that the depth-wise taper does not contribute much to
the GJ difference. The big differences in results in part, are due to the small element size
chosen and the localise effect of the load applied when working out the stiffness. The GJ
results for case studies 1 to 8 are shown in Figure 3.5 in a histogram plot. Case studies
6 and 7 show large differences between the GJ values at the tip and the root due to the
differences in the dimensions at the two respective cross sections (see Figure C6 and C7
of Appendix C).
70 Chapter 3. Wing analysis and parametric investigation
Table 3.17: The computed torsional stiffness GJ using FEMAP/NASTRAN and the cor-
responding values using classical theory for cases 1 to 8.
Figure 3.5: Distribution of torsional stiffness GJ for the first 8 case studies.
The case studies 9 and 10 (see Figures C9 and C10 of Appendix C) corresponds to Section
16 of the composite wing for which the former is representative of Section 16 of the
composite wing but has no stringer attached and also the curved surfaces of the upper and
lower skins are represented by straight surfaces. (This is in contrast to case study 10 which
has all the components of the Section 16 of the composite wing both from geometrical
and material properties points of view). The results for case studies 9(a) to 9(d) are given
3.4. Stiffness evaluation using finite element model 71
Table 3.18: The computed torsional stiffness GJ using FEMAP/NASTRAN (Case 9).
FEMAP/NASTRAN results
Case study 9
GJ (Nm2 )
(a) With manhole (aluminium) 4.03E+07
(b) Without manhole (aluminium) 4.15E+07
(c) With manhole (composite) 4.95E+07
(d) Without manhole (composite) 5.05E+07
in Table 3.18 and Figure 3.6 which clearly indicate that the effect of cut-out (manhole) is
small both for aluminium and composite sections. These sections have similar dimensions
to those of Section 16 of the composite wing as mentioned except that the stringers have
not been taken into account in the stiffness analysis.
The final set of results for the torsional stiffness GJ concerning case study 10 (which
corresponds to an exact replica of Section 16 of the composite wing including the stringers
and cut-out) are given in Table 3.19 and Figure 3.7. The case studies 10(a) and 10(b) are
for Section 16 of composite wing with elastic and rigid ribs, respectively. Clearly the
results from case study 10 reveal that the FEMAP/NASTRAN stiffness results are very
different from the classical theory results, particularly when the results from FEMAP
/NASTRAN are based on rigid rib. The exact reason for this has not identified and the
matter needs further investigation. The probable cause could be in part due to the use of
72 Chapter 3. Wing analysis and parametric investigation
Table 3.19: The computed torsional stiffness GJ using FEMAP/NASTRAN and the cor-
responding values using classical theory for case study 10.
exceptionally small element and the application of the loading point to produce a pure
torque about the centre of twist or shear centre.
Figure 3.7: Torsional stiffness GJ for section 16 of the composite wing using FEMAP
/NASTRAN analysis and classical theory.
3.4. Stiffness evaluation using finite element model 73
Table 3.20: Bending, torsional and bending-torsion coupling stiffness data used in CAL-
FUN.
2
Bending and torsional stiffnesses Coupling stiffness K and α = EI×GJ
K
, (0<α<1)
2 2
EI (Nm ) GJ (m ) K(α=0) K(α=0.01) K(α=0.25) K(α=0.50) K(α=0.75)
1.96E+06 1.51E+06 0.0 1.72E+05 8.61E+05 1.22E+06 1.49E+06
4.39E+06 4.73E+06 0.0 4.56E+05 2.28E+06 3.22E+06 3.95E+06
9.31E+06 1.24E+07 0.0 1.08E+06 5.38E+06 7.61E+06 9.32E+06
1.85E+07 2.75E+07 0.0 2.25E+06 1.13E+07 1.59E+07 1.95E+07
2.70E+07 6.06E+07 0.0 4.05E+06 2.02E+07 2.86E+07 3.51E+07
9.48E+07 2.19E+08 0.0 1.44E+07 7.21E+07 1.02E+08 1.25E+08
2.11E+08 3.25E+08 0.0 2.62E+07 1.31E+08 1.85E+08 2.27E+08
Table 3.21: Bending-torsion coupling effects on flutter speed and flutter frequency.
3.5 Conclusions
It can be concluded that special attention needs to be paid to validate and improve the
model of calculating a wing box stiffness using BOXMXES based on 2D thin-walled
box beam theory. For validation purpose, a modelling method has been developed for
extracting the EI and GJ values from the FEMAP/NASTRAN model of the 3D wing
box. It should be noted that 2D model does not take into account manhole in its analysis.
Based on the results obtained, the maximum difference between the EI obtained from the
2D BOXMXES model and the 3D FEMAP/NASTRAN model is around 36%. It is also
noted that such large difference only occurs in the inboard wing made of double cell wing
box with large geometric variation. For the outboard wing made of single cell boxes, the
maximum EI difference is about 21%. The maximum difference between the GJ obtained
from the 2D BOXMXES model and the 3D FEMAP/NASTRAN model remains about
55% in the inboard wing. For the outboard wing, the maximum difference is around 28%.
From the parametric study of the geometric effect on the 2D model stiffness, it is clear that
the geometric features such as cut-outs, sweep, taper ratio, flexible ribs in the 3D of the
wing box have significant effects on the torsional stiffness, GJ, especially the presence
of taper along the wing planform contributes more to the torsional stiffness than other
parameters such as manholes, rigidity of ribs and sweep angle.
Chapter 4
and torsional stiffnesses of the original wings are altered between +25% and -25% insteps
of 5% and their subsequent effects on the modal behaviour of the wings and flutter analysis
are carried out. This is followed by further investigation wherein the engine mass and its
location wherever applicable, are varied and the free vibration and flutter behaviour is re-
examined. In each case, the wing is idealised as an assembly of bending-torsional coupled
beams for which the frequency dependent dynamic stiffness matrix is well established [12,
13]. The investigation needed considerable efforts for data preparation to model each of
the aircraft wings. Once the data preparation was completed, a detailed parametric study
with the variations of bending and torsional stiffnesses, the engine mass and its location
wherever applicable was undertaken and the free vibration and flutter analysis was carried
out on each of the wings. The results shows some interesting trends which are discussed
and commented on in this chapter.
Figure 4.4: Three different categories of aircraft represented in logarithmic scale (OWE
- operating empty weight, MTOW - maximum take-off weight).
Sailplane
Parameters
Sailplane-S1 Sailplane-S2
Wing Span (m) 22 15
Wing Area (m2 ) 15.44 10.05
Aspect Ratio 31.35 22.4
Wing Root Chord (m) 1.0 0.9
Wing Tip Chord (m) 0.4 0.4
Sweep angle (deg) 0 0
Length overall (m) 7.6 6.72
Height Overall (m) 2.0 2.0
Weight Empty (kg) 390 234
Max Take-off weight (kg) 550 440
Max Wing Loading (kg/m2 ) 37 36
Max Cruising Speed (knots) 135 105
4.2. Particulars of the aircraft considered for the analysis 79
Transport airliner
Parameters
T1 T2 T3 T4
Wing Span (m) 40 30 35 60
Wing Area (m2 ) 162 93 123 362
Aspect Ratio 10 9 10 10
Wing Root Chord (m) 5.0 5.5 6.0 10.5
Wing Tip Chord (m) 2.5 1.5 1.5 2.5
Sweep angle (deg) 0 28 28 28
Length overall (m) 30 36 38 60
Height Overall (m) 12 11 12 17
Weight Empty (kg) 34,000 26,000 42,000 130,000
Max Take-off weight (kg) 70,000 46,000 74,000 275,000
Max Wing Loading (kg/m2 ) 434 511 600 760
Max Cruising Speed (knots) 348 529 516 569
Range (nmi) 2835 2400 2592 8000
80 Chapter 4. High aspect ratio aircraft wings
The first five natural frequencies for all of the aircraft wings (S 1 , S 2 , L1 , L2 , T 1 , T 2 , T 3
and T 4 ) are given in Table 4.4. The letters B and T used in the table and elsewhere in
the thesis indicate bending and torsion dominated modes respectively whereas the letter
C indicates a bending torsion coupled mode. Following the free vibration analysis, the
flutter analysis is carried out next for all of the three categories of aircraft. Table 4.5 gives
the flutter speed and flutter frequencies of the two sailplanes (S 1 and S 2 ), two light trainer
aircraft (L1 and L2 ) and four transport airliner (T 1 , T 2 , T 3 and T 4 ).
Table 4.4: First five natural frequencies of the baseline aircraft wings.
Table 4.5: Flutter speed and flutter frequency of the baseline aircraft wings.
Next, a detailed parametric study is carried out by varying the bending and torsional stiff-
nesses of each wing for all of the three categories of aircraft (sailplanes, light aircraft
trainers and transport airliners). The stiffness properties are varied between -25% to 25%
in steps of 5% and both free vibration and flutter analyses are carried out. Representative
results are given in Tables 4.6 to 4.37. Tables 4.6 to 4.9 represents the results obtained for
sailplane S 1 , Tables 4.10 to 4.13 represents the results obtained for sailplane S 2 , Tables
4.14 to 4.17 represents the results obtained for light aircraft trainer L1 , Tables 4.18 to 4.21
represents the results obtained for light aircraft trainer L2 , Tables 4.22 to 4.25 represents
the results obtained for transport airliner T 1 , Tables 4.26 to 4.29 represents the results
obtained for transport airliner T 2 , Tables 4.30 to 4.33 represents the results obtained for
transport airliner T 3 and Tables 4.34 to 4.37 represents the results obtained for transport
airliner T 4 .
For sailplane S 1 , the first five natural frequencies due to the variation of EI and GJ are
presented in Tables 4.6 and 4.7 respectively. Clearly it can be seen from the tables, that
when bending stiffness EI is varied bending dominated modes vary significantly but tor-
sional dominated modes vary only slightly. Similarly when torsional stiffness GJ is varied
torsion dominated modes vary significantly while bending dominated modes vary slightly.
This phenomenon as expected is observed for all other aircraft wings.
From Table 4.6, it is clear that for the S 1 sailplane the first, second and fifth modes are
bending modes for the original configuration. Interesting results can be observed for the
third and the fourth modes. The nature of these modes does not change when the bending
stiffness EI is decreased. However, these modes swap from bending to torsion when the
bending stiffness EI is increased as can be seen in the table. The fourth mode remains
torsional when EI is decreased and becomes bending modes when EI is increased. With
respect to the variation of GJ shown in Table 4.7, the nature of the first, second and fifth
modes remains unchanged as bending modes. By contrast, the third mode changes from
bending to torsion when the GJ is reduced while it remains bending as the GJ is increased.
The fourth mode changes from torsion to bending as the GJ is decreased while it remains
torsional as the GJ is increased. It should be noted from the results that the stiffness vari-
ation causes modal interchanges (flip over) between bending and torsion as evident from
Tables 4.6 and 4.7.
4.4. Variations of bending and torsional stiffnesses on aircraft wings 87
Table 4.6: The effects of the variation of EI on the natural frequencies of S 1 wing.
ωi (rad/s)
Variation in EI (%)
ω1 ω2 ω3 ω4 ω5
-25 9.216 (B) 36.91 (B) 94.96 (B) 111.5 (T) 174.4 (B)
-20 9.519 (B) 38.12 (B) 98.07 (B) 111.5 (T) 180.1 (B)
-15 9.812 (B) 39.29 (B) 101.1 (B) 111.5 (T) 185.7 (B)
-10 10.09 (B) 40.43 (B) 104.0 (B) 111.5 (T) 191.1 (B)
-5 10.37 (B) 41.54 (B) 106.9 (B) 111.5 (T) 196.3 (B)
0 10.64 (B) 42.62 (B) 109.6 (B) 111.5 (T) 201.4 (B)
5 10.91 (B) 43.67 (B) 111.5 (T) 112.4 (B) 206.4 (B)
10 11.16 (B) 44.69 (B) 111.5 (T) 114.9 (B) 211.2 (B)
15 11.41 (B) 45.70 (B) 111.5 (T) 117.6 (B) 215.9 (B)
20 11.66 (B) 46.69 (B) 111.5 (T) 120.1 (B) 220.6 (B)
25 11.89 (B) 47.65 (B) 111.5 (T) 122.6 (B) 225.2 (B)
Table 4.7: The effects of the variation of GJ on the natural frequencies of S 1 wing.
ωi (rad/s)
Variation in GJ (%)
ω1 ω2 ω3 ω4 ω5
-25 10.64 (B) 42.62 (B) 96.59 (T) 109.6 (B) 201.4 (B)
-20 10.64 (B) 42.62 (B) 99.76 (T) 109.6 (B) 201.4 (B)
-15 10.64 (B) 42.62 (B) 102.8 (T) 109.6 (B) 201.4 (B)
-10 10.64 (B) 42.62 (B) 105.8 (T) 109.6 (B) 201.4 (B)
-5 10.64 (B) 42.62 (B) 108.7 (T) 109.7 (B) 201.4 (B)
0 10.64 (B) 42.62 (B) 109.6 (B) 111.5 (T) 201.4 (B)
5 10.64 (B) 42.62 (B) 109.6 (B) 114.3 (T) 201.4 (B)
10 10.64 (B) 42.62 (B) 109.6 (B) 116.9 (T) 201.4 (B)
15 10.64 (B) 42.62 (B) 109.6 (B) 119.6 (T) 201.4 (B)
20 10.64 (B) 42.62 (B) 109.6 (B) 122.2 (T) 201.4 (B)
25 10.64 (B) 42.62 (B) 109.6 (B) 124.7 (T) 201.4 (B)
Table 4.8 represents the flutter speed and flutter frequency obtained due to the variation
of EI and Table 4.9 represents the flutter speed and flutter frequency obtained due to the
variation of GJ for sailplane S 1 . It can be seen from Table 4.8 that variation of EI does not
affect flutter speed much, but for Table 4.9 it can be observed that as GJ value increases
88 Chapter 4. High aspect ratio aircraft wings
flutter speed increases. However, when GJ value decreases flutter speed also decreases,
this is to be expected because torsional stiffness GJ generally plays a greater role in flutter
analysis and increasing GJ increases the flutter speed.
Table 4.8: The effects of the variation of EI on the flutter analysis of S 1 wing.
Table 4.9: The effects of the variation of GJ on the flutter analysis of S 1 wing.
For sailplane S 2 , the first five natural frequencies computed by varying EI and GJ are
shown in Tables 4.10 and 4.11 respectively. From the results in Table 4.10, it is clear that
the first, second and third modes are always bending modes. The fourth mode changes
4.4. Variations of bending and torsional stiffnesses on aircraft wings 89
from torsional to a coupled one as the EI is decreased, while this mode remains torsional
as the EI is increased. The fifth mode changes from coupled to a torsional one as the EI is
decreased, while the mode shapes remain as coupled as the EI is increased. With respect
to the variation of GJ in Table 4.11, the modes corresponding to the first, second and
third natural frequencies remain unchanged as bending modes regardless of the variations
within the range. The fourth mode remains torsional as the GJ is decreased while it
changes from torsional to coupled modes as the GJ is increased. The fifth mode remains
as coupled as the GJ is decreased while it changed from coupled to torsional as the GJ is
increased. Here again, the modal interchanges (flip over) between torsional and coupled
modes are prevalent as a result of the GJ variations, see Tables 4.10 and 4.11.
Table 4.10: The effects of the variation of EI on the natural frequencies of S 2 wing.
ωi (rad/s)
Variation in EI (%)
ω1 ω2 ω3 ω4 ω5
-25 11.59 (B) 36.49 (B) 81.05 (B) 145.3 (C) 164.3 (T)
-20 11.97 (B) 37.69 (B) 83.67 (B) 149.9 (C) 164.3 (T)
-15 12.33 (B) 38.84 (B) 86.19 (B) 154.4 (C) 164.4 (T)
-10 12.69 (B) 39.95 (B) 88.65 (B) 158.7 (C) 164.5 (T)
-5 13.04 (B) 41.04 (B) 91.03 (B) 162.6 (C) 164.9 (T)
0 13.38 (B) 42.09 (B) 93.35 (B) 164.2 (T) 167.4 (C)
5 13.71 (B) 43.12 (B) 95.61 (B) 164.4 (T) 171.3 (C)
10 14.03 (B) 44.13 (B) 97.80 (B) 164.5 (T) 175.1 (C)
15 14.34 (B) 45.11 (B) 99.95 (B) 164.6 (T) 178.8 (C)
20 14.65 (B) 46.07 (B) 102.0 (B) 164.7 (T) 182.5 (C)
25 14.95 (B) 47.01 (B) 104.1 (B) 164.7 (T) 186.1 (C)
90 Chapter 4. High aspect ratio aircraft wings
Table 4.11: The effects of the variation of GJ on the natural frequencies of S 2 wing.
ωi (rad/s)
Variation in GJ (%)
ω1 ω2 ω3 ω4 ω5
-25 13.37 (B) 42.03 (B) 93.01 (B) 142.7 (T) 166.2 (C)
-20 13.38 (B) 42.04 (B) 93.17 (B) 147.3 (T) 166.5 (C)
-15 13.38 (B) 42.06 (B) 93.17 (B) 151.8 (T) 166.7 (C)
-10 13.38 (B) 42.07 (B) 93.24 (B) 156.1 (T) 166.9 (C)
-5 13.38 (B) 42.08 (B) 93.30 (B) 160.3 (T) 167.1 (C)
0 13.38 (B) 42.09 (B) 93.35 (B) 164.2 (T) 167.4 (C)
5 13.38 (B) 42.10 (B) 93.39 (B) 166.8 (C) 169.0 (T)
10 13.38 (B) 42.11 (B) 93.44 (B) 167.2 (C) 172.6 (T)
15 13.38 (B) 42.12 (B) 93.48 (B) 167.4 (C) 176.4 (T)
20 13.38 (B) 42.13 (B) 93.51 (B) 167.5 (C) 180.1 (T)
25 13.38 (B) 42.14 (B) 93.54 (B) 167.6 (C) 183.7 (T)
Table 4.12 shows the flutter speed and flutter frequency obtained due to the variation
of EI and Table 4.13 represents the flutter speed and flutter frequency obtained due to
the variation of GJ for sailplane S 2 . It can be seen from Table 4.12 that the variation
of EI does not affect the flutter speed to any appreciable extent, but for Table 4.13 it
can be observed that as GJ value increases flutter speed increases and when GJ value
decreases flutter speed decreases, which is as expected. The same trend was observed for
the sailplane S 1 . Thus from the flutter point of view, the general characteristics of these
two sailplanes are similar.
Table 4.12: The effects of the variation of EI on the flutter analysis of S 2 wing.
Table 4.13: The effects of the variation of GJ on the flutter analysis of S 2 wing.
For light aircraft trainer L1 , the first five natural frequencies due to the variation of EI
and GJ are presented in Tables 4.14 and 4.15 respectively. In Table 4.14, it is clear
that the first, second and third are always bending modes while the fifth mode is always
torsional. However, the fourth mode remains bending as the EI is decreased, but this mode
changes from bending to coupled modes as the EI is increased except for the case when
the increase in EI is 25% for which the mode is torsional. With respect to the variation
of GJ the results shown in Table 4.15 reveal an interesting picture. The first, second and
third natural frequencies of the L1 wing remain unchanged as bending modes while the
fifth mode remains as torsional for most of the cases except when the GJ is reduced by
25% for which the mode becomes coupled as can be seen. By contrast, the fourth mode
remains to be bending dominant as the GJ is increased, but when GJ is reduced the mode
becomes coupled until the reduction is 15%, beyond which the mode becomes torsional.
Table 4.16 gives the flutter speed and flutter frequency obtained due to the variation of
EI and Table 4.17 represents the flutter speed and flutter frequency obtained due to the
variation of GJ for light aircraft trainer L1 . It can be seen from Table 4.16 that as EI value
increases flutter speed increases and when EI value decreases flutter speed decreases thus
flutter speed varies proportional to any change in bending stiffness EI, similarly for Table
4.17 it can be observed that as GJ value increases flutter speed increases and when GJ
value decreases flutter speed decreases, thus flutter speed varies proportional to torsional
stiffness GJ. It should be noted that GJ plays a major role than EI, because the variation
in flutter speed is a lot higher than EI when GJ is used.
92 Chapter 4. High aspect ratio aircraft wings
Table 4.14: The effects of the variation of EI on the natural frequencies of L1 wing.
ωi (rad/s)
Variation in EI (%)
ω1 ω2 ω3 ω4 ω5
-25 41.93 (B) 149.9 (B) 395.8 (B) 753.3 (B) 980.2 (T)
-20 43.31 (B) 154.8 (B) 408.8 (B) 778.0 (B) 980.3 (T)
-15 44.64 (B) 159.6 (B) 421.4 (B) 801.8 (B) 980.3 (T)
-10 45.93 (B) 164.2 (B) 433.6 (B) 825.0 (B) 980.4 (T)
-5 47.19 (B) 168.7 (B) 445.5 (B) 847.5 (B) 980.5 (T)
0 48.42 (B) 173.1 (B) 457.0 (B) 869.4 (B) 980.6 (T)
5 49.61 (B) 177.4 (B) 468.3 (B) 890.8 (C) 980.7 (T)
10 50.78 (B) 181.6 (B) 479.3 (B) 911.5 (C) 980.9 (T)
15 51.92 (B) 185.6 (B) 490.1 (B) 931.7 (C) 981.3 (T)
20 53.04 (B) 189.6 (B) 500.6 (B) 950.9 (C) 982.1 (T)
25 54.13 (B) 193.5 (B) 510.9 (B) 967.8 (T) 984.8 (T)
Table 4.15: The effects of the variation of GJ on the natural frequencies of L1 wing.
ωi (rad/s)
Variation in GJ (%)
ω1 ω2 ω3 ω4 ω5
-25 48.42 (B) 173.1 (B) 456.9 (B) 846.9 (T) 871.7 (C)
-20 48.42 (B) 173.1 (B) 457.0 (B) 865.6 (T) 880.8 (T)
-15 48.42 (B) 173.1 (B) 457.0 (B) 868.5 (C) 905.0 (T)
-10 48.42 (B) 173.1 (B) 457.0 (B) 869.0 (C) 930.6 (T)
-5 48.42 (B) 173.1 (B) 457.0 (B) 869.3 (C) 955.9 (T)
0 48.42 (B) 173.1 (B) 457.0 (B) 869.4 (B) 980.6 (T)
5 48.42 (B) 173.1 (B) 457.0 (B) 869.6 (B) 1005 (T)
10 48.42 (B) 173.1 (B) 457.0 (B) 869.6 (B) 1028 (T)
15 48.42 (B) 173.1 (B) 457.0 (B) 869.7 (B) 1051 (T)
20 48.42 (B) 173.1 (B) 457.1 (B) 869.7 (B) 1074 (T)
25 48.42 (B) 173.1 (B) 457.1 (B) 869.8 (B) 1096 (T)
4.4. Variations of bending and torsional stiffnesses on aircraft wings 93
Table 4.16: The effects of the variation of EI on the flutter analysis of L1 wing.
Table 4.17: The effects of the variation of GJ on the flutter analysis of L1 wing.
For the light aircraft trainer L2 , the first five natural frequencies due to the variation of
EI and GJ are presented in Tables 4.18 and 4.19 respectively. In Table 4.18, it can be
seen that the first and second modes are always bending modes while the third and fifth
modes are torsional except for the 25% increase in EI of the fifth mode for which the
94 Chapter 4. High aspect ratio aircraft wings
mode is coupled. The fourth mode remains by and large a coupled mode except for the
lone case when the increase in EI is 25%, for which the mode is torsional. With respect
to the variation of GJ the natural frequencies are shown in Table 4.19. The first and
second modes remains unchanged as bending and the third mode also remains unchanged
as torsion. However, the fourth mode changes its character from coupled to torsional only
when the GJ is reduced by 20% and 25%. The fifth mode is basically a torsional mode for
majority of the cases of GJ variations except when the GJ is reduced by 25% for which
the mode becomes a coupled one.
Table 4.18: The effects of the variation of EI on the natural frequencies of L2 wing.
ωi (rad/s)
Variation in EI (%)
ω1 ω2 ω3 ω4 ω5
-25 87.55 (B) 351.4 (B) 602.9 (T) 845.6 (C) 1082 (T)
-20 90.41 (B) 362.9 (B) 602.9 (T) 873.2 (C) 1082 (T)
-15 93.18 (B) 374.0 (B) 603.0 (T) 900.1 (C) 1083 (T)
-10 95.86 (B) 384.9 (B) 603.1 (T) 926.1 (C) 1083 (T)
-5 98.47 (B) 395.4 (B) 603.2 (T) 951.4 (C) 1083 (T)
0 101.0 (B) 405.6 (B) 603.2 (T) 975.9 (C) 1083 (T)
5 103.5 (B) 415.6 (B) 603.3 (T) 999.9 (C) 1083 (T)
10 105.9 (B) 425.3 (B) 603.4 (T) 1023 (C) 1084 (T)
15 108.3 (B) 434.8 (B) 603.4 (T) 1046 (C) 1084 (T)
20 110.6 (B) 444.1 (B) 603.5 (T) 1067 (C) 1085 (T)
25 112.8 (B) 453.2 (B) 603.6 (T) 1081 (T) 1094 (C)
Table 4.19: The effects of the variation of GJ on the natural frequencies of L2 wing.
ωi (rad/s)
Variation in GJ (%)
ω1 ω2 ω3 ω4 ω5
-25 100.9 (B) 405.3 (B) 522.8 (T) 937.9 (T) 976.6 (C)
-20 100.9 (B) 405.4 (B) 539.8 (T) 966.5 (T) 978.6 (T)
-15 100.9 (B) 405.5 (B) 556.4 (T) 975.0 (C) 999.8 (T)
-10 101.0 (B) 405.5 (B) 572.4 (T) 975.7 (C) 1028 (T)
-5 101.0 (B) 405.6 (B) 588.0 (T) 975.9 (C) 1056 (T)
0 101.0 (B) 405.6 (B) 603.2 (T) 975.9 (C) 1083 (T)
5 101.0 (B) 405.6 (B) 618.1 (T) 976.1 (C) 1109 (T)
10 101.1 (B) 405.7 (B) 632.5 (T) 976.2 (C) 1136 (T)
15 101.1 (B) 405.7 (B) 646.7 (T) 976.2 (C) 1161 (T)
20 101.1 (B) 405.7 (B) 660.5 (T) 976.3 (C) 1186 (T)
25 101.1 (B) 405.7 (B) 674.1 (T) 976.3 (C) 1210 (T)
4.4. Variations of bending and torsional stiffnesses on aircraft wings 95
Table 4.20 represents the flutter speed and flutter frequency obtained due to the variation
of EI and Table 4.21 represents the flutter speed and flutter frequency obtained due to the
variation of GJ for light aircraft trainer L2 . It can be seen from Table 4.20 that as EI value
increases flutter speed increases and when EI value decreases flutter speed decreases, thus
flutter speed varies proportional to any change in bending stiffness EI, also, it should be
noted that the change in flutter speed is very small. For Table 4.21 it can be observed
that as GJ value increases flutter speed increases and when GJ value decreases flutter
speed decreases, thus flutter speed varies proportional to torsional stiffness GJ. Also the
variation in flutter speed due to torsional stiffness GJ is more evident. The same trend was
observed for the light aircraft trainer L1 . Thus from the flutter point of view, the general
characteristics of these two light aircraft trainers are similar.
Table 4.20: The effects of the variation of EI on the flutter analysis of L2 wing.
Table 4.21: The effects of the variation of GJ on the flutter analysis of L2 wing.
For transport airliner T 1 , the first five natural frequencies computed by varying the EI
and GJ are presented in Tables 4.22 and 4.23, respectively. In Table 4.22, it can be seen
that the first and second modes are always bending modes while the third mode is always
coupled. For the fourth mode, it can be seen that as EI decreases the mode remains as
bending, but it changes from bending to coupled and finally to torsion as EI is increased.
The fifth mode changes from torsion to coupled for the entire range of EI variation. In
Table 4.23, it can be seen that the first and second modes are bending while the third mode
is coupled for all the cases. However, for the fourth mode, it can be seen that as GJ is
decreased the mode changes from bending to coupled and finally towards torsion, but it
remains as a bending mode when GJ is increased. By contrast, the fifth mode changes
from torsion to coupled modes when GJ is increased or decreased except for the lone case
when the reduction in GJ is -25% for which the mode becomes bending. The interesting
phenomena of modal interchanges are again observed in the results of Tables 4.22 and
4.23.
4.4. Variations of bending and torsional stiffnesses on aircraft wings 97
Table 4.22: The effects of the variation of EI on the natural frequencies of T 1 wing.
ωi (rad/s)
Variation in EI (%)
ω1 ω2 ω3 ω4 ω5
-25 9.983 (B) 28.69 (B) 45.37 (C) 76.27 (B) 97.64 (C)
-20 10.31 (B) 29.63 (B) 45.38 (C) 78.74 (B) 97.66 (C)
-15 10.63 (B) 30.53 (B) 45.38 (C) 81.13 (B) 97.68 (C)
-10 10.93 (B) 31.41 (B) 45.39 (C) 83.45 (B) 97.71 (C)
-5 11.23 (B) 32.26 (B) 45.39 (C) 85.68 (B) 97.73 (C)
0 11.52 (B) 33.09 (B) 45.40 (C) 87.86 (B) 97.75 (T)
5 11.81 (B) 33.89 (B) 45.40 (C) 89.98 (B) 97.77 (C)
10 12.09 (B) 34.68 (B) 45.41 (C) 92.04 (B) 97.79 (C)
15 12.36 (B) 35.45 (B) 45.42 (C) 94.04 (C) 97.82 (C)
20 12.62 (B) 36.19 (B) 45.43 (C) 95.97 (C) 97.86 (C)
25 12.88 (B) 36.93 (B) 45.44 (C) 97.61 (T) 98.14 (C)
Table 4.23: The effects of the variation of GJ on the natural frequencies of T 1 wing.
ωi (rad/s)
Variation in GJ (%)
ω1 ω2 ω3 ω4 ω5
-25 11.52 (B) 33.01 (B) 39.36 (C) 84.74 (T) 87.46 (B)
-20 11.52 (B) 33.03 (B) 40.64 (C) 87.31 (T) 87.78 (C)
-15 11.52 (B) 33.05 (B) 41.88 (C) 87.65 (C) 90.20 (C)
-10 11.53 (B) 33.06 (B) 43.08 (C) 87.74 (B) 92.78 (C)
-5 11.52 (B) 33.08 (B) 44.26 (C) 87.81 (B) 95.29 (C)
0 11.52 (B) 33.09 (B) 45.40 (C) 87.86 (B) 97.75 (T)
5 11.53 (B) 33.09 (B) 46.51 (C) 87.90 (B) 100.1 (C)
10 11.53 (B) 33.10 (B) 47.60 (C) 87.94 (B) 102.5 (C)
15 11.53 (B) 33.11 (B) 48.67 (C) 87.98 (B) 104.8 (C)
20 11.53 (B) 33.12 (B) 49.71 (C) 88.01 (B) 106.9 (C)
25 11.53 (B) 33.12 (B) 50.73 (C) 88.03 (B) 109.2 (C)
Table 4.24 represents the flutter speed and flutter frequency obtained due to the variation
of EI and Table 4.25 represents the flutter speed and flutter frequency obtained due to
the variation of GJ for transport airliner T 1 . It can be seen from Table 4.24 that as EI
value varies flutter speed fluctuates, the flutter speed tends to increase for minor change
in EI values for both cases of either increase or decrease but the flutter speed decreases
98 Chapter 4. High aspect ratio aircraft wings
for major variation in EI values. For Table 4.25 it can be observed that as GJ changes,
the flutter speed tends to increase for minor change in GJ values but the flutter speed
decreases for major variation in GJ values.
Table 4.24: The effects of the variation of EI on the flutter analysis of T 1 wing.
Table 4.25: The effects of the variation of GJ on the flutter analysis of T 1 wing.
For transport airliner T 2 , the first five natural frequencies due to the variation of EI and
GJ are presented in Tables 4.26 and 4.27 respectively. In Table 4.26, It can be seen that
the first, second and third modes are always bending while the fourth and fifth modes
4.4. Variations of bending and torsional stiffnesses on aircraft wings 99
are always coupled modes. Surprisingly, there is no change in the characterisation of the
modes due to the variation of EI for the entire range. Table 4.27, which shows the effects
of the variation of GJ on the first five natural frequencies reveal similar characteristics in
that the nature of the modes remains unchanged due to GJ variations.
Table 4.26: The effects of the variation of EI on the natural frequencies of T 2 wing.
ωi (rad/s)
Variation in EI (%)
ω1 ω2 ω3 ω4 ω5
-25 17.09 (B) 48.15 (B) 87.41 (B) 119.5 (C) 176.7 (C)
-20 17.64 (B) 49.64 (B) 90.14 (B) 119.8 (C) 181.4 (C)
-15 18.18 (B) 51.13 (B) 92.81 (B) 120.1 (C) 185.9 (C)
-10 18.71 (B) 52.56 (B) 95.37 (B) 120.4 (C) 190.1 (C)
-5 19.21 (B) 53.95 (B) 97.85 (B) 120.7 (C) 194.1 (C)
0 19.71 (B) 55.29 (B) 100.3 (B) 120.9 (C) 197.7 (C)
5 20.19 (B) 56.60 (B) 102.6 (B) 121.1 (C) 201.3 (C)
10 20.66 (B) 57.86 (B) 104.8 (B) 121.4 (C) 204.5 (C)
15 21.13 (B) 59.11 (B) 107.0 (B) 121.6 (C) 207.6 (C)
20 21.57 (B) 60.31 (B) 109.2 (B) 121.9 (C) 210.3 (C)
25 22.02 (B) 61.51 (B) 111.3 (B) 122.1 (C) 212.9 (C)
Table 4.27: The effects of the variation of GJ on the natural frequencies of T 2 wing.
ωi (rad/s)
Variation in GJ (%)
ω1 ω2 ω3 ω4 ω5
-25 19.68 (B) 54.87 (B) 99.14 (B) 106.1 (C) 187.8 (C)
-20 19.68 (B) 54.98 (B) 99.47 (B) 109.2 (C) 190.5 (C)
-15 19.69 (B) 55.07 (B) 99.72 (B) 112.3 (C) 192.7 (C)
-10 19.70 (B) 55.15 (B) 99.93 (B) 115.2 (C) 194.7 (C)
-5 19.70 (B) 55.23 (B) 100.1 (B) 118.1 (C) 196.3 (C)
0 19.71 (B) 55.29 (B) 100.3 (B) 120.9 (C) 197.7 (C)
5 19.71 (B) 55.34 (B) 100.4 (B) 123.6 (C) 199.0 (C)
10 19.71 (B) 55.39 (B) 100.5 (B) 126.3 (C) 200.1 (C)
15 19.72 (B) 55.44 (B) 100.6 (B) 128.9 (C) 201.1 (C)
20 19.72 (B) 55.48 (B) 100.7 (B) 131.5 (C) 201.9 (C)
25 19.72 (B) 55.52 (B) 100.8 (B) 134.1 (C) 202.8 (C)
100 Chapter 4. High aspect ratio aircraft wings
Table 4.28 represents the flutter speed and flutter frequency obtained due to the variation
of EI and Table 4.29 represents the flutter speed and flutter frequency obtained due to the
variation of GJ for transport airliner T 2 . It can be seen from Table 4.28 that as EI value
increases flutter speed decreases gradually and when EI value decreases flutter speed in-
creases gradually, it can be said that the flutter speed varies inversely proportional to any
change in bending stiffness EI, also, it should be noted that the change in flutter speed
is very small. For Table 4.29 it can be observed that as GJ value increases flutter speed
increases and when GJ value decreases flutter speed decreases, thus flutter speed varies
proportional to torsional stiffness GJ. It should be noted that T 2 consists of single engine
on its wing.
Table 4.28: The effects of the variation of EI on the flutter analysis of T 2 wing.
Table 4.29: The effects of the variation of GJ on the flutter analysis of T 2 wing.
For the transport airliner T 3 , the first five natural frequencies due to the variation of EI and
GJ are shown in Tables 4.30 and 4.31 respectively. In Table 4.30, it can be seen that the
first and second modes are always bending modes. However, for the third mode when EI
is decreased, the mode always remains bending whereas when EI is increased, the modes
remain bending until the increase is 10%, beyond which the modes become coupled. The
fourth mode remains torsional for decrease in the EI values but it becomes a coupled
mode when EI is increased on or above the 10% value. Now the fifth mode is always the
torsional mode for any increase of the EI value but it becomes a coupled mode when EI
is reduced by 10% or more. The results for the GJ variation are shown in Table 4.31. It
can be observed that the first and second modes are always bending modes. For the third
mode, if the GJ is decreased the mode changes from bending to coupled and then finally
to torsional for -15%, -20% and -25% variations, respectively. The character of this mode
remains unchanged when GJ is increased. For the fourth mode when GJ decreases, the
mode becomes coupled from 10% reduction and beyond. However, when GJ increases
the mode remains always torsional. With regard to the fifth mode, any decrease in GJ
does not alter the basic torsional nature of the mode but when GJ is increased beyond 5%
the mode becomes coupled.
102 Chapter 4. High aspect ratio aircraft wings
Table 4.30: The effects of the variation of EI on the natural frequencies of T 3 wing.
ωi (rad/s)
Variation in EI (%)
ω1 ω2 ω3 ω4 ω5
-25 10.39 (B) 30.03 (B) 58.77 (B) 72.51 (T) 107.9 (C)
-20 10.73 (B) 31.00 (B) 60.63 (B) 72.55 (T) 110.1 (C)
-15 11.06 (B) 31.94 (B) 62.43 (B) 72.59 (T) 110.9 (C)
-10 11.38 (B) 32.85 (B) 64.17 (B) 72.63 (T) 111.3 (C)
-5 11.69 (B) 33.73 (B) 65.85 (B) 72.68 (T) 111.5 (T)
0 11.99 (B) 34.59 (B) 67.47 (B) 72.74 (T) 111.6 (T)
5 12.28 (B) 35.42 (B) 69.02 (B) 72.82 (T) 111.8 (T)
10 12.57 (B) 36.23 (B) 70.46 (B) 72.95 (C) 111.9 (T)
15 12.85 (B) 37.03 (B) 71.68 (C) 73.28 (C) 112.0 (T)
20 13.13 (B) 37.79 (B) 72.33 (C) 74.13 (C) 112.1 (T)
25 13.39 (B) 38.55 (B) 72.56 (C) 75.36 (C) 112.2 (T)
Table 4.31: The effects of the variation of GJ on the natural frequencies of T 3 wing.
ωi (rad/s)
Variation in GJ (%)
ω1 ω2 ω3 ω4 ω5
-25 11.98 (B) 34.45 (B) 62.98 (T) 67.16 (C) 97.38 (T)
-20 11.98 (B) 34.48 (B) 64.90 (C) 67.40 (C) 100.4 (T)
-15 11.98 (B) 34.52 (B) 66.47 (C) 67.91 (C) 103.3 (T)
-10 11.99 (B) 34.54 (B) 67.13 (B) 69.26 (T) 106.2 (T)
-5 11.99 (B) 34.57 (B) 67.35 (B) 70.98 (T) 108.9 (T)
0 11.99 (B) 34.59 (B) 67.47 (B) 72.74 (T) 111.6 (T)
5 11.99 (B) 34.60 (B) 67.56 (B) 74.48 (T) 114.3 (T)
10 11.99 (B) 34.62 (B) 67.63 (B) 74.19 (T) 116.8 (C)
15 11.99 (B) 34.63 (B) 67.69 (B) 77.86 (T) 119.1 (C)
20 11.99 (B) 34.65 (B) 67.74 (B) 79.50 (T) 121.3 (C)
25 11.99 (B) 34.66 (B) 67.79 (B) 81.11 (T) 123.0 (C)
Table 4.32 represents the flutter speed and flutter frequency obtained due to the variation
of EI and Table 4.33 represents the flutter speed and flutter frequency obtained due to the
variation of GJ for transport airliner T 3 . It can be seen from Table 4.32 that as EI value
increases flutter speed decreases and when EI value decreases flutter speed increases, it
can be said that the flutter speed varies inversely proportional to any change in bending
stiffness EI, also, it should be noted that the change in flutter speed is very small. For Table
4.33 it can be observed that as GJ value increases flutter speed increases and when GJ
4.4. Variations of bending and torsional stiffnesses on aircraft wings 103
value decreases flutter speed decreases, thus flutter speed varies proportional to torsional
stiffness GJ. It should be noted that T 3 consists of single engine on its wing. This is
similar to the characteristics observed in transport airliner T 2 .
Table 4.32: The effects of the variation of EI on the flutter analysis of T 3 wing.
Table 4.33: The effects of the variation of GJ on the flutter analysis of T 3 wing.
For transport airliner T 4 , the first five natural frequencies due to the variation of EI and GJ
are presented in Tables 4.34 and 4.35 respectively. In Table 4.34, it can be seen from the
mode shapes analysed that the first, second and fourth modes are always bending while
third and fifth modes are always torsional. There is virtually no changes in the character
of the modes due to the variation of EI. In Table 4.35, similar observations are made when
GJ is varied. Clearly, the first, second and fourth modes are always bending while third
and fifth modes are always torsional.
Table 4.34: The effects of the variation of EI on the natural frequencies of T 4 wing.
ωi (rad/s)
Variation in EI (%)
ω1 ω2 ω3 ω4 ω5
-25 7.784 (B) 22.90 (B) 45.26 (T) 62.41 (B) 94.08 (T)
-20 8.039 (B) 23.65 (B) 45.26 (T) 64.46 (B) 94.09 (T)
-15 8.286 (B) 24.38 (B) 45.26 (T) 66.44 (B) 94.09 (T)
-10 8.526 (B) 25.09 (B) 45.26 (T) 68.37 (B) 94.09 (T)
-5 8.760 (B) 25.78 (B) 45.26 (T) 70.24 (B) 94.09 (T)
0 8.988 (B) 26.45 (B) 45.26 (T) 72.06 (B) 94.09 (T)
5 9.210 (B) 27.09 (B) 45.26 (T) 73.84 (B) 94.09 (T)
10 9.426 (B) 27.74 (B) 45.26 (T) 75.58 (B) 94.09 (T)
15 9.638 (B) 28.36 (B) 45.26 (T) 77.28 (B) 94.09 (T)
20 9.846 (B) 28.97 (B) 45.26 (T) 78.94 (B) 94.09 (T)
25 10.05 (B) 29.57 (B) 45.26 (T) 80.57 (B) 94.09 (T)
Table 4.35: The effects of the variation of GJ on the natural frequencies of T 4 wing.
ωi (rad/s)
Variation in GJ (%)
ω1 ω2 ω3 ω4 ω5
-25 8.988 (B) 26.45 (B) 39.19 (T) 72.06 (B) 81.49 (T)
-20 8.988 (B) 26.45 (B) 40.48 (T) 72.06 (B) 84.16 (T)
-15 8.988 (B) 26.45 (B) 41.73 (T) 72.06 (B) 86.75 (T)
-10 8.988 (B) 26.45 (B) 42.94 (T) 72.06 (B) 89.26 (T)
-5 8.988 (B) 26.45 (B) 44.12 (T) 72.06 (B) 91.71 (T)
0 8.988 (B) 26.45 (B) 45.26 (T) 72.06 (B) 94.09 (T)
5 8.988 (B) 26.45 (B) 46.38 (T) 72.06 (B) 96.41 (T)
10 8.988 (B) 26.45 (B) 47.47 (T) 72.06 (B) 98.68 (T)
15 8.988 (B) 26.45 (B) 48.54 (T) 72.07 (B) 100.9 (T)
20 8.988 (B) 26.45 (B) 49.58 (T) 72.07 (B) 103.1 (T)
25 8.988 (B) 26.45 (B) 50.60 (T) 72.07 (B) 105.2 (T)
4.4. Variations of bending and torsional stiffnesses on aircraft wings 105
Table 4.36 represents the flutter speed and flutter frequency obtained due to the variation
of EI and Table 4.37 represents the flutter speed and flutter frequency obtained due to the
variation of GJ for transport airliner T 4 . It can be seen from Table 4.36 that as EI value
increases flutter speed decreases and when EI value decreases flutter speed increases and
it can be said that the flutter speed varies inversely proportional to bending stiffness EI.
For Table 4.37 it can be observed that as GJ value increases flutter speed increases and
when GJ value decreases flutter speed decreases, thus flutter speed varies proportional to
torsional stiffness GJ. This is similar with respect to the characteristics observed in other
transport airliners.
Table 4.36: The effects of the variation of EI on the flutter analysis of T 4 wing.
Table 4.37: The effects of the variation of GJ on the flutter analysis of T 4 wing.
4.5 Effects due to the variation of engine mass and its lo-
cation
Next, a parametric study is carried out for transport airliner wings by varying the engine
mass and its location. Both sailplanes and the light aircraft trainer wings have been ex-
cluded because they do not have any engines on them. Thus, only the transport airliner
wings are considered in the parametric study based on the changes of engine masses and
their locations. It should be noted that among the four transport airliners, T 1 and T 4 have
two engines on each of their wings whereas T 2 and T 3 , each has a single engine mounted
on each of their wings. A variation between -25% to 25% in steps of 5% is allowed and for
the engine locations, some realistic distances from the wing root and the current locations
of the engines are considered. The first five natural frequencies with the identifications of
the modes and the flutter analysis are presented in Tables 4.38 to 4.45 which are discussed
next. The results corresponding to the original mass and its location of the baseline wing
are shown in bold.
For the transport airliner T 1 , the first five natural frequencies due to the variation of engine
masses and its locations are presented in Tables 4.38 and 4.39 respectively. It should be
noted that the transport airliner T 1 has two engines on its wing and the masses of each of
the two engines are equally increased by the same proportion as indicated in Table 4.38.
In Table 4.38, it can also be seen that the first, second and fourth modes are always bend-
ing while third mode is always coupled and fifth mode is always torsional. There is no
change in their modal characteristics due to the variation of engine mass. Flutter speed
increases as engine mass is increased and the flutter speed decreases as engine mass is
decreased, thus flutter speed varies directly proportional to change in engine mass. Table
4.39 shows that the first, second and fourth modes are always bending while the third
mode is coupled and the fifth mode is torsional. Here also, there is no change in their
modal characteristics due to the variation of engine location. It can be seen as the two
engines are brought closer (9.53, 6.93), lowest flutter speed occurred and when the two
engines are separated apart (12.07, 4.32), highest flutter speed occurred. Overall it can be
noted that the right balance in engine mass and location provides the best flutter speed,
this can be observed in other aircraft wings analysed too.
4.5. Effects due to the variation of engine mass and its location 107
Table 4.38: The effects of the variation of engine mass on the natural frequencies, flutter
speeds and flutter frequencies of T 1 .
Table 4.39: The effects of the variation of engine location on the natural frequencies,
flutter speeds and flutter frequencies of T 1 .
Location Location
Flutter Flutter
of engine of engine
ωi (rad/s) speed, frequency,
1 from 2 from
Vs ω
wing root wing root
(m/s) (rad/s)
(m) (m) ω1 ω2 ω3 ω4 ω5
12.07 5.59 10.73 (B) 34.40 (B) 42.45 (C) 79.79 (B) 92.83 (T) 254 27.63
12.07 4.32 10.80 (B) 35.10 (B) 42.96 (C) 86.25 (B) 95.67 (T) 256 27.92
10.80 6.93 11.36 (B) 32.13 (B) 44.05 (C) 86.87 (B) 97.10 (T) 252 28.84
10.80 5.59 11.52 (B) 33.09 (B) 45.40 (C) 87.86 (B) 97.75 (T) 249 28.77
10.80 4.32 11.61 (B) 33.78 (B) 46.17 (C) 93.66 (B) 98.95 (T) 246 29.49
9.53 6.93 12.00 (B) 31.79 (B) 45.62 (C) 91.88 (B) 96.55 (T) 234 26.47
9.53 5.59 12.18 (B) 32.85 (B) 47.28 (C) 91.66 (B) 96.91 (T) 252 27.81
108 Chapter 4. High aspect ratio aircraft wings
For the transport airliner T 2 , the first five natural frequencies due to the variation of en-
gine mass and its location are presented in Tables 4.40 and 4.41 respectively. It should
be noted that the transport airliner T 2 has single engine on its wing. In Table 4.40, it can
be seen that the first, second and third modes are always bending while the fifth mode
is always coupled. Also, the fourth mode is coupled but as its engine mass is reduced
to 25% it change to torsional. Flutter speed increases as engine mass is increased and
the flutter speed decreases as engine mass is decreased, thus flutter speed varies directly
proportional to change in engine mass. In Table 4.41, It can be seen that the first and
second modes are always bending while the third mode changes from bending to coupled
as the engine is moved towards the wing tip but remains bending as the engine is moved
towards wing root. The fourth mode changes from coupled to torsional as the engine is
moved towards the wing root. Also the fourth mode changes from coupled to bending
as the engine is move towards wing tip. The fifth mode remains coupled as the engine is
moved towards the root but changed to bending when it is moved towards wing tip. It can
be seen as the engine is brought closer to the wing root, the flutter speed is increased and
as the engine is moved away from the root, the flutter speed is decreased.
Table 4.40: The effects of the variation of engine mass on the natural frequencies, flutter
speeds and flutter frequencies of T 2 .
Table 4.41: The effects of the variation of engine location on the natural frequencies,
flutter speeds and flutter frequencies of T 2 .
For the transport airliner T 3 , the first five natural frequencies due to the variation of engine
mass and its location are presented in Tables 4.42 and 4.43 respectively. It should be noted
that the transport airliner T 3 has single engine on its wing. Clearly it can be seen from
the Table 4.42, all the natural frequencies are virtually unaltered due to the variation of
engine mass. It can also be seen that the first, second and third modes are always bending
while fourth and fifth modes are always torsional. It can be seen any change in engine
mass reduces the flutter speed. In Table 4.43, it can be seen that the first and the second
modes are always bending and the fifth mode is always torsional. The third mode changes
from bending to coupled and then to torsional as the engine is moved towards the wing
root and it remains bending when it is moved towards the wing tip. The fourth modes
change from torsional to coupled and then to bending as the engine is moved towards the
wing root and it remains torsional when it is moved towards the wing tip. It can be seen
as the engine is brought closer to the wing root, the flutter speed is increased and as the
engine is moved away from the root, the flutter speed is decreased.
110 Chapter 4. High aspect ratio aircraft wings
Table 4.42: The effects of the variation of engine mass on the natural frequencies, flutter
speeds and flutter frequencies of T 3 .
Table 4.43: The effects of the variation of engine location on the natural frequencies,
flutter speeds and flutter frequencies of T 3 .
For the transport airliner T 4 , the first five natural frequencies due to the variation of engine
masses and its locations are presented in Tables 4.44 and 4.45 respectively. It should be
noted that the transport airliner T 4 has two engines on its wing. Note that, the masses of
each of the two engines are equally increased by the same proportion as indicated in Table
4.44. Clearly it can be seen from the Table 4.44, all the natural frequencies variation is
quite small due to the variation of the engine masses from -25% to +25%. It can also be
seen from the mode shapes analysed that the first, second and fourth modes are always
bending while third and fifth modes are torsional. There is no change in their modes due
4.5. Effects due to the variation of engine mass and its location 111
to the variation of engine mass. It can be seen that as the engine mass is decreased, flutter
speed increases and as the engine mass is increased, the flutter speed decreases.
For Table 4.45, It can be seen that the first, second and fourth modes are always bending
while third and fifth modes are always torsional. There is no appreciable change in the
modal characteristics due to the variation of the engine positions. It can be seen as the
two engines are brought closer, low flutter speed occurs not the lowest though and when
the two engines are separated apart, high value for flutter speed occurs and it can be
noticed that the lowest or highest flutter speed does not occur solely on the basis of engine
location.
Table 4.44: The effects of the variation of engine mass on the natural frequencies, flutter
speeds and flutter frequencies of T 4 .
Table 4.45: The effects of the variation of engine location on the natural frequencies,
flutter speeds and flutter frequencies of T 4 .
Location Location
Flutter Flutter
of engine of engine
ωi (rad/s) speed, frequency,
1 from 2 from
Vs ω
wing root wing root
(m/s) (rad/s)
(m) (m) ω1 ω2 ω3 ω4 ω5
17.6 5.6 8.811 (B) 26.69 (B) 45.07 (T) 70.65 (B) 94.31 (T) 479 29.14
17.6 4.6 8.814 (B) 26.73 (B) 45.08 (T) 72.67 (B) 94.35 (T) 477 29.20
17.1 7.4 8.972 (B) 26.30 (B) 45.22 (T) 66.62 (B) 94.00 (T) 390 46.51
17.1 5.6 8.988 (B) 26.45 (B) 45.26 (T) 72.06 (B) 94.09 (T) 386 46.50
17.1 4.6 8.991 (B) 26.48 (B) 45.28 (T) 74.15 (B) 94.13 (T) 386 46.50
15.8 7.4 9.400 (B) 25.92 (B) 45.82 (T) 68.11 (B) 93.73 (T) 433 43.00
15.8 4.6 9.418 (B) 26.08 (B) 45.87 (T) 72.88 (B) 93.82 (T) 428 42.99
It should be noted that unlike the cases for T 2 and T 3 wings which carry single engine,
the analysis of the aircraft with two engines on each of its wing reveals that there would
be no appreciable change in their modal characteristics as a consequence of reasonably
realistic variations of engine masses and their locations.
4.6 Conclusions
Using the dynamic stiffness method for the theoretical development and the Wittrick-
Williams algorithm as the solution technique, the modal behaviour of a wide range of
aircraft wings has been investigated which includes two sailplane, two light aircraft trainer
and four transport airliner wings. Natural frequencies, mode shapes and flutter analysis
for the wings are presented and the results are discussed in detail. The bending and
torsional stiffnesses of each wing are varied between +25% and -25% in steps of 5%
and their subsequent effects on the natural frequencies, mode shapes and flutter speed
are investigated which showed significant changes in results and some interesting trends.
As an overview of the results observed in the analysis it can be said that for sailplanes,
the variation of EI does not affect the flutter speed to any appreciable extent while the
flutter speed varies proportional to the change in GJ. For light aircraft trainers, the flutter
speed varies proportional to both EI and GJ. For transport airliner, the flutter speed varies
inversely proportional to change in EI and directly proportional to change in GJ. Overall
GJ variation affects flutter speed more than EI variation. An interesting feature of this
particular investigation reveals that the modal interchanges or flip-overs between bending,
4.6. Conclusions 113
torsional and/or coupled modes can occur as a result of the variations which can have
profound influence in aeroelastic studies. Subsequent parametric studies on the effects
of variation of engine mass and its location that are applicable to the transport airliners
showed that the natural frequencies, mode shapes and flutter speed can be changed to
some appreciable extent to make provision for the avoidance of aeroelastic problems. It
can be concluded that the right balance in engine mass and location provides the best
flutter speed.
Chapter 5
that of a transport airliner for which some pertinent details are given in Table 5.1.
5.2 Results for the free vibration and flutter analysis us-
ing model Type I
As explained above, the model Type I considers the whole aircraft configuration by lump-
ing the fuselage, tailplane, fin and rudder masses and their inertias at the nodal point of
the wing and fuselage centreline intersection, as shown in Figure 5.1, without represent-
ing the fuselage, tailplane, fin and rudder by individual bending-torsion coupled beam
elements unlike the case for the wing. The initial values of the data for the lumped mass
and pitching moment of inertia for the rest of the aircraft other than the wing, for the orig-
inal aircraft were taken to be 15,000 kg and 70,000 kgm2 , respectively. These values were
halved because only one symmetric half of the aircraft was considered for the analysis.
The bending (EI) and torsional (GJ) stiffnesses along the span for this transport aircraft
are illustrated in Figure 5.2.
The first ten natural frequencies and mode shapes of the aircraft were computed. As
the aircraft was completely free and exhibiting natural vibrations in symmetric motion,
the first two natural frequencies turned out to be zero (from a numerical standpoint, the
value was exceedingly small and close to zero) corresponding to the rigid body modes in
heave and pitch. The first five elastic modes ω3 to ω7 are illustrated in Figure 5.3, where
116 Chapter 5. Whole aircraft analysis
Figure 5.1: Representation of lumped masses and inertias of fuselage, tailplane, fin and
rudder for half of an aircraft.
Transport airliner
Cantilevered aircraft wing
Whole aircraft configuration
configuration
Flutter Speed (m/s) 203.3 249.0
Flutter Frequency (rad/s) 51.63 28.77
the bending and torsional displacements are respectively shown by solid black lines and
broken red lines. The first three elastic modes are basically bending dominated modes,
whereas the fourth elastic one is a coupled mode. By contrast the fifth elastic mode is
more or less a torsional mode. In the subsequent flutter analysis, these elastic modes as
well as the two rigid body modes in heave and pitch were included. The results of the
analysis for the whole aircraft configuration showing flutter speed and flutter frequency
are shown in Table 5.2 alongside the results computed by using the cantilever model.
5.3. Effects of fuselage mass and inertia on transport airliner 117
Variation Moment
Flutter Flutter
of of ωi (rad/s)
speed frequency
fuselage inertia
(m/s) (rad/s)
mass (%) (%) ω3 ω4 ω5 ω6 ω7 ω8 ω9 ω10
-25 -25 14.02 23.15 43.40 52.35 99.81 106.6 188.4 209.9 202.8 52.13
-20 -25 13.94 23.10 43.39 52.13 99.81 106.4 188.2 209.9 203.6 52.13
-15 -25 13.87 23.07 43.38 51.93 99.81 106.3 188.0 209.9 205.3 52.13
-10 -25 13.80 23.03 43.38 51.75 99.80 106.1 187.9 209.9 204.9 51.63
-5 -25 13.73 23.00 43.37 51.58 99.80 106.0 187.7 209.9 203.5 51.63
0 -25 13.67 22.97 43.36 51.43 99.80 105.9 187.6 209.9 202.8 51.63
5 -25 13.62 22.94 43.35 51.29 99.80 105.7 187.4 209.9 204.2 51.63
10 -25 13.56 22.91 43.35 51.16 99.80 105.6 187.3 209.9 205.5 51.63
15 -25 13.51 22.89 43.34 51.04 99.80 105.6 187.2 209.9 206.5 51.13
20 -25 13.46 22.87 43.33 50.93 99.80 105.5 187.1 209.9 205.7 51.13
Moment
Flutter Flutter
Fuselage of ωi (rad/s)
speed frequency
mass inertia
(m/s) (rad/s)
(%) (%) ω3 ω4 ω5 ω6 ω7 ω8 ω9 ω10
-25 25 14.01 22.74 43.19 52.29 99.81 106.6 188.4 209.9 202.6 52.13
-20 25 13.93 22.69 43.18 52.07 99.81 106.4 188.2 209.9 204.4 52.13
-15 25 13.85 22.65 43.17 51.87 99.81 106.3 188.0 209.9 204.9 51.63
-10 25 13.78 22.62 43.17 51.68 99.80 106.1 187.9 209.9 203.3 51.63
-5 25 13.72 22.59 43.16 51.52 99.80 106.0 187.7 209.9 202.1 51.63
0 25 13.66 22.56 43.15 51.36 99.80 105.9 187.6 209.9 203.6 51.63
5 25 13.60 22.53 43.15 51.22 99.80 105.7 187.4 209.9 205.0 51.63
10 25 13.55 22.50 43.14 51.09 99.80 105.6 187.3 209.9 205.9 51.13
15 25 13.49 22.48 43.13 50.97 99.80 105.5 187.2 209.9 204.9 51.13
20 25 13.45 22.46 43.13 50.86 99.80 105.5 187.1 209.9 204.0 51.13
25 25 13.40 22.43 43.12 50.76 99.80 105.4 187.0 209.9 203.1 51.13
119
120 Chapter 5. Whole aircraft analysis
5.4 Results for the free vibration and flutter analysis us-
ing model Type II
The Type II model for the whole aircraft configuration includes bending-torsion coupled
beams elements for all the components of the aircraft namely, fuselage, tail plane, fin and
rudder. The model also includes the effect of the lumped mass or inertia such as that of an
engine. A schematic diagram of such a model is shown in Figure 5.4 for one symmetric
half of the aircraft.
It should be noted that for the symmetric case, all out of plane displacements and ro-
tations (∆y , Θ x , Θz ) on the enforced plane of symmetry will be zero whereas for the
anti-symmetric case, all in plane displacements and rotations (∆ x , ∆z , Θy ) on the enforced
plane of symmetry will be zero. The stick model configuration having 16 elements for the
wing, 18 elements for the fuselage, 8 elements for the tailplane and 6 elements for the fin
and rudder are used in the analysis for both the symmetric and antisymmetric cases. When
computing the flutter speed and flutter frequency only the wing modes which included the
fundamental bending and fundamental torsion were considered, but the rigid body modes
were implemented in the analysis for both symmetric and antisymmetric cases. In addi-
tion to the rigid body modes of the whole aircraft, five elastic modes of the wings were
included in the flutter analysis. The results for the flutter speed and flutter frequencies for
both symmetric and antisymmetric cases are shown in Table 5.5.
Figure 5.4: Bending-Torsion coupled beam idealisation for one symmetric half of an
aircraft.
Clearly, the antisymmetric flutter speed for the whole aircraft configuration is significantly
lower than the symmetric case as shown in Table 5.5 which is not unusual. It should be
noted that the unsteady aerodynamic forces arising from the tailplane for the symmetric
5.5. Conclusions 121
Table 5.5: Flutter speed and flutter frequency for the whole aircraft configuration using
Type II.
Transport airliner
Symmetric whole Antisymmetric whole Cantilevered aircraft
aircraft configuration aircraft configuration wing configuration
Flutter speed (m/s) 248.8 174.5 249.0
Flutter frequency (rad/s) 48.12 32.45 28.77
case and fin and rudder for the antisymmetric case were not taken into account in the
flutter analysis.
For the symmetric flutter analysis of whole aircraft configuration, the rigid body modes
in heave and pitch have been included along with five elastic modes whereas for the
antisymmetric flutter analysis for the whole aircraft configuration, the rigid body modes
in rolling motion is accounted for during the flutter analysis. Clearly, when a wing only
analysis is carried out with cantilever end condition, the symmetric and antisymmetric
cases cannot be properly covered, let alone the non-inclusion of the rigid body modes due
to the fixity of the built-in end of the aircraft wing.
5.5 Conclusions
Using the dynamic stiffness method for the theoretical development and the Wittrick-
Williams algorithm as the solution technique, the free vibration and flutter analysis for
whole aircraft configuration has been carried out. Two model types have been considered
for whole aircraft analysis. In the first type of analysis masses and inertia are lumped at
the wing and fuselage centreline intersection and for the second type of analysis both sym-
metric and antisymmetric motions are considered for whole aircraft modelled as bending-
torsion coupled beams. Also a parametric study is carried out by varying lumped masses
and inertias between +25% and -25% in steps of 5% and their subsequent effects on the
natural frequencies, mode shapes and flutter speed are investigated for the first type of air-
craft model. It can be seen from the results that their contribution to natural frequencies,
flutter speeds and flutter frequencies are very minimal.
Sect Section B
An overview and layout of Section B
Free vibration analysis is of great importance in the design of aeronautical, civil, auto-
mobile and marine engineering structures, amongst the others. These structures are often
idealised by beam elements such as the case with high aspect ratio aircraft wings. The
traditional finite element method (FEM) is generally used when carrying out the static and
dynamic analysis of beam and other structures. Clearly the order of the mass and stiffness
matrices in the FEM decides the number of natural frequencies that can be meaningfully
computed. The higher order natural frequencies and mode shapes, will of course, be con-
siderably less accurate. It should be noted that there is a powerful alternative to FEM as
well the classical method, which has no restriction on higher order natural frequency com-
putation and yet it retains the exactness of results. The alternative is that of the dynamic
stiffness method (DSM) which is elegant and versatile and can be used to analyse the
free vibration behaviour of complex structures. The DSM is different, but in many ways
similar to the FEM in that it has analogous procedure for assembling structural properties
of individual structural elements. However, a major difference exists between the DSM
and the FEM which is that the former is unaffected by the number of elements used in
the analysis and always gives exact results whereas the latter is mesh dependent and the
accuracy of results depends on the number of elements used in the analysis. For instance,
one single structural element can be used in the DSM to compute any number of natural
frequencies without any loss of accuracy which, of course, is impossible in the FEM. The
uncompromising accuracy of the DSM stems from the fact that the frequency dependent
shape function used to derive the element dynamic stiffness matrix of a structural ele-
ment comes from the exact solution of the governing differential equation of motion of
the element undergoing free natural vibration. The overall frequency dependent dynamic
stiffness matrix of the final structure is obtained by assembling the individual dynamic
stiffness matrices of all constituent elements in the structure, in the usual way as in the
case of the FEM, but the formulation leads to a non-linear eigenvalue problem and the
natural frequencies are generally extracted by applying the well-established algorithm of
Wittrick and William. Section B deals with fundamental research in dynamic stiffness
formulation and its application. Each chapter is self-contained and consists of a brief in-
troduction, literature review, theory, discussion of computed results and conclusions.
The layout of this Section B is as follows. Chapter 6 deals with the free vibration anal-
ysis of functionally graded beams and frameworks which is followed by Chapter 7 for
which the subject matter is free vibration of cracked beam. Next in Chapter 8, analyti-
cal development of advanced beam model incorporating Rayleigh-Love and Timoshenko
beam theories is given precedence. Finally Chapter 9 deals with the theoretical, analyti-
cal and computational development of the free vibration analysis of axial-bending coupled
beams.
Chapter 6
The free vibration analysis of functionally graded beams (FGBs) and frameworks con-
taining FGBs is carried out by applying the dynamic stiffness method and deriving the
elements of the dynamic stiffness matrix in explicit algebraic form. The stated rule that
the material properties of the FGBs vary continuously through the thickness according
to the power law forms the fundamental basis of the governing differential equations of
motion in free vibration. The differential equations are solved in closed analytical form
when the free vibratory motion is harmonic. The dynamic stiffness matrix is then for-
mulated by relating the amplitudes of forces to those of the displacements at the two
ends of the beam. Next, the explicit algebraic expressions for the dynamic stiffness ele-
ments are derived with the help of symbolic computation. Finally the Wittrick-Williams
algorithm is applied as solution technique to solve the free vibration problems of FGBs
with uniform cross-section, stepped FGBs and frameworks consisting of FGBs. Some
numerical results are validated against published results, but in the absence of published
results for frameworks containing FGBs, consistency checks on the reliability of results
are performed.
where Et and Eb are the Young’s moduli, and ρt and ρb are the densities at the top and
bottom surfaces of the beam, respectively.
In Equation 6.1, k (k ≥ 0) is the power law index parameter which indicates the material
property variation through the beam thickness. The parameter k has been extensively dis-
cussed in the literature [27, 28] and hence it is not elaborated here. However, three special
cases maybe observed. Clearly k = 1 indicates a linear variation of the composition be-
tween the top and bottom surfaces of the beam, k = 0 represents the case when the beam
is made of full material of the top surface whereas infinite k represents the case when the
beam is made of full material of the bottom surface.
where v and w are the corresponding displacements of a point on the neutral axis of the
beam. It should be noted that due to the variation of the material properties through the
thickness, the neutral axis would no longer be at the central line of the beam cross-section
[30, 31].
Using the displacement field given by Equation 6.2 and through the application of Hamil-
ton’s principle, the governing differential equations of motion in free vibration of the FGB
are given by [17, 19]
where Z Z
Ai = i
z E(z)dA, Bi = zi ρ(z)dA (i = 0, 1, 2) (6.4)
The natural boundary conditions from the Hamiltonian formulation [17, 19] give the fol-
lowing expressions for axial force F, shear force S and bending moment M as follows:
Clearly, due to the use of FGM, the axial (v) and bending motions (w) are coupled as
evident from Equations 6.3 and 6.5.
Assuming harmonic oscillation so that
where ω is the angular or circular frequency, V(y) and W(y) are the amplitudes of (v) and
(w) , respectively.
128 Chapter 6. Functionally graded beams
y
ξ= (6.7)
L
The differential equations of motion in Equations 6.3 can now be written as:
(B0 ω2 L3 + A0 LD2 )V(ξ) − (B1 ω2 L2 D + A1 D3 )W(ξ) = 0
(6.8)
(B1 ω2 L3 D + A1 LD3 )V(ξ) + (B0 ω2 L4 − B2 ω2 L2 D2 − A2 D4 )W(ξ) = 0
By combining the above two differential equations, it is possible to obtain a sixth order
ordinary differential equation satisfying both V(ξ) and W(ξ) to give:
where
H = V(ξ) or W(ξ) (6.10)
and
µ3 + aµ2 − bµ − c = 0 (6.13)
where
µ = λ2 (6.14)
Using an approach similar to the one described in [32, 33], it can be shown that the three
roots of the cubic equation (Equation 6.13) are real and the solution for H in Equation 6.9
can be expressed in terms of trigonometric and hyperbolic functions. This is advantageous
when deriving the explicit expressions for the dynamic stiffness elements of the FGB.
Explicit expressions are particularly useful when some, but not all of the stiffness elements
are needed, e.g. sensitivity analysis in optimisation studies. Thus if the roots [34] of
Equation 6.13 are α, β and γ, the solution for H is given by
where
" 1 #1 " 1 #1 " 1 #1
q 2 φ a 2 q 2 π−φ a 2 q 2 π+φ a 2
α= 2 cos − β= 2 cos + γ= 2 cos +
3 3 3 3 3 3 3 3 3
(6.16)
with
a2
q=b+ (6.17)
3
and
3
27c − 9ab − 2a
φ = cos−1 (6.18)
3
2 a2 + 3b 2
H (ξ) of Equation 6.15 represents the solution for both axial displacement V(ξ) and bend-
ing displacement W(ξ), containing different sets of constants as follows
The two different sets of constants Q1 -Q6 and R1 -R6 can be related with the help of any
one of the two of Equations 6.8 to give
k kβ kγ
Q1 = α
L
R2 , Q 3 = L
R4 , Q 5 = L
R6 ,
k kβ kγ
Q2 = L
α
R1 , Q 4 = − L
R3 , Q 6 = − L
R5 (6.21)
where
α(A α2 +B ω2 L2 ) β(A β2 −B ω2 L2 ) γ(A γ2 −B ω2 L2 )
kα = (A 1α2 +B 1ω2 L2 ) , kβ = (A 1β2 −B 1ω2 L2 ) , kγ = (A 1γ2 −B 1ω2 L2 ) (6.22)
0 0 0 0 0 0
With the help of Equations 6.5, 6.19 and 6.20, the expressions for bending rotation θ(ξ),
axial force F(ξ), bending moment M(ξ), and shear force S (ξ) can be obtained after some
simplification, as
′
W (ξ) 1
θ (ξ) = = {R1 α sinh αξ + R2 α cosh αξ − R3 β sin βξ + R4 β cos βξ
L L (6.23)
−R5 γ sin γξ + R6 γ cos γξ}
e
′ ′′
F (ξ) = − AL0 V − A1
A0 L
W = AL0 {−R1 eLα cosh αξ − R2 eLα sinh αξ + R3 Lβ cos βξ
e e e
+R4 Lβ sin βξ + R5 Lγ cos γξ + R6 Lγ sin γξ}
′′ A (6.24)
A1 L ′
M (ξ) = − AL22 W − A2
=−V {−R1 gα cosh αξ − R2 gα sinh αξ + R3 gβ cos βξ
2
L2
+R4 gβ sin βξ + R5 gγ cos γξ + R6 gγ sin γξ}
(6.25)
B2 L2 ω2 B1 L3 ω2
S (ξ) = L3 W + A2 W − A2 V − A2 V
A2 ′′′ ′ A1 L ′′
n o
= AL32 R1 fα sinh αξ + R2 fα cosh αξ + R3 fβ sin βξ − R4 cos βξ + R5 sin γξ − R6 cos γξ
(6.26)
130 Chapter 6. Functionally graded beams
where
At ξ = 0 : V = V1 , W = W1 , θ = θ1 , F = F1 , S = S 1 , M = M1 (6.30)
Figure 6.2: Sign convention for positive axial force F, shear force S and bending moment
M.
The displacement vector δ and the force vector P can be expressed as:
δ = [V1 W1 θ1 V2 W2 θ2 ]T , P = [ F1 S 1 M1 F2 S 2 M2 ] T (6.32)
6.3. Theoretical formulation 131
δ=BR (6.33)
where
kβ kγ
0 0 0
kα
L L L
1 0 1 0 1 0
0 β γ
α
0 0
(6.34)
L L L
B = kα S hα kα Chα k S kβ C β k S kγ C γ
− βL β − γL γ
L
L L L
Chα
S hα Cβ Sβ Cγ S γ
αS hα αChα βS βCβ γS γCγ
− Lβ − Lγ
L L L L
with
P = AR (6.36)
where
W1 eβ W1 eγ
0 0 0
− W1Leα L L
0 W3 f α 0 −W3 fβ 0 −W3 fγ
−W2 gα 0 W2 gβ 0 W 2 gγ 0
A = W1 eαChα W1 eα S hα W1 eβ Cβ W1 eβ S β W1 eγ Cγ W 1 eγ S γ
(6.37)
L L
− L − L − L − L
−W3 fα S hα −W3 fαChα −W3 fβ S β W3 fβCβ −W3 fγ S γ W3 fγCγ
W2 gαChα W2 gα S hα −W2 gβCβ −W2 gβ S β −W2 gγCγ −W2 gγ S γ
By eliminating the constant vector R from Equations 6.33 and 6.36, P and δ can be related
to give the dynamic stiffness matrix relationship of the FGB as
P=Kδ (6.38)
where
K = A B−1 (6.39)
is the required dynamic stiffness matrix with elements ki j (i = 1, 2, 3..6; j =1,2, 3,..6). With
the help of symbolic computation [20, 22], the matrix B of Equation 6.34 was inverted
algebraically and the inverted matrix was pre-multiplied by the matrix A of Equation
6.37 in order to generate the explicit expressions for each of the elements of the dynamic
132 Chapter 6. Functionally graded beams
stiffness matrix K. The stiffness expressions are simplified very considerably by means of
symbolic computations.
The above 6×6 frequency dependent dynamic stiffness matrix K of Equation 6.39 can now
be used to compute the natural frequencies and mode shapes of either an individual FGB
or an assembly of FGBs for different boundary conditions. A reliable and accurate method
of solving the problem is to apply the Wittrick-Williams algorithm [23] which is well
suited for the application of DSM. The algorithm uses the Sturm sequence property of the
dynamic stiffness matrix to ensure that no natural frequencies of the structure analysed
are missed. Basically the algorithm [23] gives the number of natural frequencies of a
structure that lie below an arbitrarily chosen trial frequency specified by the user. As
successive trial frequencies can be chosen by the user, this simple feature of the algorithm
can be exploited to bracket any natural frequency between its upper and lower bounds to
any desired accuracy. The results given in the next section were computed by applying
the Wittrick-Williams algorithm as the customary solution technique.
L2
r
ρb
λi = ωi (6.40)
h Eb
where ωi is the ith angular natural frequency in rad/s, ρb and Eb are density and Young’s
modulus of the bottom surface of the FGB.
6.4. Results and discussion
Table 6.1: Non-dimensional natural frequencies (λi ) of a uniform FGB with L/h = 10 and k = 0.5 for different boundary conditions.
q
ωi L 2 ρb
Boundary Non-dimensional natural frequencies λi = h Eb
condi- λ1 λ2 λ3 λ4 λ5
tions Present Ref.[30] Present Ref.[30] Present Ref.[30] Present Ref.[30] Present Ref.[30]
C-F 1.721 1.660 10.657 10.284 27.280 27.265 29.280 28.303 55.876 54.033
S-S 4.820 4.788 19.036 18.328 41.954 40.603 54.561 53.830 72.554 70.782
C-S 7.521 7.295 24.036 23.228 49.074 47.264 54.559 54.580 81.560 79.070
C-C 10.903 10.530 29.598 28.624 54.561 54.565 56.723 54.944 91.046 88.349
133
134 Chapter 6. Functionally graded beams
Figure 6.4: Natural frequencies and mode shapes of FGB with L/h = 10, k = 0.5 for
different boundary conditions (C-Clamped, F-Free, S-Simple support).
Table 6.1 shows the first five natural frequencies of the FGB with L/h = 10 and k = 0.5 for
C-F, S-S, C-S and C-C boundary conditions alongside the results reported in a recently
published paper [30]. The close agreement between the results from the current investi-
gation and the published ones is clearly evident. The maximum discrepancy between the
two sets of results is less than 4%. The corresponding mode shapes shown in Figure 6.4,
reveal that for the C-F and C-C boundary conditions, the first, second, fourth and fifth
6.4. Results and discussion 135
Table 6.2: Dimensionless fundamental natural frequency of a stepped FGB with k = 0.5
for different step ratios and boundary conditions.
modes are essentially bending (solid line) modes whereas the third one is axial (dashed
line). By contrast, for the S-S and C-S boundary conditions, the first, second, third and
fifth modes are basically bending modes and the fourth one axial.
The next set of results was obtained for a stepped beam (see Figure 6.5) made of FGM,
for which some comparative results are available in the literature. The DSM theory devel-
oped in this chapter can easily account for such problems with any step location, thickness
variation and boundary conditions. However, for brevity only the results for stepped FGB
with step locations L1 = 0.25L, L1 = 0.5L and L1 = 0.75L (see Figure 6.5) and the power
law index k = 0.5 are presented in Table 6.2 together with the published results [38]. Note
that for consistency, FGM type-II of Ref. [25] is used so that the results are directly com-
parable. Clearly, the results from the current investigation are in close agreement with
those of [25].
The final set of results was obtained for a portal frame consisting of three beam mem-
bers AB, BC and CD as shown in Figure 6.6. The natural frequencies of this frame are
available in the literature [35] when all the three beam members of the frame are made of
isotropic material and the supports at both points A and D are either clamped (built-in) or
pinned (simply-supported). These results are based on exact analytical theory. Using the
136 Chapter 6. Functionally graded beams
current theory, results are obtained for both the boundary conditions C-C and S-S at A and
D, respectively and making (i) all the three members AB, BC and CD isotropic (which
is achieved by substituting the power-law index parameter k to zero and using both the
top and bottom surface material properties to be the same and isotropic), (ii) AB and CD
isotropic, but BC made of FGM, (iii) BC isotropic and AB and CD made of FGM and
(iv) AB, BC and CD are all made of FGM. When computing numerical results, all three
members of the portal frame are assumed to have the same rectangular cross-section and
length. The width and depth (height) of the cross-section are taken to be 0.04 m and 0.02
m, respectively and length of each member is set to 1m. When any of the beam mem-
bers is isotropic, it is considered to be made of steel with Young’s modulus 200 GPa and
density 7500 kg/m3 whereas if it is made of FGM, the bottom surface is considered to
be steel with the above properties and the top surface ceramic with Young’s modulus 380
GPa and density 3960 kg/m3 . The computed natural frequencies are non-dimensionalised
with respect to the metallic properties to give
r
ρAL4
λi = ωi (6.41)
EI
The results of the investigation when all members of the portal frame (see Figure 6.6)
are metallic, are given in Table 6.3 showing the first three non-dimensional natural fre-
quencies of the frame when the points A and D are clamped (C-C) or simply-supported
(S-S), together with the results reported in Ref. [35]. The agreement between the sets
of results in Table 6.3 is excellent. Table 6.4 shows the three non-dimensional natural
frequencies when the vertical members AB and CD and the horizontal member BC of
the frame are made of either isotropic metal or FGM in turn, as indicated, and the points
A and D are clamped (C-C). Similar results are obtained for the case when the points A
and D of Figure 6.6 are simply-supported (S-S). The results for the S-S case are shown
in Table 6.5. Clearly the results shown in Tables 6.4 and 6.5 when compared to Table 6.3
indicate that significant changes in natural frequencies can occur as a result of using func-
tionally graded beams. This can have profound influence in the design of fire-resistant
multi-storey and multi-bay building structures.
6.4. Results and discussion 137
Table 6.3: Non-dimensional natural frequencies of a portal frame made of metallic mem-
bers.
q
4
Boundary Non-dimensional natural frequencies (λi = ωi ρAL EI
)
conditions λ1 λ2 λ3
at A and D Present Ref.[44] Present Ref.[44] Present Ref.[44]
C-C 3.204 3.205 12.639 12.648 20.627 20.629
S-S 1.463 1.463 9.866 9.870 14.854 14.856
Table 6.4: Non-dimensional natural frequencies of portal frame made of metal and FGM
with built-in boundary conditions at A and D for different k values.
q
4
Non-dimensional natural frequency (λi = ωi ρAL EI
)
Frame type with different categories of constituent members AB, BC and CD
k
AB, CD: metallic; BC: FGM BC: metallic; AB, CD: FGM AB, BC, CD: FGM
λ1 λ2 λ3 λ1 λ2 λ3 λ1 λ2 λ3
0.5 3.711 15.226 22.470 4.046 14.879 28.255 4.841 19.089 31.137
1.0 3.585 14.677 22.065 3.869 14.461 26.323 4.412 17.390 28.355
5.0 3.360 13.602 21.344 3.553 13.634 23.203 3.741 14.749 24.054
138 Chapter 6. Functionally graded beams
Table 6.5: Non-dimensional natural frequencies of portal frame made of metal and FGM
with simple-support boundary conditions at A and D for different k values.
q
4
Non-dimensional natural frequency (λi = ωi ρAL EI
)
Frame type with different categories of constituent members AB, BC and CD
k
AB, CD: metallic; BC: FGM BC: metallic; AB, CD: FGM AB, BC, CD: FGM
λ1 λ2 λ3 λ1 λ2 λ3 λ1 λ2 λ3
0.5 1.694 11.083 16.282 1.831 12.617 20.301 2.211 14.906 22.434
1.0 1.640 10.848 15.961 1.752 12.025 18.917 2.015 13.583 20.438
5.0 1.541 10.375 15.394 1.610 10.921 16.694 1.708 11.518 17.334
Figure 6.7: Natural frequencies and mode shapes of portal frame for various cases with
points A and D built-in. (a) AB, BC and CD are all metallic, (b) AB and CD are metallic,
but BC is FGM, (c) AB and CD are FGM, but BC is metallic, (d) AB, BC and CD are all
FGM.
6.4. Results and discussion 139
Figure 6.8: Natural frequencies and mode shapes of portal frame for various cases with
points A and D simply-supported. (a) AB, BC and CD are all metallic, (b) AB and CD
are metallic, but BC is FGM, (c) AB and CD are FGM, but BC is metallic, (d) AB, BC
and CD are all FGM.
For illustrative purposes, representative mode shapes for the portal frame of Figure 6.6
are presented in Figures 6.7 and 6.8 when the points A and D of the frame are built-in
(clamped) and simply-supported, respectively. The power law index parameter k is set
to 0.5 when computing the mode shapes. Figures 6.7(a) and 6.8(a) represent the mode
shapes when all three members of the portal frame are metallic whereas Figures 6.7(b)
and 6.8(b) shows the mode shapes when the columns AB and CD are metallic, but the
beam BC is made of FGM. By contrast Figures 6.7(c) and 6.8(c) show the mode shapes
for the case when the beam BC of the portal frame is metallic, but its columns AB and
CD are made of FGM. Finally, Figures 6.7(d) and 6.8(d) show the mode shapes when all
three members of the portal frame are made of FGM. Although the basic nature of the
mode shapes for the portal frame remains the same depending on the order of the natural
frequency on a case to case basis, significant changes in the natural frequencies are found
to occur when using FGM as evident from Figures 6.7 and 6.8. As expected the first mode
of the portal frame in each case is a sway mode with the frame oscillating between left and
right with virtually no elastic displacement of the central beam. The second mode shows
140 Chapter 6. Functionally graded beams
elastic deformations of all three members with no nodal point or any point of inflection
within any member. By contrast, the third mode reveals somehow a different picture in
that a node with zero displacement appears towards the top end of each columns whereas
a node for the beam appears near its centre. The mode shapes shown in Figures 6.7 and
6.8 are typical, as expected from the modal analysis of a portal frame and they are in
accord with similar mode shapes reported by other investigators [35, 36].
6.5 Conclusions
The dynamic stiffness matrix of a functionally graded beam is developed by deriving ex-
plicit expressions for the individual stiffness elements in explicit algebraic form. This is
achieved through the application of symbolic computation. The dynamic stiffness theory
is applied by using the Wittrick-Williams algorithm as solution technique to compute the
natural frequencies and mode shapes of some representative problems of uniform func-
tionally gradient beams, for which the material properties are considered to vary contin-
uously in the thickness direction according to a power law distribution. A stepped beam
made of functionally graded material is also investigated for its free vibration characteris-
tics. The results show good agreement with published results. Importantly, the theory has
been applied to study the free vibration behaviour of a portal frame with its constituent
members made of both isotropic and functionally graded material (FGM). The investi-
gation has shown that significant changes in the free vibration behaviour are possible by
using FGM. The literature on the free vibration of frameworks containing FGBs is virtu-
ally non-existent and this chapter addresses this problem and fills a gap in the literature by
analysing a portal frame. The developed theory can be applied to analyse high-rise build-
ing structures made of FGM which has advantageous mechanical properties of metal and
virtuous fire-resistant characteristics of ceramic and it is in this context, the investigation
carried out is expected to be most useful.
Chapter 7
Investigation into the static, dynamic and buckling behaviour of cracked beams has aroused
continuing interest amongst researchers. This is because the subject matter is of consid-
erable practical importance in engineering design to ensure safety and integrity of load
carrying structures that are vulnerable to cracks and other damages. The reduction in the
strength and stiffness properties of a structure due to the presence of single or multiple
cracks can be dangerous and may lead to catastrophic structural failures under both static
and dynamic loads. This has inspired the current research to carry out a parametric in-
vestigation into the free vibration characteristics of a cracked beam using the dynamic
stiffness method. As a fundamental prerequisite, first the dynamic stiffness matrix of a
cracked beam is formulated through the application of the compliance properties of the
crack. Basically, the cracked beam is idealised by connecting two Bernoulli-Euler beams
at the crack location where a local flexibility matrix representing the crack is introduced.
The dynamic stiffness matrix of the combined system is then developed. Once the dy-
namic stiffness matrix of the cracked beam is identified, the free vibration problem is
then formulated. The formulation leads to a non-linear eigenvalue problem for which
the Wittrick-Williams algorithm being ideally suited as solution technique, is applied to
yield natural frequencies and mode shapes of the cracked beam for different boundary
conditions. The first set of numerical results are obtained to validate the theory. This is
achieved by comparing results from the current theory with those available in the pub-
lished literature. Next a detailed parametric study is carried out by changing the crack
location, crack depth and boundary conditions of the beam. The results are discussed and
the chapter concludes with some significant remarks.
142 Chapter 7. Free vibration of cracked beam
for which some comparative results are available in the literature is analysed to confirm
the validity and accuracy of the proposed theory.
The cracked beam shown in Figure 7.1 is idealised into three structural elements denoted
by letters (a), (b) and (c) having lengths L1 , L2 and L3 , respectively, as shown in Figure
144 Chapter 7. Free vibration of cracked beam
Figure 7.2: Node numbering and member (element) lettering of a cracked beam.
7.2. Thus, the cracked beam is represented by two Bernoulli-Euler beam elements (a)
and (b), connected together by a crack element (c). In addition to the element (member)
lettering, the figure also shows the node numbering of the entire cracked beam. As shown
in the Figure 7.2, nodes 1 and 2 are connected by a beam element (a), nodes 2 and 3 are
connected by a crack element (c) whereas nodes 3 and 4 are connected by a beam element
(b). For clarity and subsequent development of the theory, the (fictitious) length of the
crack element is assumed to be L3 instead of zero. This fictitious length is irrelevant and
have no consequence on the theory. This is because the crack is located at a point on the
axis of the beam and the flexibility matrix introduced at that point is not considered to
be an explicit function of length. The associated stiffness matrix of the element (c) will
be computed from the compliance matrix at the crack, which will eventually turn out to
be a nodal stiffness matrix rather than a conventional element stiffness matrix connecting
two nodes and separated by a distance. The crack element’s compliance (flexibility) and
subsequent stiffness properties are established using fracture mechanics theory. At the
crack location, the construction of a nodal stiffness matrix with the help of the compliance
properties resulting from the crack is an essential part of the present theory. An in house
program is developed by using dynamic stiffness matrix for the crack element and then by
using the Wittrick-Williams Algorithm, it became possible to ascertain how many natural
frequencies of a structure lie below an arbitrarily chosen trial frequency. This simple
feature of the algorithm is exploited to advantage to converge upon any required natural
frequency to any desired accuracy.
The dynamic stiffness matrix of a freely vibrating beam element (e) connecting nodes i
and j, see Figure 7.3, relates the amplitudes of forces pi and pj to those of the displace-
ments δi and δj at the ends i and j of the element. In the usual matrix notation and referring
7.3. Theoretical formulation 145
Figure 7.3: Forces and displacements at the ends of a beam element e connecting nodes
i and j.
where
n oT n oT
pi = p xi , pyi , mi p j = p x j , py j , m j (7.2)
n oT n oT
δi = δ xi , δyi , θi δ j = δ x j , δy j , θ j
p xi and p x j are axial force, pyi and py j are shear forces and mi and m j are bending moments,
T denotes a transpose.
The dynamic stiffness elements of the sub-matrices ke 11 and ke 12 of Equation 7.1 are
frequency dependent and are given by [60]
a1 = (EA/L)ν cot ν; b1 = (EIλ3 /L3 )(sin λ cosh λ + cos λ sinh λ)/Z (7.4)
c1 = (EIλ/L)(sin λ cosh λ − cos λ sinh λ)/Z; d1 = (EIλ2 /L2 ) sin λ sinh λ/Z (7.5)
with
ν = ωL(m/EA)1/2 ; λ = (mω2 L4 /EI)1/4 ; Z = 1 − cos λ cosh λ (7.8)
where ω is the circular (or angular) frequency, and L, m, EA and EI are the element
(member) length, mass per unit length, extensional (axial) rigidity and flexural rigidity,
146 Chapter 7. Free vibration of cracked beam
respectively.
ke 22 in Equation 7.1 can be obtained from ke 11 by writing –d1 for d1 in Equation 7.3 and
ke 21 is the transpose of ke 12 .
For a given trial frequency, it is thus possible to compute the dynamic stiffness elements
of members (a) and (b) in Figure 7.2 by substituting the appropriate lengths L1 and L2 in
Equations 7.1-7.8, respectively.
The stiffness matrix for the crack element (c) in Figure 7.2 can be obtained from the 3×3
flexibility matrix C given as
C11 0 0
C = 0 C22 0
(7.9)
0 0 C33
There are literally dozens of papers which deal with the derivation of the flexibility matrix
by integrating the stress intensity factor. In particular, Zheng and Kessissoglou [53] have
given explicit expressions for the elements C11 , C22 and C33 for both rectangular and
circular cross-section beams in terms of the cross-sectional dimensions (width b and depth
h for a rectangle, and radius r for a circle) and the crack length a through the depth.
For plane strain problems and for a rectangular cross-section shown in Figure 7.4, these
expressions for crack-length over depth ratios within the range 0≤ a/h ≥ 0.5 are taken
from [53] and written as follows
1 − µ2 1 − µ2 1 − µ2
c11 = F (1, 1) c22 = F (2, 2) c33 = F (3, 3) (7.10)
Eb Eb Eb
where µ is the Poisson’s ratio, E is the Young’s modulus and F(1,1), F(2,2) and F(3,3) are
non-dimensional functions given by [53]
−34.954632ξ9 + 9.054062ξ10
(7.11)
−0.326018 × 10−6 ξ + 1.454954ξ2 − 1.455784ξ3 − 0.421981ξ4
1
−0.279522ξ5 + 0.455399ξ6 − 2.432830ξ7 + 5.427219ξ8
F (2, 2) = e (1−ξ)
−6.643057ξ9 + 4.466758ξ10
(7.12)
−0.219628 × 10−4 ξ + 52.37903ξ2 − 130.2483ξ3 + 308.442769ξ4
1
F (3, 3) = e (1−ξ) −602.445544ξ5 + 939.044538ξ6 − 1310.95029ξ7 + 1406.52368ξ8
−1067.4998ξ9 + 391.536356ξ10
(7.13)
where
ξ = a/h (7.14)
7.3. Theoretical formulation 147
where
Note that the stiffness matrix kc corresponds to a crack element of zero length. Thus the
length L3 appearing in Figure 7.2 for the crack element is irrelevant and hence can be
disregarded.
Referring to Figure 7.2 and using Equations 7.1-7.13, the frequency dependent dynamic
stiffness matrix K(ω) of the cracked beam can now be formed in the usual way by assem-
bling the element stiffness matrices of two Bernoulli-Euler beams (a) and (b) (connecting
nodes 1 and 2, and 3 and 4, respectively) and one crack element (c), connecting nodes 2
and 3). In matrix notation, the assembled stiffness matrix K(ω) that will lead to an eigen-
value problem can be expressed symbolically as follows.
148 Chapter 7. Free vibration of cracked beam
0 0
a a
k11 k12
k a a
k +k c
k c
0
K(ω) = 21 22 c 11 c 12 b
(7.19)
0 k21 b
k22 + k11 k12
0 0
b b
k21 k22
Appropriate boundary conditions at nodes 1 and 4 in Figure 7.2, can be applied by deleting
the particular rows and columns of K(ω), corresponding to zero displacements in order to
compute the natural frequencies and mode shapes of individual cases such as cantilever,
simply-supported and clamped-clamped cracked beams. (Note that for a free-free cracked
beam, the whole K(ω) matrix must be used for the eigenvalue problem.) The dynamic
stiffness matrix of Equation 7.19 can now be used to compute the natural frequencies and
mode shapes of cracked beams with various end conditions. A non-uniform cracked beam
can also be analysed for its free vibration characteristics by idealising it as an assemblage
of many uniform cracked beams. An accurate and reliable method of calculating the
natural frequencies and mode shapes of a structure using the dynamic stiffness method is
to apply the well-known algorithm of Wittrick and Williams [23] which has been featured
in numerous papers.
gives higher value than Simple-simple and Cantilever cases and Simple-Simple support
gives higher value than Cantilever for all non-dimensional crack lengths as can be seen
from Table 7.1. Similar results are observed for all other Tables 7.2, 7.3 and 7.4 and it can
be seen as crack location ζ increases the value of natural frequencies increases gradually
too with highest value occurs at crack location ζ = 0.8.
150
Table 7.1: Natural frequencies with various boundary conditions at crack location, ζ = 0.2.
Table 7.4: Natural frequencies with various boundary conditions at crack location, ζ = 0.8.
151
152 Chapter 7. Free vibration of cracked beam
Table 7.5: Natural frequencies with various boundary conditions for intact beam.
Table 7.5 shows the value of first three natural frequencies for various boundary conditions
of an intact (no cracks) beam. The result of the cantilever case is compared with Ref [51]
where it can be observed that the results obtained using the current theory are very close
to the published results (The maximum discrepancy is around 1.8%). It should be noted
that the results obtained by using the dynamic stiffness method is exact and the small
difference that exists with the reference [51] may be attributed to the fact that the present
theory is based on an exact dynamic stiffness method whereas the theory used in Ref [51]
is an approximate one (i.e. FEM).
Figure 7.5: The natural frequency ratio between the cracked beam and the intact beam for
the fundamental mode of a cantilever cracked beam having non-dimensional crack length
ξ = 0.4.
Figure 7.5 shows the natural frequency ratio between the crack beam and the intact beam
7.4. Results and discussion 153
for the fundamental natural frequency ω1 relative to that of the intact beam ω01 of a can-
tilever cracked beam when the non-dimensional crack length is ξ = 0.4. It can be seen that
the fundamental frequency drops as non-dimensional crack length increases, which is ob-
served in Table 7.2 too. So it can be said that the effect of non-dimensional crack length is
significant on the fundamental natural frequency of the cantilever cracked beam. Similar
observations were made for higher natural frequencies and for other boundary conditions
which was not shown in the form of figures for brevity but can be observed from Tables
7.1-7.4. These observations were also reported by a number of earlier investigators [51,
52, 61].
The final set of results was obtained to illustrate the mode shapes of the cracked beam.
The first five natural frequencies and mode shapes of the cantilever are shown in Figure
154 Chapter 7. Free vibration of cracked beam
7.6 when the crack is located at ζ = 0.6 and the non-dimensional crack length is ξ = 0.4.
The bending displacements for the cracked beam and the intact beam are shown such that
bending is represented by black solid lines and torsion mode is shown as red dashed lines.
7.5 Conclusions
Using the dynamic stiffness method and Bernoulli-Euler beam theory, the free vibration
behaviour of a cracked beam is investigated. In the development of the theory, the cracked
beam is idealised by two intact beams joined together and by introducing a flexibility ma-
trix of the crack at the joint. The formulation leads to a non-linear eigenvalue problem
and was solved by applying the Wittrick-Williams algorithm. A parametric study by
varying the crack location, crack depth and boundary conditions is carried out and the
results showing the natural frequencies and mode shapes are illustrated. The investigation
has shown that the crack location, crack depth and boundary conditions have significant
effects on the natural frequencies and modes shapes of a cracked beam. The natural fre-
quencies of intact beam for cantilever boundary condition are readily available in the
existing literature and in this research natural frequencies of intact beam for simply sup-
ported, clamped-pinned and clamped-clamped boundary conditions have been obtained
for cracked beams and contrasted against the results of intact beams. The theory devel-
oped can be extended to frame works and other building structures and within the context
of structural health monitoring purposes, the developed theory is expected to be most
useful.
Chapter 8
In this chapter, an exact dynamic stiffness matrix for a beam is developed by integrating
the Rayleigh-Love theory for longitudinal vibration into the Timoshenko theory for bend-
ing vibration. It should be noted that no one appears to have made any attempt to combine
Rayleigh-Love bar and Timoshenko beam theories particularly when investigating the free
vibration characteristics of frameworks. This will be important within the high frequency
range when using the SEA technique. In the formulation, the Rayleigh-Love theory ac-
counted for the transverse inertia in longitudinal vibration whereas the Timoshenko beam
theory accounted for the effects of shear deformation and rotating inertia in bending vi-
bration.
The dynamic stiffness matrix is developed by solving the governing differential equations
of motion in free vibration of a Rayleigh-Love bar and a Timoshenko beam and then
imposing the boundary conditions for displacements and forces. Next the two dynamic
stiffness theories are combined using a unified notation. The ensuing dynamic stiffness
matrix is subsequently used for free vibration analysis of uniform and stepped bars as
well as frameworks through the application of the Wittrick-Williams algorithm as solu-
tion technique. Illustrative examples are given to demonstrate the usefulness of the theory
and some of the computed results are compared with published ones and this chapter
closes with some concluding remarks.
156 Chapter 8. Combined Rayleigh-Love and Timoshenko theories
of the Rayleigh-Love bar is integrated with the dynamic stiffness matrix of a Timoshenko
beam [76–78] which accounts for the effects of shear deformation and rotatory inertia to
allow for the free vibration analysis of individual members and plane frames in the low,
medium and high frequency range through the application of the Wittrick-Williams algo-
rithm [23] as solution technique. Using the developed theory, a wide range of problems
is solved and some of the computed results are compared with published literature.
It should be noted that there are no specific hard boundaries between the regimes of low,
medium and high frequencies, but a useful descriptor which gives an indicative guidance
to frequency range is the vibrational wavelength when compared to the overall length
of the structure. Thus an engineering judgement can be reasonably made based on the
product of the wave number and a typical length of the structure, which is essentially
the Helmholtz number. Large values of this number represent the high frequency range
whereas lower values determine the low to medium frequency range. For the type of prob-
lems investigated in this chapter, the low to medium range of frequencies is characterised
to be below 1500 Hz whereas frequencies above this value constitute the high frequency
range.
bar in free axial (or longitudinal) vibration can be derived by using Hamilton’s principle
as the first step towards the dynamic stiffness formulation. The focus area of the derivation
in this section is, of course, on the axial stiffnesses only.
Figure 8.1: Coordinate system and notation for a Rayleigh-Love bar and a Timoshenko
beam.
Referring to Figure 8.1 and noting that if u is the axial displacement at a distance x from
the origin, the kinetic and potential energies of the bar T bar and Vbar are respectively given
by [71, 79]
Z !2 2
!2
1 L ∂u 2 ∂ u
T bar = ρA + ρIP ν dx (8.1)
2 0 ∂t ∂x∂t
and !2
1 L
Z
∂u
Vbar = EA dx (8.2)
2 0 ∂x
where ρ is the density of the bar material, A is the cross-sectional area of the bar so that
ρA represents the mass per unit length, IP is the polar second moment of area so that ρIP
represents the polar mass moment of inertia per unit length, E is the Young’s modulus of
the bar material so that EA represents the axial or extensional rigidity of the bar and ν is
the Poisson’s ratio of the bar material.
Hamilton’s principle states Z t2
δ (T bar − Vbar )dt = 0 (8.3)
t1
where t1 and t2 are the time interval in the dynamic trajectory, and δ is the usual variational
operator.
The governing differential equations of motion of the Rayleigh-Love bar and the associ-
ated boundary condition in free vibration can now be derived by substituting the kinetic
(T bar ) and potential (Vbar ) energy expressions of Equations 8.1 and 8.2 into Equation 8.3,
8.3. Development of dynamic stiffness formulation 159
using the δ operator, integrating by parts and then collecting terms. In an earlier publica-
tion, the entire procedure to generate the governing differential equations of motion and
natural boundary conditions for bar or beam type structures was automated by Banerjee et
al [22] by applying symbolic computation. In this way, the governing differential equation
of motion of the Rayleigh-Love bar is obtained as [71, 79]
∂2 u ∂2 u 4
2 ∂ u
EA − ρA + ρI P ν =0 (8.4)
∂x2 ∂t2 ∂x2 ∂t2
As a by-product of the Hamiltonian formulation, the expression for the axial force f (x, t)
follows from the natural boundary condition to give [71, 79]
∂u ∂3 u
f (x, t) = −EA − ρIP ν2 (8.5)
∂x ∂x∂t2
If harmonic oscillation is assumed, then
where ω is the angular or circular frequency, and U(x) are the amplitudes of u.
Substituting Equation 8.6 into Equation 8.5 gives
d2 U
(EA − ρIP ν2 ω2 ) 2
+ ρAω2 U = 0 (8.7)
dx
As a result of the harmonic oscillation assumption, the amplitude F(x) of the force f (x, t)
in Equation 8.5 becomes
dU
F (x) = −(EA − ρIP ν2 ω2 ) (8.8)
dx
Introducing the differential operator D = d/dξ and the non-dimensional length ξ as
x
ξ= (8.9)
L
Equation 8.7 becomes
(D2 + γ2 )U = 0 (8.10)
where
α2
γ2 = (8.11)
1 − β2
with
ρAω2 L2 ρIP ν2 ω2
α2 = ; β2 = (8.12)
EA EA
The expression for the amplitude of the axial force in Equation 8.8 using Equations 8.9
and 8.12 becomes
EA dU
F (ξ) = − (1 − β2 ) (8.13)
L dξ
160 Chapter 8. Combined Rayleigh-Love and Timoshenko theories
Figure 8.2: Boundary conditions for displacements and forces in axial vibration for a
Rayleigh-Love bar.
Substituting Equations 8.16 and 8.17 into Equations 8.14 and 8.15, the following matrix
relationships can be obtained
∆ x1 0 1 C1
= (8.18)
∆ x2 sinγ cosγ C2
and
F x1 EA 2
−1 0 C1
= γ(1 − β ) (8.19)
F x2 L cosγ −sinγ C2
The constants C1 and C2 can now be eliminated from Equations 8.18 and 8.19 to give the
dynamic stiffness matrix of an axially vibrating Rayleigh-Love bar relating amplitudes of
the forces and displacements at its ends as follows:
F x1 a1 a2 ∆ x1
= (8.20)
F x2 a2 a1 ∆ x2
where the elements of the 2×2 dynamic stiffness matrix are given by
EA EA
γ 1 − β2 cotγ , γ 1 − β2 cosecγ
a1 = a2 = − (8.21)
L L
8.3. Development of dynamic stiffness formulation 161
It should be noted that the Rayleigh-Love theory has a limitation that β2 in Equations
8.11 and 8.12 must be less than one which is usually the case, otherwise, the solution of
Equation 8.10 would not be harmonic and hence no oscillatory motion will take place.
This limitation has been pointed out in the literature, e.g. see Equation 13 of [80].
∂w
=θ+γ (8.24)
∂x
or
∂w
−θ γ= (8.25)
∂x
Thus, the potential energy Vbeam of Equation 8.22 becomes
!2 !2
1 L
1 L
Z Z
∂θ ∂w
Vbeam = EI dx + kAG − θ dx (8.26)
2 0 ∂x 2 0 ∂x
Substituting the expressions for the kinetic and potential energies T beam and Vbeam from
Equations 8.22 and 8.26 into Hamilton’s principle expressed in Equation 8.3 and then
integrating by parts and collecting terms yield the governing differential equations of
162 Chapter 8. Combined Rayleigh-Love and Timoshenko theories
motion and the associated boundary conditions providing the expressions for bending
moment (M) and shear force (S) as follows [81]. Governing differential equations
∂2 w
!
∂ ∂w
− ρA 2 + kAG −θ =0 (8.27)
∂t ∂x ∂x
∂2 θ ∂2 θ
!
∂w
− ρI 2 + EI 2 + kAG −θ (8.28)
∂t ∂x ∂x
Natural boundary conditions Shear force:
∂2 θ ∂2 θ
!
∂w
v = −kAG − θ = EI 2 − ρI 2 (8.29)
∂x ∂x ∂t
Bending moment:
∂θ
m = −EI (8.30)
∂x
Introducing the non-dimensional length ξ = x/L and assuming harmonic oscillation so
that
w (x, t) = W(ξ)eiωt (8.31)
θ (x, t) = Θ(ξ)eiωt (8.32)
where W(ξ) and Q(ξ) are the amplitudes of the bending displacement and bending rotation
of the harmonically vibrating Timoshenko beam.
Equations 8.27 and 8.28 can now be combined to give a fourth order ordinary differential
equation as follows which is identically satisfied by both W(ξ) and Q(ξ)
D4 + b2 r2 + s2 D2 − b2 1 − b2 r2 s2 H = 0
h i
(8.33)
where
d 1 d
D= = (8.34)
dξ L dx
4
ρAω2 L I EI
b2 = ; r2 = ; s2 = (8.35)
EI AL2 kAGL2
and
H = W or Θ (8.36)
If a trial solution H = eλξ is assumed where λ is a constant, yet to be known, the auxiliary
or characteristic equation of the differential Equation 8.33 is given by
λ4 + b2 r2 + s2 λ2 − b2 1 − b2 r2 s2 = 0
(8.37)
q
b2 r2 − s2 2 + 4b2
−b2 r2 + s2 ±
= (8.38)
2
Clearly λ2 will be always real and for the negative value of the expression under the square
root sign of Equation 8.38, one of the two values of λ2 will be always negative, resulting in
two imaginary roots of λ which will lead to part of the solution of Equation 8.33 in terms
of trigonometric functions whereas the other value of λ2 when using the positive value
before the square root sign can be either positive or negative depending on whether the
square root expression in Equation 8.38 is bigger than or smaller than b2 (r2 + s2 ). If this
second value of λ2 is positive which is usually the case, the two roots of λ2 will be real,
yielding the remaining solution of Equation 8.33 in terms of hyperbolic functions so that
the two of the four integration constants in the solution will be connected to trigonomet-
ric functions and the other two to hyperbolic functions. However, for exceptionally high
frequencies or for exceptionally squat beams, the second value of λ2 can be negative like
the first one which will give the entire solution of Equation 8.33 in terms of trigonometric
functions only. The two sets of solutions and their conditionality are explained below.
The expression for λ2 in Equation 8.38 can be expressed in the following alternative form
b2
r
2 + s2 2 +
4
λ2 = 2 2
1 2 r 2 s2 (8.39)
− r + s ± r − b
2 b2
It is clear from Equation 8.39 that if b2 r2 s2 < 1, one of the values of λ2 will be negative
and the other value will be positive whereas if b2 r2 s2 > 1, they both will be negative.
Thus the solutions for bending displacement W and bending rotation Θ for these two
cases resulting from the differential equation of Equation 8.33 are given by
1. b2 r2 s2 < 1
2. b2 r2 s2 > 1
W (ξ) = A1 cosΦ + A2 sinΦ + A3 cosΛ+ A4 sinΛ (8.42)
Θ (ξ) = B1 cosΦ + B2 sinΦ + B3 cosΛ + B4 sinΛ (8.43)
where
b2 r2 + s2 b2
r
4
Φ2 = r2 + s2 2 + 2 1 − b2 r2 s2 (8.44)
+
2 2 b
and
b2 r2 + s2 b2
r
4
jΛ2 = − r2 + s2 2 + 2 1 − b2 r2 s2 (8.45)
+
2 2 b
164 Chapter 8. Combined Rayleigh-Love and Timoshenko theories
with
j = 1 f or b2 r2 s2 < 1; j = −1 f or b2 r2 s2 > 1 (8.46)
and A1 - A4 and B1 - B4 are two different sets of constants.
It should be noted from Equation 8.35 that
ρIω2
b2 r 2 s 2 = (8.47)
kAG
For most of the practical problems, b2 r2 s2 will be less than one unless ω is exceptionally
large. This is because the shear rigidity kAG is generally much bigger than the rotatory
inertia per unit length ρI for any realistic cross-section and beam material, but neverthe-
less, the solutions given by Equations 8.42 and 8.43 are included in the theory to cover
the exceptional case when b2 r2 s2 is greater than one.
With the help of Equation 8.27 or 8.28 and the solution given by Equations 8.40-8.43, it
can be shown that the two sets of constants A1 - A4 and B1 - B4 are related. Using Equation
8.27, the following relationships between B1 - B4 and A1 - A4 are obtained.
kΦ kΦ kΛ kΛ
B1 = A2 ; B2 = − A1 ; B3 = A4 ; B4 = j A3 (8.48)
L L L L
where
Φ2 − b2 s2 Λ2 + jb2 s2
! !
kΦ = ; kΛ = (8.49)
Φ Λ
Because of the harmonic oscillation hypothesis adopted for the freely vibrating Timo-
shenko beam as indicated by Equations 8.31 and 8.32 and also by the introduction of the
non-dimensional length ξ = x/L, the expressions for the amplitudes of the shear force
(V) and bending moment (M) arising from Equations 8.29, 8.30 and 8.48 will take the
following form.
EI d2 Θ
!
2 2 EI
V= 2 − b r Θ = (A1 eΦ sinΦξ − A2 eΦ cosΦ + jA3 eΛ sinΛξ + A4 eΛ cosΛξ )
L dξ2 L3
(8.50)
EI dΘ EI
M=− 2 = − 2 (−A1 ΦkΦ cosΦξ − A2 ΦkΦ sinΦξ + jA3 ΛkΛ cosΛξ + jA4 ΛkΛ sinΛξ )
L dξ L
(8.51)
where
eΦ = Φ2 − b2 r2 kΦ ; eΛ = j Λ2 + jb2 r2 kΛ
(8.52)
and j and kφ , kΛ have already been defined in Equations 8.46 and 8.49, respectively.
Now from the expressions for the amplitudes of displacements W and Θ given by Equa-
tions 8.40 - 8.43 and the corresponding forces V and M given by Equations 8.50 and 8.51,
the dynamic stiffness matrix of the Timoshenko beam can be formulated by applying
the boundary conditions in algebraic form relating the amplitudes of forces and displace-
ments.
8.3. Development of dynamic stiffness formulation 165
Referring to Figure 8.3, the boundary conditions for the displacements and forces can be
applied as follows
Figure 8.3: Boundary conditions for displacements and forces for a Timoshenko beam.
Substituting Equations 8.53 and 8.54 into Equations 8.40 - 8.43 and Equations 8.50 and
8.51, the following two matrix equations are obtained for displacements and forces, re-
spectively, in terms of the constants A1 - A4 .
0 0
Fy1 −W 3 eΦ W 3 eΛ A1
M1 W2 ΦkΦ 0 − jW 2 ΛkΛ 0 A2
Fy2
=
−W3 eΦ S
A3 (8.57)
W 3 eΦ C − jW3 eΛ S −W 3 eΛC
M2 −W 2 ΦkΦC −W 2 ΦkΦ S jW2 ΛkΛC jW 2 ΛkΛ S A4
or
F = RA (8.58)
where
S = sinΦ; C = cosΦ (8.59)
S = sinh Λ; C = cosh Λ b2 r2 s2 < 1( j = 1)
166 Chapter 8. Combined Rayleigh-Love and Timoshenko theories
and
W1 , W2 and W3 are defined as follows
EI EI EI
W1 = ; W2 = 2 ; W3 = (8.61)
L L L3
The constants A1 - A4 can now be eliminated from Equations 8.55 and 8.57 to give the
4×4 dynamic stiffness matrix of the Timoshenko beam. This can be achieved by inverting
the square matrix of Equation 8.55, i.e. Q matrix of Equation 8.56 and pre-multiplying it
with the square matrix of Equation 8.57, i.e. R matrix and performing the matrix operation
RQ−1 numerically to give the dynamic stiffness matrix. Alternatively, the matrix inversion
and matrix multiplication procedures can be carried out symbolically (algebraically) to
generate explicit expressions for each of the stiffness elements of the dynamic stiffness
matrix to give.
Fy1 d1 d2 d4 d5 ∆y1
M1 = d2 d3 −d5 d6 Θ1 (8.62)
F d
y2 4 d5 d1 −d2 ∆y2
M2 −d5 d6 −d2 d3 Θ2
where
d1 = W3 b2 Γ(CS + ηS C)/(ΛΦ)
n o
d2 = W2 ZΓ (Φ + jηΛ) S S − (Λ − ηΦ) 1 − CC /(Λ + ηΦ)
d3 = W1 Γ S C − jηCS
d5 = W 2 ZΓ(C − C)
d6 = W1 Γ jηS − S
with
Z = Φ − b2 s2 /Φ; η = Z/( jΛ + b2 s2 /Λ);
Figure 8.4: Amplitudes of displacements and forces at the ends of a combined Rayleigh-
Love bar and a Timoshenko beam.
Timoshenko beam incorporating the bending stiffness to enable the free vibration analysis
of plane frames to be made.
Referring to Figure 8.4 and Equations 8.20 and 8.62, the resulting dynamic stiffness ma-
trix is given by
j = j0 + s{K f } (8.67)
where Kf ,the overall dynamic stiffness matrix of the wing whose elements depend on
ω and is evaluated at ω = ω∗ ; s{Kf } is the number of negative elements on the leading
diagonal of K∆f , K∆f is the upper triangular matrix obtained by applying the usual form of
Gauss elimination to Kf , and j0 is the number of natural frequencies of the wing still lying
between ω = 0 and ω = ω∗ when the displacement components to which Kf corresponds
are all zeros. (Note that the structure can still have natural frequencies when all its nodes
are clamped, because exact member equations allow each individual member to displace
between nodes with an infinite number of degrees of freedom, and hence infinite number
of natural frequencies between nodes.) Thus
j0 = Σ jm (8.68)
any frequency of practical interest. Nevertheless, j0 count of the algorithm is not really a
peripheral issue, particularly for achieving computational efficiency and avoiding further
unnecessary discretisation of the structure. The need to compute j0 stems from the fact
that the DSM allows infinite number of natural frequencies to be accounted for when all
the nodes of the structure are fully restrained and yet one or more structural members can
vibrate freely on their own between the nodes resulting in δ = 0 modes in the eigenvalue
equation [KD ]{δ} = 0.
Thus, proceeding in the same way as in the case of classical Bernoulli-Euler bar [60] the
number of clamped-clamped natural frequencies jR of a Rayleigh-Love bar lying below
an arbitrarily chosen trial frequency ω∗ is given by
γ
jR = highest integer < (8.70)
π
jc = jd f or b2 r2 s2 < 1
(8.72)
jc = jd + je f or b2 r2 s2 ≥ 1
with
φ
jd = highest integer <
π (8.73)
Λ
je = highest integer < +1
π
In Equation 8.73, φ and Λ have already been defined in Equations 8.44 and 8.45 re-
spectively. Thus the number of clamped-clamped natural frequencies jm exceeded by an
170 Chapter 8. Combined Rayleigh-Love and Timoshenko theories
individual member by the trial frequency ω∗ with the inclusion of the Rayleigh-Love bar
and the Timoshenko beam theories is given by
jm = jR + jT (8.74)
Now the root count j0 of Equation 8.68 can be computed using the Equation 8.68 where
the summation Σ over m is extended to include all elements in the structure.
The ratio between the natural frequencies for the clamped-clamped bar obtained from the
Rayleigh-Love and classical Bernoulli-Euler theories can be expressed with the help of
Equations 8.75 and 8.76 to give
ωn 1
= r (8.77)
ωn0 ν2 n2 π2
1+
( Lr )2
8.5. Results and discussion 171
r
Ip
r= (8.78)
A
Proceeding in a similar way and imposing appropriate boundary conditions, the natural
frequency ratio for a cantilever bar using the Rayleigh-Love and classical Bernoulli-Euler
theories can be expressed as
ωn 1
= r (8.79)
ωn0 (2n−1)2 π2 ν2
1+ 2
4( Lr )
The validity of the Equations 8.77 and 8.79 has been further confirmed by using the de-
veloped dynamic stiffness matrix of a Rayleigh-Love bar shown in Equation 8.20.
Clearly Equations 8.77 and 8.79 indicate that the natural frequency ratio ωωnn is dependent
0
on the Poisson’s ratio ν of the bar material as well as the slenderness ratio L/r of the bar.
The Poisson’s ratio ν for an isotropic material is generally constant and maybe assumed
to be 0.3 which is used here in the analysis.
Figures 8.5 and 8.6 shows the variation of the ratio of the first five natural frequencies us-
ing the Rayleigh-Love and classical Bernoulli-Euler theories against the slenderness ratio
L/r for the clamped-clamped and cantilever bar respectively. Clearly for smaller values
of slenderness ratios and for higher natural frequencies, the classical Bernoulli-Euler the-
ory is considerably less accurate. The errors incurred in the fifth natural frequency when
using the classical Bernoulli-Euler theory are 27% and 24% for the clamped-clamped and
cantilever bar respectively when the slenderness ratio is 5. It should be noted that in the
Statistical Energy Analysis (SEA) for which modal density in the high frequency range is
required, the classical Bernoulli-Euler theory can be inadequate.
172 Chapter 8. Combined Rayleigh-Love and Timoshenko theories
Figure 8.5: The first five natural frequency ratios using the Rayleigh-Love and classi-
cal Bernoulli-Euler theories for a clamped-clamped bar in axial vibration. ωn = natural
frequency using Rayleigh-Love theory; ωn0 = natural frequency using classical Bernoulli-
Euler theory.
Figure 8.6: The first five natural frequency ratios using the Rayleigh-Love and clas-
sical Bernoulli-Euler theories for a cantilever bar in axial vibration mode. ωn = natural
frequency using Rayleigh-Love theory; ωn0 = natural frequency using classical Bernoulli-
Euler theory.
8.5. Results and discussion 173
The numerical values for the data taken from [92] are: r1 = 0.05m, r2 = 0.03m, r3
= 0.075m, l1 = 0.05m, l2 = 0.17m, l3 = 0.13m, E1 = 200×109 Pa, E2 = 70×109 Pa,
E3 =100×109 Pa, ρ1 = 7.85×103 kg/m3 , ρ2 = 2.7×103 kg/m3 , ρ3 = 8.4×103 kg/m3 , ν1 =
0.30, ν2 = 0.33, ν3 = 0.34.
The first four natural frequencies computed using the Rayleigh-Love dynamic stiffness
theory are shown in column 2 of Table 8.1 alongside the results reported in [82] shown
in column 3. The corresponding natural frequencies computed using classical Bernoulli-
Euler dynamic stiffness theory [60] are also shown in the parenthesis in column 2. Al-
though the agreement of the results between the present theory and those of [82] are good
for the second and fourth natural frequencies (the differences are well within 3%), but
for the first and third natural frequencies there are some discrepancies which are around
13% and 15% respectively. The fundamental natural frequency of the bar quoted in [82]
is well above the corresponding natural frequency obtained from the classical Bernoulli-
Euler theory. This is surely in error because the effect of the transverse inertia presumably
accounted for in [82] is expected to diminish the natural frequency and not increase it. The
mode shapes corresponding to the four natural frequencies using the present theory are
174 Chapter 8. Combined Rayleigh-Love and Timoshenko theories
Table 8.1: Natural frequencies of a stepped bar in longitudinal vibration (results from the
conventional classical theory are shown in the parenthesis in column 2).
shown in Figure 8.8 by black solid lines which are in broad agreement with the ones re-
ported in [82]. The mode shapes shown by red dashed lines are those computed using
the classical Bernoulli-Euler theory. Clearly, the first three mode shapes have undergone
very little change as the result of using the present theory as opposed to the classical
Bernoulli-Euler theory, but the fourth mode being a higher order mode has turned out to
be significantly different, as expected. The exact reason for the discrepancies in the first
and third natural frequencies were unable to be pinpointed when the results are compared
with those of [82], but it should be recognised that the series solution approach used in
[82] is different from the dynamic stiffness methodology developed in this chapter. It is
to be noted that both the Rayleigh-Love and the classical Bernoulli-Euler theories give
almost the same results for the fundamental natural frequency, but the differences in the
second, third and fourth frequencies are 7%, 4% and 21% respectively. Understandably,
the classical Bernoulli-Euler theory overestimates the natural frequencies whereas the
more refined Rayleigh-Love theory which accounts for the added transverse and lateral
inertia of the bar yields lower values of the natural frequencies which is apparently con-
tradicted by the result for the fundamental natural frequency reported in [82].
8.5. Results and discussion 175
Figure 8.8: Natural frequencies and mode shapes of the three-stepped bar of Figure 8.7.
The final set of results was obtained using example-3 which is that of a plane frame
shown in Figure 8.9. Each element of the frame has the same uniform geometrical, cross
sectional and material properties and the data used in the analysis are as follows: EI =
4.0×106 Nm2 , EA = 8.0×108 N, kAG = 2.0×108 N, ρA = 30 kg/m, ρI p = 0.157 kgm, ν =
1/3, k = 2/3.
176 Chapter 8. Combined Rayleigh-Love and Timoshenko theories
Figure 8.9: A plane frame for free vibration analysis using Rayleigh-Love and Timo-
shenko theories.
A wide range of the natural frequencies of the frame was computed using the present
theory as well as the classical Bernoulli-Euler theory. Apart from the computation of
the first five natural frequencies which were sequentially chosen, the higher order natural
frequencies were sparingly and sparsely chosen so as to cover the order of the natural
frequencies between 50th and 400th . The results are shown in Table 8.2. Clearly, higher
the order of the frequency, higher the incurred error due to using the classical Bernoulli-
Euler theory. The first five natural frequencies of the frame are virtually unaltered. As
expected, the classical Bernoulli-Euler theory overestimates the natural frequencies.
One of the potential application areas of the theory developed in this chapter is the Sta-
tistical Energy Analysis (SEA) for which accurate natural frequency predictions in the
low, medium and high frequency range are essential. To this end, the uncompromis-
ing accuracy of the dynamic stiffness method developed in this chapter by applying the
Rayleigh-Love and Timoshenko theories is further demonstrated by computing the num-
ber of natural frequencies of the frame (see Figure 8.10) which lies within the frequency
ranges of 0 < fi ≤ 1kHz, 0 < fi ≤ 2kHz, 0 < fi ≤ 3kHz and up to 0 < fi ≤ 10kHz
which cover low, medium and high frequency bands. Figure 8.10 shows the frequency
distribution, i.e. the modal density of the frame. It will be difficult to obtain these results
with such accuracy using conventional finite element method.
8.5. Results and discussion 177
8.6 Conclusions
Starting from the derivations of the governing differential equations of motion in free
vibration, the dynamic stiffness matrix of a beam using both the Rayleigh-Love and Tim-
oshenko theories has been developed. With the help of the Wittrick-Williams algorithm
as solution technique, the theory is applied to investigate the free vibration behaviour of
a uniform Rayleigh-Love bar, a stepped Rayleigh-Love bar, and a framework for which
the modal density distribution is presented by capturing its natural frequencies in the low,
medium and high frequency range. Some representative mode shapes of the stepped bar
are also illustrated. The theory developed is particularly helpful when carrying out high
frequency free vibration analysis of skeletal structures. A potential application of the re-
search described in this chapter falls within the area of statistical energy analysis for which
the knowledge of modal density distributions in the high frequency range is essential.
Chapter 9
formation. No one appears to have made any attempt to solve the free vibration problem
of axial-bending coupled beams by using the dynamic stiffness method (DSM) which is
well known for its accuracy and computational efficiency [32, 93, 94]. The purpose of this
chapter is to fill this gap in the literature. Starting from the derivation of the governing
differential equations of motion, the dynamic stiffness matrix of an axial-bending cou-
pled beam is developed. The resulting dynamic stiffness matrix is exploited through the
application of the Wittrick-Williams algorithm [23] as solution technique to compute the
natural frequencies and mode shapes of an illustrative example of axial-bending coupled
beam for various boundary conditions. The results are contrasted with those obtained
from the classical beam theory which ignores the axial-bending coupling effects.
sections is shown in Figure 9.2. The method developed below is focused on axial-bending
coupling only and other forms of couplings arising from shear, torsion and warping effects
are not included in the theory. The principal assumptions made are those of linear small
deflection theory and also that the cross section of the beam is singly symmetric. In the
formulation, the contributions of shear stress and transverse normal stress are assumed to
be small and hence neglected in the analysis. The beam deforms only in one plane and
any form of non-linearity arising from large deflections and/or geometric stiffness due to
the presence of any axial load is not considered in this chapter, but interested readers are
referred to [83–85] which provide useful information on the development of non-linear
beam models.
Figure 9.1: Coordinate system and notation for an axial-bending coupled beam. Gc :
Centroid, E s : Shear centre.
Figure 9.2: Samples of beam cross sections with non-coincident centroid and shear cen-
tre.
If v, w and θ are axial displacement, bending displacement and bending rotation of a point
at a distance y from the origin and at a height z from the elastic axis, i.e. the point (y, z) in
the coordinate system (Figure 9.1), one can write
′
v = v0 − zw0 , w = w0 (9.1)
9.3. Theoretical development 183
where v0 and w0 are the corresponding displacement components of the point (y, 0) on the
Y-axis (i.e. the elastic axis) and a prime denotes differentiation with respect to y. (Note
that v and v0 represent the displacement components and must not be confused with the
velocity.)
Using linear, small deflection elasticity theory, the expression for the strain εy in the Y-
direction can be expressed as
′ ′′
εy = v0 − zw0 (9.2)
1 L
Z Z
U= Eε2y dAdy (9.3)
2 0 A
where E is the Young’s modulus of the beam material and the integrations are carried out
over the beam cross-sectional area A and length L. Note that the effect of shear deforma-
tion is assumed to be small and hence, neglected in the analysis.
Substituting εy from Equation 9.2 into Equation 9.3, and integrating over the beam cross-
section, we obtain
1 L
Z
EA[v0 ]2 − 2EAzα v0 w0 + EIe [w0 ]2 dy
′ ′ ′′ ′′
U= (9.4)
2 0
where A and Ie are the area of cross-section and second moment of area about the elastic
axis so that EA and EIe are the extensional and bending stiffnesses of the beam, respec-
tively.
The kinetic energy of the beam is given by
1 L
Z Z
ρ (v̇)2 + (ẇ)2 dy
n o
T= (9.5)
2 0 A
where ρ is the density of the beam material and an over dot represents differentiation with
respect to time t.
Equation 9.5 with the help of Equation 9.1 becomes
1 L
Z ′ 2
ρA(v̇0 )2 − 2ρAzα v̇0 ẇ0 + ρIe ẇ0 + ρA(ẇ0 )2 dy
′
T= (9.6)
2 0
where t1 and t2 are the time interval in the dynamic trajectory, and δ is the usual variational
operator.
The governing differential equations of motion for the axial-bending coupled beam and
the associated boundary condition in free vibration can now be derived by substituting
184 Chapter 9. Coupled axial-bending DSM for beam elements
the potential (U) and kinetic (T) energy expressions of Equations 9.4 and 9.6 into Equa-
tion 9.7, using the δ operator, integrating by parts and then collecting terms. In an earlier
publication, the entire procedure to generate the governing differential equations of mo-
tion and natural boundary conditions for bar or beam type structures was automated by
Banerjee et al [22] by applying symbolic computation. In this way, the governing differ-
ential equations of motion of the axial-bending coupled beam and the associated natural
boundary conditions are obtained as follows. Governing differential equations:
′′ ′′′ ′
EAv0 − EAzα w0 − ρAv̈0 + ρAzα ẅ0 = 0 (9.8)
′′′′ ′′′ ′ ′′
EIe w0 − EAzα v0 + ρAẅ0 + ρAzα v̈0 − ρIe ẅ0 = 0 (9.9)
Bending moment:
′′ ′
M = −EIe w0 + EAzα v0 (9.11)
Shear force:
′′′ ′′ ′
S = EIe w0 − EAzα v0 + ρAzα v̈0 − ρIe ẅ0 (9.12)
Assuming harmonic oscillation with circular or angular frequency ω rad/s, one can write
where V and W are the amplitudes of the axial and bending displacements, respectively.
Substituting Equation 9.13 into Equations 9.8 and 9.9 and introducing the non-dimensional
length ξ = y/L and the differential operator D = dξd , yield the following ordinary differen-
tial equations in V and W
ω2 ρAzα
!
2 EA 2 EAzα 3
ω ρA + 2 D V − D+ 3 D W =0 (9.14)
L L L
where
H = V or W (9.17)
9.3. Theoretical development 185
and
n o n o
a2 + b2 r02 − 2µ2 b2 a2 r02 − µ2 − 1 a2 b2
a1 = b2 µ2
; b1 = b2 µ2
; c1 = b2 µ2
(9.18)
1− a2
1− a2
1− a2
with
2 ρAω2 L2 2 ρAω2 L4 2 IE a2 2 z2α
a = ; b = ; r0 = = ; µ = 2 (9.19)
EA EIe AL2 b2 L
By assuming the solution in the form H = eλξ where λ is a constant, yet to be determined,
the characteristic or auxiliary equation of the differential equation, Equation 9.16 can be
expressed as
λ6 + a1 λ4 + b1 λ2 − c1 = 0 (9.20)
The polynomial equation, Equation 9.20 can be reduced to a cubic and solved analytically
using standard procedure [107]. By taking the square root of the three roots of the cubic,
which could be real or complex, the six roots r j ( j = 1, 2, · · · , 6) of the characteris-
tic or auxiliary equation Equation 9.20 can be computed leading to the solutions of the
differential equation, Equation 9.16 as:
6
X 6
X
V (ξ) = R je ;
r jξ
W (ξ) = Q j er j ξ (9.21)
j=1 j=1
where R j and Q j ( j = 1, 2, · · · , 6) are two sets of different constants which can be related
to each other by using Equations 9.14 and 9.15. The relationship between R j and Q j is
obtained as:
Q j = α jR j (9.22)
where
µb2 r j a2 + r2j
αj = n o (9.23)
a2 r4j − b2 1 − r02 r2j
The constant vectors Q and R can be written as:
The amplitudes of the axial force (F), shear force (S) and bending moment (M) are ob-
tained in terms of R j using Equations 9.10-9.13 and Equations 9.18 and 9.19 as
6
d2 W
!
EA dV EA X
F=− −µ 2 =− r j 1 − µα j r j R j er j ξ (9.26)
L dξ dξ L j=1
186 Chapter 9. Coupled axial-bending DSM for beam elements
EI d3 W µb2 dV
!
2 2 dW 2
S = 3 + b r0 − 2 + µb V
L dξ3 dξ a dξ
6
r 2
(9.27)
EI X
2 2 2
j
b2 +
R er j ξ
= 3 α r r + b r − µ
j j j 0 2 j
L
j=1
r
0
6 ! 6
EI d2 W µ dV
!
EI X µ X
M=− 2 2
− 2 =− 2 r j α jr j − 2 R j er j ξ (9.28)
L dξ r0 dξ L j=1 r0 j=1
At ξ = 0 : V = V1 ; W = W1 ; Θ = Θ1 ; F = F1 ; S = S 1 ; and M = M1 (9.29)
Figure 9.3: Sign convention for positive axial force F, shear force S and bending moment
M.
The displacement vector δ and the force vector P of the beam connecting the ends 1 and
2, see Figure 9.4, can be expressed as:
The relationship between the displacement δ and the constant vector R is now obtained
using Equations 9.21-9.23, 9.25 and Equations 9.29-9.30 to give
δ=BR (9.32)
9.3. Theoretical development 187
Figure 9.4: Boundary condition for displacements and forces for an axial-bending cou-
pled beam.
where
1 1 1 1 1 1
α1 α2 α3 α4 α5 α6
r1 α1 /L r2 α2 /L r3 α3 /L r4 α4 /L r5 α5 /L r6 α6 /L
B = r r r r r
(9.33)
e 1
e 2
e 3
e 4
e 5
er6
α1 er1 α2 er2 α3 er3 α4 er4 α5 er5 α6 er6
r1 α1 er1 /L r2 α2 er2 /L r3 α3 er3 /L r4 α4 er4 /L r5 α5 er5 /L r6 α6 er6 /L
Similarly, the relationship between the force vector P and the constant vector R is obtained
using Equations 9.26-9.30 to give
P=AR (9.34)
where elements of each row of the A matrix are indicated by the first of the two subscripts
as given below with j representing the column number.
2
EA n o EI
2 2 2
r j
b2 +
r j 1 − µα j r j ; ;
a1 j = − a2 j = 3 α r r + b r − µ
j j j 0 2
L L
r0
( )
EI µ EA n o
a3 j = − 2 r j (α j r j − 2 ) ; a4 j = r j 1 − µα j r j er j ; (9.35)
L r0 L
2
r
( )
EI
2 2 2
2 j EI µ
e ; a6 j = 2 r j (α j r j − 2 ) er j ;
r
a5 j = − 3 α r r + b r0 − µ b + 2 j
L j j j r0 L r0
By eliminating the constant vector R from Equations 9.32 and 9.34, P and δ can now
be related to give the dynamic stiffness matrix relationship of the axial-bending coupled
beam as
P=Kδ (9.36)
188 Chapter 9. Coupled axial-bending DSM for beam elements
where
K = A B−1 (9.37)
is the required frequency-dependent dynamic stiffness matrix. It should be noted that the
resulting dynamic stiffness matrix K of Equation 9.37 will be always symmetric and real
[95] with imaginary part of each element being zero although the matrices A and B are
complex and asymmetric. The expanded dynamic stiffness matrix giving the relationship
between the amplitudes of the forces to those of the displacements can be expressed in
the following way.
F1 k11 k12 k13 k14 k15 k16 V1
S 1 k12 k22 k23 k24 k25 k26 W1
M1 k13 k23 k33 k34 k35 k36 Θ1
F = k
V (9.38)
k k k k k
2 14 24 34 44 45 46 2
S k
2 15 k25 k35 k45 k55 k56 W2
M2 k16 k26 k36 k46 k56 k66 Θ2
or
P1 K11 K12 ∆1
= (9.39)
P2 K21 K22 ∆2
where K11 , K12 , K21 and K22 are all submatrices of order 3×3 each and P1 and P2 are
force vectors of node 1 (left-hand end) and node 2 (right-hand end) and ∆1 and ∆2 are the
corresponding displacement vectors, respectively.
The above frequency dependent dynamic stiffness matrix K can now be used to compute
the natural frequencies and mode shapes of either an individual axial-bending coupled
beam, or an assembly of axial-bending coupled beams for different boundary conditions.
A reliable and accurate method of solving the eigenvalue problem is to apply the Wittrick-
Williams algorithm [23] which is well suited for the DSM applications. The algorithm
uses the Sturm sequence property of the dynamic stiffness matrix and ensures that no
natural frequencies of the structure analysed are missed.
1 m. The material used is steel with Young’s modulus E = 200 GPa and density ρ =7850
kg/m3 . When preparing data, the stiffness, mass and other geometrical properties of the
channel section are calculated as follows:
(i) Axial stiffness (EA) = 2.9×108 N, (ii) Bending stiffness (EIe ) = 1.965×106 Nm2 , (iii)
Mass per unit length (ρA) = 11.383 kg/m, (iv) Rotatory inertia per unit length (ρIe ) =
0.07713 kgm and (iv) Elastic axis off-set from the mass axis (zα ) = 0.075616 m. Based
on these data and referring to Figure 9.2(a), it can be ascertained that the elastic axis and
mass axis are respectively 0.090356 m and 0.01474 m below the mid-length of the sides
with height h. It should be noted that the nature of the coupling in inertial only and the
stiffness coupling is ignored in the model.
The first five natural frequencies of the beam using the present theory are shown in Table
9.1 together with the corresponding results computed by using the classical beam theory
for C-F, P-P, C-P and C-C boundary conditions. It should be noted that when computing
results to simulate the classical beam theory as a degenerate case of the present theory, the
elastic axis off-set from the mass axis (zα ) was set to a small number close to zero (typi-
cally of the order of 10−6 ) in order to avoid numerical overflow and the bending stiffness
for the input data was recalculated about the centroidal axis to give EIg = 3.0689×105
Nm2 . The results for the classical beam theory case were further checked up to the ac-
curacy quoted in Table 9.1 by using the results quoted in standard text (for example, see
page 278 of [96]). These were also checked with the help of the computer program pub-
lished by Williams and Howson [60] who used the traditional uncoupled classical beam
theory when applying the dynamic stiffness method. The discrepancies in the result for
the five natural frequencies between the classical beam theory and the present theory are
shown in Table 9.1 for each set of the boundary conditions. The percentage difference is
given relative to the present theory. Clearly, unacceptably large errors can incur due to
using the classical beam theory as shown in the table.
In order to confirm the correctness of the natural frequencies shown in Table 9.1 additional
checks has been performed due to the absence of comparative results in the literature.
This was carried out by using the individual elements of the dynamic stiffness matrix of
Equations 9.38 and 9.39 and imposing the necessary boundary conditions. For instance,
the determinant formed by the matrix K11 (or K22 ) of Equation 9.38 was set to zero to
represent the cantilever boundary condition (C-F) and the determinant value was then
computed for a wide range of frequencies. The zeroes of the determinant established the
natural frequencies of the cantilever beam as an alternative method without resorting to
the Wittrick-Williams algorithm [23] as solution technique. The latter of course, is robust
and has a much wider applicability. For instance, a non-uniform axial-bending coupled
beam or a framework consisting of several axial-bending coupled beams can be modelled
as an assembly of many uniform axial-bending coupled beams to form the overall dy-
190 Chapter 9. Coupled axial-bending DSM for beam elements
namic stiffness matrix of the final structure to which the Wittrick-Williams algorithm can
be applied to compute the natural frequencies and mode shapes in a straightforward man-
ner. This cannot be easily accomplished by using the determinant method. An illustrative
example of the determinant plot for the above C-F beam is shown in Figure 9.5 where
its first two natural frequencies are identified at 575.85 rad/s (91.65 Hz) and 3557.2 rad/s
(566.15 Hz), respectively by tracking the zeroes of the determinant |K11 | when it crosses
the horizontal axis representing the frequency. These two natural frequencies agreed com-
pletely with the C-F results reported in column 2 of Table 9.1. These results were further
checked by imposing the boundary conditions in Equations 9.32-9.33 and Equations 9.34-
9.35 and tracking the zeroes of the 6×6 determinant formed by the first three rows of B
matrix and the last three rows of A matrix which together essentially represent the can-
tilever boundary conditions with the left hand end built-in and the right hand end free.
Further checks were performed for other boundary conditions, the details of which are
not reported here for brevity.
9.4. Results and discussion
Table 9.1: Natural frequencies of a channel section beam for different boundary conditions using present theory and classical beam theory.
191
192 Chapter 9. Coupled axial-bending DSM for beam elements
Figure 9.5: Determinant plot of K11 to locate the first two natural frequencies of the
cantilever channel section beam.
The mode shapes corresponding to the natural frequencies of Table 9.1 were computed
using the present theory and the classical beam theory. Results for the C-F, P-P, C-P
and C-C boundary conditions are illustrated in Figures 9.6-9.9, respectively. Clearly the
modes generated by the classical beam theory are uncoupled for all cases, as expected.
By contrast, the present theory shows substantial coupling between the axial and bending
deformation in most of the cases. For the C-F case, the first mode is predominantly bend-
ing with very little axial deformation whereas the second, fourth and the fifth modes are
heavily coupled. However, the third mode is essentially an axial mode with no bending
displacement present. It should be noted that pure axial mode is possible for the C-F
boundary condition. This is because the axial inertial forces can be possibly balanced by
the elastic forces arising from the axial deformations only without involving any bending
9.4. Results and discussion 193
deformation. For this C-F boundary conditions, a direct comparison between the results
computed by the present theory and the classical beam theory indicates that a relatively
small change in the natural frequency due to the application of the two theories, can cause
substantial changes in the mode shapes. In particular, the fourth and fifth modes shown
in Figure 9.6 are heavily coupled in axial and bending motions when using the present
theory, but the percentage differences in the corresponding natural frequencies for these
two modes when compared with the classical uncoupled classical beam theory are only
around 4% and 7%, respectively. Similar observations were made for other boundary
conditions. Interestingly, in an earlier investigation on the free vibration behaviour of
twisted beams, it was shown that even a substantial amount of twist (up to 30 degrees)
caused very little difference to the natural frequencies, but significant changes to the mode
shapes (see Figure 4 of [97]).
Figure 9.6: The first five natural frequencies and mode shapes for the channel-section
beam using present theory and classical beam theory for cantilever boundary condition.
The results shown in Table 9.1 clearly indicate that the errors incurred in the natural fre-
quencies for the pinned-pinned (P-P) boundary condition of the channel section beam can
be very large when applying the classical beam theory as opposed to the present theory.
194 Chapter 9. Coupled axial-bending DSM for beam elements
Figure 9.7: The first five natural frequencies and mode shapes for the channel-section
beam using present theory and classical beam theory for pinned-pinned boundary condi-
tion.
The next set of results comprising the natural frequencies for the clamped-pinned (C-P)
boundary conditions using the present theory and the classical beam theory are shown in
columns 8 and 9 of Table 9.1, respectively with the percentage differences given in the
10th column. Clearly, the errors due to using the classical beam theory are significant, but
not as great as was in the case of P-P boundary condition. The mode shapes for the five
natural frequencies for this C-P case are shown in Figure 9.8 which illustrate the pres-
ence of substantial coupling between the axial and bending deformation when using the
present theory whereas there is no coupling present when using classical beam theory, as
expected. Pure axial mode arising from the classical beam theory (third mode of Figure
9.8) is essentially a coupled mode which the present theory is capable to detect. As ex-
plained in the previous paragraph, the presence of a pinned support at any end of the beam
at its shear centre prevents the occurrence of any pure axial mode when using the present
theory.
196 Chapter 9. Coupled axial-bending DSM for beam elements
Figure 9.8: The first five natural frequencies and mode shapes for the channel-section
beam using present theory and classical beam theory for clamped-pinned boundary con-
dition.
The final set of results was obtained for the clamped-clamped (C-C) boundary conditions
for which the first five natural frequencies have already been shown in columns 11 and 12
of Table 9.1 using both the present and classical beam theories. The percentage difference
relative to the present theory is shown in the 13th column of the table. The correspond-
ing mode shapes are shown in Figure 9.9 which illustrate substantial coupling between
the axial and bending deformations in the first, second, fourth and fifth mode captured
by the present theory. The classical beam theory naturally fails to predict such couplings
as expected, which is demonstrated by the mode shapes shown in Figure 9.9. The third
mode is however, a pure axial mode captured by both theories. It should be recognised
that the C-C boundary condition can produce a pure axial mode as was the case with C-F
boundary condition. This is because the clamped support at the ends of the beam allows
neither bending displacement nor bending rotation which makes it possible for the beam
to freely vibrate in the axial direction only as a distinctive case.
9.4. Results and discussion 197
Figure 9.9: The first five natural frequencies and mode shapes for the channel-section
beam using present theory and classical beam theory for clamped-clamped boundary con-
dition.
free-free (F-F) boundary condition of the example problem (channel section beam) given
in section 9.4 would be much more convincing than other boundary conditions when val-
idating the theory and results. This is because the F-F boundary condition will make use
of all the dynamic stiffness terms unlike other boundary conditions for which some of the
stiffness terms will not show up in the calculation as they will be constrained due to the
support conditions.
Table 9.2: Natural frequencies of a channel section beam for free-free boundary condi-
tions using present theory and CUF.
9.5 Conclusions
Axial-bending coupled dynamic stiffness matrix for a beam with non-coincident centroid
and shear centre has been developed by deriving the governing differential equations us-
ing Hamilton’s principle, then solving the equations and finally imposing the boundary
conditions. The resulting dynamic stiffness matrix is applied with particular reference
to the Wittrick-Williams algorithm to solve the free vibration problems of an illustrative
example with substantial coupling between the axial and bending deformation. The nat-
ural frequencies and mode shapes for different boundary conditions are compared and
contrasted with those obtained from the classical beam theory. It has been shown that
the classical beam theory may give unacceptably large errors when investigating the free
vibration characteristics of axial-bending coupled beams. The theory and results are also
validated using Carrera unified formulation (CUF). The developed dynamic stiffness ma-
trix can be applied to complex structures, including frameworks.
Chapter 10
• Both metallic and composite wing has been analysed. Bending–torsional material cou-
pling, which is not possible in metallic wing, can be varied on the composite wing to
achieve desired aeroelastic effect. Also, by varying ply orientations and by using double
cells improved results can be achieved.
• The geometric features such as cut-outs, sweep, taper ratio, flexible ribs in the 3D of
the wing box have significant effect on the GJ results, especially taper along the wing
planform contribute more to the torsional stiffness rather than the rest such as manholes,
rigidity of ribs and sweep angle.
• It can be ascertained from the computed results that for sailplanes, the variation of EI
does not affect the flutter speed to any appreciable extent while the flutter speed varies
proportional to the change in GJ. For light aircraft trainers, the flutter speed varies pro-
portional to both EI and GJ. For transport airliner, the flutter speed varies inversely pro-
portional to change in EI and directly proportional to change in GJ. Overall GJ variation
affects flutter speed more than EI variation.
• The literature on the free vibration of frameworks containing FGBs is virtually non-
existent and chapter 6 addresses this by analysing a portal frame.
• The natural frequencies of intact beam for cantilever boundary condition was available
in the existing literature and in chapter 7 natural frequencies of intact beam for simply sup-
ported, clamped-pinned and clamped-clamped boundary conditions have been obtained.
• A potential application of the research described in chapter 8 falls within the area of
statistical energy analysis for which the knowledge of modal density distributions in the
high frequency range is essential.
• It has been shown in chapter 9 that the classical beam theory may give unacceptably
large errors when investigating the free vibration characteristics of axial-bending coupled
beams.
10 1 2 Line 1
• Number of laminate parts
• 1 – SI units
• 2– Circumferential asymmetric configuration
20 20 20 20 20 20 20 20 Line 2
• 20 – Number of Plies
204 Appendix A. Programs used in Section A
45.0 0.0 -45.0 90.0 -45.0 -45.0 90.0 -45.0 0.0 45.0 Line 4
• Ply orientation
2 23 1 2 Line 1
• 2 – Double cell
• 23 – Number of laminates
• 1 – SI units
• 2 – Circumferential asymmetric configuration
42 42 42 42 42 42 42 Line 2
• 42 – Number of plies
45.0 0.0 45.0 90.0 0.0 0.0 90.0 45.0 0.0 45.0 Line 4
• Ply orientation
1 18 Line 10
• The upper and lower node number connecting the mid - wall
206 Appendix A. Programs used in Section A
1 2 3 4 5 6 Line 2
• Mode numbers
23.5 Line 3
• Sweep angle in degrees
13 Line 4
• Number of beam elements used to idealise the wing
13 1 Line 5
• 13 – Non – uniform wing, number of beam elements
• 1 – SI unit
1 Line 7
• Number of nodes with lumped mass
1 Line 10
• Modal analysis as well as flutter analysis
1 1 Line 13
• 1 - Data group given in lines 11 and 12 are used in flutter analysis
• 1 – user specified accuracy needed for flutter speed only
2 Line 15
• Print flutter speed, flutter frequency, natural frequency, normal modes, real and
imaginary parts of the flutter determinants
1 2 3 4 5 6 Line 2
• Mode numbers
23.5 Line 3
• Sweep angle in degrees
13 Line 4
• Number of beam elements used to idealise the wing
13 1 Line 5
• 13 – Non – uniform wing, number of beam elements
• 1 – SI unit
1 Line 7
208 Appendix A. Programs used in Section A
1 Line 10
1 1 Line 13
2 Line 15
• Print flutter speed, flutter frequency, natural frequency, normal modes, real and
imaginary parts of the flutter determinants
A.5. BIGCALFUN input (example data file) 209
-1 -2 -3 -4 -5 -6 Line 2
• Mode numbers
1 2 1 Line 3
• Connection list, showing the topology of the structure
2.36E07 5.44E08 7.29E+05 27.26 1.68 0.0001 1.0 0.0001 2.10E+06 0.04
Line 4
• EIZZ , AA, GJ, m/L, Ialpha , Zalpha , K, Xre f , EIXX , Xalpha
15 1 0 15 2 0 15 3 0 15 4 0 15 5 0 15 6 0 Line 7
• Degrees of freedom affected
0000 Line 10
• Tail data
210 Appendix A. Programs used in Section A
0000 Line 11
• Fin data
0 0 0 0 1 1 0 Line 13
• IFLAG, IFLAG2, IFLAG3, NZI, IDSL, NUNIT, NCM
1 2 2 1 0 Line 14
• IDWU, IQU2, IPOUT, IDVG, IQS
Table B.1: Laminate layup and stacking sequences for sections 1-19
LAMINATE
SECTION
LAYUP
SECTIONS
1-3 UPPER SKIN [45/45/0/45/0/-45/90/-45/90/90]s
LOWER SKIN [45/45/0/-45/0/45/90/-45/90/90]s
FRONT SPAR [45/0/45/45/-45/0/-45/90/90/-45]s
REAR SPAR [45/45/0/45/0/-45/90/90-45/90]s
SECTIONS
4-7 UPPER SKIN [45/45/0/45/0/0/-45/-45/90/-45/90/90]s
LOWER SKIN [45/45/0/45/-45/0/45/0/90/45/90/90]s
FRONT SPAR [45/0/45/45/-45/0/0/-45/90/90/-45/90]s
REAR SPAR [45/45/0/45/0/-45/0/-45/90/90/-45/90]s
SECTIONS
8-11 UPPER SKIN [45/45/0/45/45/0/0/-45/-45/90/0/-45/90/90]s
LOWER SKIN [45/45/0/45/-45/0/45/0/0/90/-45/-45/90/90]s
FRONT SPAR [45/0/45/45/45/-45/0/0/-45/90/90/90/-45/90]s
REAR SPAR [45/45/0/45/45/0/-45/0/-45/90/90/90-45/90]s
SECTIONS
12-15 UPPER SKIN [45/45/0/45/45/0/0/-45/-45/90/0/-45/90/90/90/90/-45/90]s
LOWER SKIN [45/45/0/45/45/-45/0/45/0/0/90/-45/-45/90/90/90/-45]s
FRONT SPAR [45/0/45/45/45/-45/45/0/0/0/-45/90/90/90/-45/90]s
REAR SPAR [45/45/0/45/45/0/45/0/-45/0/-45/-45/90/90/90/-45/90]s
SECTIONS
16-19 UPPER SKIN [45/45/45/0/45/45/0/0/-45-45/90/0/-45/-45-90/90/90/-45/90/90/90/-45/90]s
LOWER SKIN [45/45/45/0/0/45/45/-45/0/45/0/0/90/-45/45/90/90/90-45]s
FRONT SPAR [45/0/45/45/45/-45/45/0/0/0/-45/90/-45/90/-45/90/90/-45/90]s
REAR SPAR [45/45/45/0/45/45/0/45/0/0/-45/0/-45/-45/90/90/90/-45/90]s
213
Table B.2: Laminate layup and stacking sequences for sections 20-26.
Figure C.1: Case study 1 showing principal dimension of the cross section.
C.2. Case study of uniform wing box – refined mesh 215
Figure C.2: Case study 2 showing principal dimension of the cross section.
Figure C.3: Case study 3 showing principal dimension of the cross section.
216 Appendix C. Geometric representation of parametric case studies
Figure C.4: Case study 4 showing principal dimension of the cross section.
Figure C.5: Case study 5 showing principal dimension of the cross section.
C.6. Case study of trapezoidal wing box – symmetry taper 217
Figure C.6: Case study 6 showing principal dimension of the cross section.
Figure C.7: Case study 7 showing principal dimension of the cross section.
218 Appendix C. Geometric representation of parametric case studies
Figure C.8: Case study 8 showing principal dimension of the cross section.
C.9 Case study of wing box model with and without the
manhole
Figure C.9: Illustrations of case study 9a (Box 1), 9b (Box 2), 9c (Box 3), 9d (Box 4).
C.10. Case study of composite wing box showing section 16 219
Figure C.10: Case study 10 showing principal dimension of the cross section.
Bibliography
[11] JR Banerjee. “The dynamic stiffness method: theory, practice and promise”. In:
Computational Technology Reviews 11.1 (2015), pp. 31–57.
[12] JR Banerjee. “Coupled bending–torsional dynamic stiffness matrix for beam ele-
ments”. In: International journal for numerical methods in engineering 28.6 (1989),
pp. 1283–1298.
[13] JR Banerjee. “A FORTRAN routine for computation of coupled bending-torsional
dynamic stiffness matrix of beam elements”. In: Advances in Engineering Software
and Workstations 13.1 (1991), pp. 17–24.
[14] JR Banerjee. “Use and capability of Calfun—A program for calculation of flutter
speed using normal modes”. In: Proceedings of the International AMSE Confer-
ence on Modelling and Simulation, Athens, Greece. 1984, pp. 121–131.
[15] JR Banerjee. “User’s guide to the computer program CALFUN (CALculation of
Flutter speed Using Normal modes)”. In: MEAD/AERO Report No. 164. Depart-
ment of Mechanical Engineering and Aeronautics, The City University London,
1989.
[16] JR Banerjee and FW Williams. “Free vibration of composite beams-an exact method
using symbolic computation”. In: Journal of Aircraft 32.3 (1995), pp. 636–642.
[17] JR Bannerjee and FW Williams. “Exact dynamic stiffness matrix for composite
Timoshenko beams with applications”. In: Journal of sound and vibration 194.4
(1996), pp. 573–585.
[18] Theodore Theodorsen. “General theory of aerodynamic instability and the mecha-
nism of flutter”. In: (1979).
[19] Samuel J Loring. Use of generalized coordinates in flutter analysis. Tech. rep. SAE
Technical Paper, 1944.
[20] Thomas Henry Gordon Megson. Aircraft structures for engineering students. Butterworth-
Heinemann, 2016.
[21] Liviu Librescu and Ohseop Song. Thin-walled composite beams: theory and appli-
cation. Vol. 131. Springer Science & Business Media, 2005.
[22] Frank L Matthews and Rees D Rawlings. Composite materials: engineering and
science. CRC press, 1999.
[23] JR Banerjee, Huijuan Su, and C Jayatunga. “A dynamic stiffness element for free
vibration analysis of composite beams and its application to aircraft wings”. In:
Computers & Structures 86.6 (2008), pp. 573–579.
[24] GA Georghiades and JR Banerjee. “Flutter prediction for composite wings using
parametric studies”. In: AIAA journal 35.4 (1997), pp. 746–748.
222 Bibliography
[25] JR Banerjee and FW Williams. “Free vibration of composite beams-an exact method
using symbolic computation”. In: Journal of Aircraft 32.3 (1995), pp. 636–642.
[26] Terrence A Weisshaar and RJ Ryan. “Control of aeroelastic instabilities through
stiffness cross-coupling”. In: Journal of Aircraft 23.2 (1986), pp. 148–155.
[27] Erian A Armanios and Ashraf M Badir. “Free vibration analysis of anisotropic
thin-walled closed-section beams”. In: AIAA journal 33.10 (1995), pp. 1905–1910.
[28] Ashraf Badir. “Analysis of two-cell composite beams”. In: 36th Structures, Struc-
tural Dynamics and Materials Conference. 1995, p. 1208.
[29] Marthinus C Van Schoor and Andreas H von Flotow. “Aeroelastic characteristics
of a highly flexible aircraft”. In: Journal of Aircraft 27.10 (1990), pp. 901–908.
[30] Deman Tang and Earl H Dowell. “Experimental and theoretical study on aeroelas-
tic response of high-aspect-ratio wings”. In: AIAA journal 39.8 (2001), pp. 1430–
1441.
[21] Anthony C Hearn. REDUCE user’s manual. Version 3.2. Rand Corporation, 1985.
[22] JR Banerjee et al. “Use of computer algebra in Hamiltonian calculations”. In: Ad-
vances in Engineering Software 39.6 (2008), pp. 521–525.
[23] W H Wittrick and FW Williams. “A general algorithm for computing natural fre-
quencies of elastic structures”. In: The Quarterly Journal of Mechanics and Applied
Mathematics 24.3 (1971), pp. 263–284.
[24] Diana Virginia Bambill, Carlos Adolfo Rossit, and Daniel Horacio Felix. “Free vi-
brations of stepped axially functionally graded Timoshenko beams”. In: Meccanica
50.4 (2015), pp. 1073–1087.
[25] Nuttawit Wattanasakulpong and Jarruwat Charoensuk. “Vibration characteristics of
stepped beams made of FGM using differential transformation method”. In: Mec-
canica 50.4 (2015), pp. 1089–1101.
[26] Metin Aydogdu and Vedat Taskin. “Free vibration analysis of functionally graded
beams with simply supported edges”. In: Materials & design 28.5 (2007), pp. 1651–
1656.
[27] Zheng Zhong and Tao Yu. “Analytical solution of a cantilever functionally graded
beam”. In: Composites Science and Technology 67.3-4 (2007), pp. 481–488.
[28] SA Sina, HM Navazi, and H Haddadpour. “An analytical method for free vibra-
tion analysis of functionally graded beams”. In: Materials & Design 30.3 (2009),
pp. 741–747.
[29] Shi-rong Li, Ze-qing Wan, and Jing-hua Zhang. “Free vibration of functionally
graded beams based on both classical and first-order shear deformation beam the-
ories”. In: Applied Mathematics and Mechanics 35.5 (2014), pp. 591–606.
[30] MA Eltaher, AE Alshorbagy, and FF Mahmoud. “Determination of neutral axis po-
sition and its effect on natural frequencies of functionally graded macro/nanobeams”.
In: Composite Structures 99 (2013), pp. 193–201.
[31] KS Al-Basyouni, Abdelouahed Tounsi, and SR Mahmoud. “Size dependent bend-
ing and vibration analysis of functionally graded micro beams based on modified
couple stress theory and neutral surface position”. In: Composite Structures 125
(2015), pp. 621–630.
[32] JR Banerjee. “Coupled bending–torsional dynamic stiffness matrix for beam ele-
ments”. In: International journal for numerical methods in engineering 28.6 (1989),
pp. 1283–1298.
References for Section B 225
[33] JR Banerjee. “Explicit analytical expressions for frequency equation and mode
shapes of composite beams”. In: International Journal of Solids and Structures
38.14 (2001), pp. 2415–2426.
[34] Louis A Pipes and Lawrence R Harvill. Applied mathematics for engineers and
physicists. Courier Corporation, 2014.
[35] CP Filipich and PAA Laura. “In-plane vibrations of portal frames with end sup-
ports elastically restrained against rotation and translation”. In: Journal of Sound
Vibration 117 (1987), pp. 467–473.
[36] İlhan Tatar. “Vibration characteristics of portal frames”. MA thesis. Izmir Institute
of Technology, 2013.
[37] Peter Gudmundson. “Eigenfrequency changes of structures due to cracks, notches
or other geometrical changes”. In: Journal of the Mechanics and Physics of Solids
30.5 (1982), pp. 339–353.
[38] Fr D Ju and ME Mimovich. “Experimental diagnosis of fracture damage in struc-
tures by the modal frequency method”. In: Journal of Vibration, Acoustics, Stress,
and Reliability in Design 110.4 (1988), pp. 456–463.
[39] Peter Gudmundson. “The dynamic behaviour of slender structures with cross-sectional
cracks”. In: Journal of the Mechanics and Physics of Solids 31.4 (1983), pp. 329–
345.
[40] TG Chondros, AD Dimarogonas, and J Yao. “A continuous cracked beam vibration
theory”. In: Journal of sound and vibration 215.1 (1998), pp. 17–34.
[41] S Christides and ADS Barr. “One-dimensional theory of cracked Bernoulli-Euler
beams”. In: International Journal of Mechanical Sciences 26.11-12 (1984), pp. 639–
648.
[42] AD Dimarogonas and CA Papadopoulos. “Vibration of cracked shafts in bending”.
In: Journal of sound and vibration 91.4 (1983), pp. 583–593.
[43] CA Papadopoulos and AD Dimarogonas. “Coupling of bending and torsional vi-
bration of a cracked Timoshenko shaft”. In: Ingenieur-Archiv 57.4 (1987), pp. 257–
266.
[44] CA Papadopoulos and AD Dimarogonas. “Coupled longitudinal and bending vi-
brations of a rotating shaft with an open crack”. In: Journal of sound and vibration
117.1 (1987), pp. 81–93.
[45] CA Papadopoulos and AD Dimarogonas. “Coupled longitudinal and bending vi-
brations of a cracked shaft”. In: Journal of vibration, acoustics, stress, and relia-
bility in design 110.1 (1988), pp. 1–8.
226 Bibliography
[46] Thomas M Tharp. “A finite element for edge-cracked beam columns”. In: Inter-
national Journal for Numerical Methods in Engineering 24.10 (1987), pp. 1941–
1950.
[47] N Miyazaki. “Application of line-spring model to dynamic stress intensity factor
analysis of pre-cracked bending specimen”. In: Engineering fracture mechanics
38.4-5 (1991), pp. 321–326.
[48] Robert Y Liang, Jialou Hu, and Fred Choy. “Theoretical study of crack-induced
eigenfrequency changes on beam structures”. In: Journal of Engineering Mechan-
ics 118.2 (1992), pp. 384–396.
[49] HP Lee and TY Ng. “Natural frequencies and modes for the flexural vibration of a
cracked beam”. In: Applied Acoustics 42.2 (1994), pp. 151–163.
[50] G Bamnios and A Trochides. “Dynamic behaviour of a cracked cantilever beam”.
In: Applied Acoustics 45.2 (1995), pp. 97–112.
[51] M Kisa, J Brandon, and M Topcu. “Free vibration analysis of cracked beams by a
combination of finite elements and component mode synthesis methods”. In: Com-
puters & structures 67.4 (1998), pp. 215–223.
[52] M Kisa and J Brandon. “The effects of closure of cracks on the dynamics of a
cracked cantilever beam”. In: Journal of sound and vibration 238.1 (2000), pp. 1–
18.
[53] Ding Yang Zheng and NJ Kessissoglou. “Free vibration analysis of a cracked
beam by finite element method”. In: Journal of Sound and vibration 273.3 (2004),
pp. 457–475.
[54] JA Loya, L Rubio, and J Fernández-Sáez. “Natural frequencies for bending vibra-
tions of Timoshenko cracked beams”. In: Journal of Sound and Vibration 290.3-5
(2006), pp. 640–653.
[55] E Viola, P Ricci, and MH Aliabadi. “Free vibration analysis of axially loaded
cracked Timoshenko beam structures using the dynamic stiffness method”. In:
Journal of Sound and Vibration 304.1-2 (2007), pp. 124–153.
[56] AS Bouboulas and NK Anifantis. “Formulation of cracked beam element for analy-
sis of fractured skeletal structures”. In: Engineering Structures 30.4 (2008), pp. 894–
901.
[57] Hiroyuki Okamura, Katsuhiko Watanabe, and Tachio Takano. “Deformation and
strength of cracked member under bending moment and axial force”. In: Engineer-
ing Fracture Mechanics 7.3 (1975), pp. 531–539.
References for Section B 227
[58] G Gounaris and A Dimarogonas. “A finite element of a cracked prismatic beam for
structural analysis”. In: Computers & Structures 28.3 (1988), pp. 309–313.
[59] N Papaeconomou and A Dimarogonas. “Vibration of cracked beams”. In: Compu-
tational Mechanics 5.2-3 (1989), pp. 88–94.
[60] FW Williams and WP Howson. “Compact computation of natural frequencies and
buckling loads for plane frames”. In: International Journal for Numerical Methods
in Engineering 11.7 (1977), pp. 1067–1081.
[61] ML Kikidis and CA Papadopoulos. “Slenderness ratio effect on cracked beam”. In:
Journal of Sound and Vibration 155.1 (1992), pp. 1–11.
[62] Andrew J Keane and WG Price. Statistical energy analysis: an overview, with ap-
plications in structural dynamics. Cambridge University Press, 1997.
[63] Richard H Lyon, Richard G DeJong, and Manfred Heckl. Theory and application
of statistical energy analysis. 1995.
[64] JC Wohlever and RJ Bernhard. “Mechanical energy flow models of rods and beams”.
In: Journal of sound and vibration 153.1 (1992), pp. 1–19.
[65] Y Lase, MN Ichchou, and L Jezequel. “Energy flow analysis of bars and beams:
theoretical formulations”. In: Journal of Sound and Vibration 192.1 (1996), pp. 281–
305.
[66] OM Bouthier and RJ Bernhard. “Simple models of energy flow in vibrating mem-
branes”. In: Journal of sound and vibration 182.1 (1995), pp. 129–147.
[67] OM Bouthier and RJ Bernhard. “Simple models of the energetics of transversely
vibrating plates”. In: Journal of Sound and Vibration 182.1 (1995), pp. 149–164.
[68] Bengt Å Åkesson. “PFVIBAT—a computer program for plane frame vibration
analysis by an exact method”. In: International Journal for Numerical Methods
in Engineering 10.6 (1976), pp. 1221–1231.
[69] WP Howson, JR Banerjee, and FW Williams. “Concise equations and program for
exact eigensolutions of plane frames including member shear”. In: Engineering
Software III. Springer, 1983, pp. 443–452.
[70] Ju-Bum Han et al. “Vibrational energy flow models for the Rayleigh–Love and
Rayleigh–Bishop rods”. In: Journal of Sound and Vibration 333.2 (2014), pp. 520–
540.
[71] Michael Shatalov et al. “Longitudinal vibration of isotropic solid rods: from clas-
sical to modern theories”. In: Advances in Computer Science and Engineering, M.
Schmidt, ed., InTech Open, Rijeka, Croatia (2011), pp. 187–214.
[72] John William Strutt and Baron Rayleigh. The theory of sound. Dover, 1945.
228 Bibliography
[73] Augustus Edward Hough Love. A treatise on the mathematical theory of elasticity.
Cambridge university press, 2013.
[74] VD Belov, SA Rybak, and BD Tartakovskii. “Propagation of vibrational energy
in absorbing structures”. In: Soviet Physics Acoustics-USSR 23.2 (1977), pp. 115–
119.
[75] DJ Nefske and SH Sung. “Power flow finite element analysis of dynamic systems:
basic theory and application to beams”. In: Journal of Vibration, Acoustics, Stress,
and Reliability in Design 111.1 (1989), pp. 94–100.
[76] Franklin Y Cheng. “Vibrations of Timoshenko beams and frameworks”. In: Journal
of the structural division 96.3 (1970), pp. 551–571.
[77] TM Wang and TA Kinsman. “Vibrations of frame structures according to the Tim-
oshenko theory”. In: Journal of Sound and Vibration 14.2 (1971), pp. 215–227.
[78] WP Howson and FW Williams. “Natural frequencies of frames with axially loaded
Timoshenko members”. In: Journal of Sound and Vibration 26.4 (1973), pp. 503–
515.
[79] Singiresu S Rao. Vibration of continuous systems. Vol. 464. Wiley Online Library,
2007.
[80] Mihai Valentin Predoi et al. “High frequency longitudinal damped vibrations of a
cylindrical ultrasonic transducer”. In: Shock and Vibration 2014 (2014).
[81] Walter C Hurty and Moshe F Rubinstein. “Dynamics of structures”. In: Prentice-
Hall Series in Engineering of the Physical Sciences, Englewood Cliffs: Prentice-
Hall, 1964 (1964).
[82] Igor Fedotov et al. “Analysis for an N-stepped Rayleigh bar with sections of com-
plex geometry”. In: Applied Mathematical Modelling 32.1 (2008), pp. 1–11.
[83] AS Yigit and AP Christoforou. “Coupled axial and transverse vibrations of oilwell
drillstrings”. In: Journal of sound and vibration 195.4 (1996), pp. 617–627.
[84] Seon Han and Haym Benaroya. “Coupled transverse and axial vibration of a com-
pliant tower-Comparison of linear and nonlinear models”. In: 41st Structures, Struc-
tural Dynamics, and Materials Conference and Exhibit. 2000, p. 1347.
[85] Marcelo A Trindade, Claudio Wolter, and Rubens Sampaio. “Karhunen–Loeve de-
composition of coupled axial/bending vibrations of beams subject to impacts”. In:
Journal of sound and vibration 279.3-5 (2005), pp. 1015–1036.
[86] Rubens Sampaio, Marcelo Tulio Piovan, and G Venero Lozano. “Coupled ax-
ial/torsional vibrations of drill-strings by means of non-linear model”. In: Mechan-
ics Research Communications 34.5-6 (2007), pp. 497–502.
References for Section B 229
[87] Jerry H Ginsberg. “Coupling of axial and transverse displacement fields in a straight
beam due to boundary conditions”. In: The Journal of the Acoustical Society of
America 126.3 (2009), pp. 1120–1124.
[88] Stefano Lenci and Giuseppe Rega. “Axial–transversal coupling in the free nonlin-
ear vibrations of Timoshenko beams with arbitrary slenderness and axial boundary
conditions”. In: Proceedings of the Royal Society A: Mathematical, Physical and
Engineering Sciences 472.2190 (2016), p. 20160057.
[89] Zhiyang Lei, Jinpeng Su, and Hongxing Hua. “Longitudinal and transverse cou-
pling dynamic properties of a Timoshenko beam with mass eccentricity”. In: Inter-
national Journal of Structural Stability and Dynamics 17.07 (2017), p. 1750077.
[90] Zhen Ni and Hongxing Hua. “Axial-bending coupled vibration analysis of an axially-
loaded stepped multi-layered beam with arbitrary boundary conditions”. In: Inter-
national Journal of Mechanical Sciences 138 (2018), pp. 187–198.
[91] Erasmo Carrera and Marco Petrolo. “On the effectiveness of higher-order terms in
refined beam theories”. In: Journal of Applied Mechanics 78.2 (2011), p. 021013.
[92] Phill-Seung Lee and Ghyslaine McClure. “A general three-dimensional L-section
beam finite element for elastoplastic large deformation analysis”. In: Computers &
structures 84.3-4 (2006), pp. 215–229.
[93] JR Banerjee. “Modal analysis of sailplane and transport aircraft wings using the
dynamic stiffness method”. In: Journal of Physics: Conference Series. Vol. 721. 1.
IOP Publishing. 2016, p. 012005.
[94] R Butler and JR Banerjee. “Optimum design of bending-torsion coupled beams
with frequency or aeroelastic constraints”. In: Computers & structures 60.5 (1996),
pp. 715–724.
[95] JR Banerjee. “Free vibration of sandwich beams using the dynamic stiffness method”.
In: Computers & structures 81.18-19 (2003), pp. 1915–1922.
[96] FS Tse, IE Morse, and RT Hinkle. “Mechanical Vibrations, Theory and Applica-
tions, Allyn and Becon”. In: Inc.(Pearson Publishers) (1978).
[97] JR Banerjee. “Free vibration analysis of a twisted beam using the dynamic stiff-
ness method”. In: International Journal of Solids and Structures 38.38-39 (2001),
pp. 6703–6722.
[98] RED Bishop, SM Cannon, and S Miao. “On coupled bending and torsional vibra-
tion of uniform beams”. In: Journal of sound and vibration 131.3 (1989), pp. 457–
464.