A THEORETICAL FORMULATION FOR FLUTTER

ANALYSIS OF A TYPICAL SUBSONIC AIRCRAFT

WING (SARAS) USING QUASI-STEADY
AERODYNAMIC THEORY
 ABSTRACT

A theoretical formulation for flutter analysis has been utilized to

develop a working method for determining flutter speed of a typical
subsonic aircraft wing. A Galerkin type of analysis has been used to derive
the matrix form of equations from the differential equations of motion of
the subsonic wing. Quasi-steady aerodynamic theory has been used to
model the aerodynamic forces.

A computer code in FORTRAN has been prepared for generation of

matrices while the eigen value analysis is performed through MATLAB.
The code is benchmarked through the flutter of a rectangular wing. The
results from the code agree reasonably with those obtained from the
industrial code NASTRAN.

The method is then extended to the flutter analysis of the actual

clean wing with no control surface effects. The tapered wing is modeled
as a stepped assembly of constant section beam elements. Results indicate
that the aircraft wing taken is very stiff and therefore is not flutter prone at
all in the subsonic regime. To simulate subsonic flutter conditions, a
hypothetically reduced stiffness analysis is performed.

In all the cases, the agreement of the results with those of

NASTRAN (that uses the Doublet Lattice Method, DLM) indicates the
validity of the present method of analysis using the quasi-steady
aerodynamic theory, the present work can be extended to study more
complicated cases of flutter in the aircraft wing with control surface effects
and the T-Tail assembly of aircraft which are expectedly quite prone to
subsonic flutter.
 ACKNOWLEDGEMENT

I express my sincere thanks to Dr. A.R.Upadhya, Director, Dr.

M.N.Sathyanarayana, Head, Technical Secretariat and Dr. S.Vishwanath, Head
of Structures Division, NAL, for providing the necessary facilities to carry out my
dissertation work at NAL.

I take this opportunity to express my gratitude and indebtedness to

Dr.Somenath Mukherjee, Scientist, NAL, Bangalore for suggesting this project
and keen interest through out the dissertation work. I learnt a lot of physics and
mathematics and also how to work independently from him.

I wish to express heartfelt gratitude to Dr. D.L.Prabhakara, Professor,

Department of Civil Engineering, Malnad College of Engineering, Hassan for his
excellent guidance during the course of this dissertation work.

I wish to thank Smt Manju and Smt Swarnalatha, Scientists, NAL,

Bangalore for their extended help during this project work.

I express my profound sense of gratitude to Dr.D.L.Manjunath, Professor

and Head, Department of Civil Engineering, Malnad College of Engineering,
Hassan for allowing me to do this dissertation work in NAL, Bangalore and for his
valuable suggestions and advice.

I am thankful to Dr. Karisiddappa, Professor, PG coordinator M.Tech

[CADS] Dept. of Civil Engineering, Dr. K.Munjunath, Assistant Professor and
Dr. M.S.Kalappa, Assistant Professor, Dept. of Civil Engineering Department,
Malnad College of Engineering, Hassan, for their constant support through the
course.

I am also grateful to my Sister and my Parents for their constant support.

 CHAPTER-1.

INTRODUCTION

1.1 Aeroelasticity

Aero elasticity is the study of the effect of the aerodynamic

forces on elastic bodies. If in the analysis of any structural dynamic
systems aerodynamic loading is included then the resulting dynamic
phenomenon may be classified as Aeroelastic.

The classical theory of elasticity deals with the stress and

deformation of an elastic body under prescribed external forces or
displacements. The external loading acting on the body is, in general,
independent of the deformation of the body. It is usually assumed that
the deformation is small and does not substantially affect the action of
external forces. In such a case, we often neglect the changes in
dimensions of the body and base our calculations on the initial shape.

The situation is different, however, in most problems of aero

elasticity. The aerodynamic forces depend critically on the attitude of
the body relative to the flow. The elastic deformation plays an
important role in determining the external loading itself. The
magnitude of the aerodynamic force is not known until the elastic
deformation is determined. In general, therefore, the external load is
not known until the problem is solved. Aero elastic phenomena have a
significant influence on the design of flight vehicles.

1.2 Aeroelastic flutter

One of the interesting problems in aeroelasticity is the stability

(or rather instability) of structure in wind. Since, for a given
configuration of the elastic body, the aerodynamic force increases
rapidly with the wind speed, while elastic stiffness is independent of
the wind, there may exist a critical wind speed at which the structure
becomes dynamically unstable. Such dynamic instability may cause
excessive oscillatory deformations that increase in amplitude
exponentially with time, and may lead to the destruction of the
structure.

A major problem is the flutter of structures such as airplanes or

suspension bridges, when small disturbances of an incidental nature
induce more or less violent oscillations. It is characterized by the
interplay of aerodynamic, elastic and inertia forces and is called a
problem of dynamic aeroelastic instability. The particular case of an
oscillation with zero frequency, in which in general the inertia force is
neglected, is called the steady state, or static aeroelastic instability.

There is a close relationship between the stability problems and

the response problems. Mathematically, most stability problems can
be described by a system of homogenous equations, which are
satisfied by a trivial solution of zero displacement or zero motion,
meaning that nothing happens at all. On the other hand, a response
problem is represented by a nonhomogenous system; i.e., the initial
conditions and the external forces are such as to cause the governing
equations to be nonhomogenous, and to admit a solution not
vanishing identically.

A response problem generally associates with a stability

problem. As an example, consider the response of an airplane wing to
atmospheric turbulences. We can formulate the problem of flutter by
asking the following questions: is there a critical speed of flight at
which the airplane structure becomes exceedingly sensitive to the
atmospheric turbulence; i.e., does there exist a speed at which the
structure may have a motion of finite amplitude, even in the limiting
case of an atmospheric turbulence of zero intensity? This is
equivalent to the following formulation, which is usually made in
flutter analysis: Is there a critical speed at which the aeroelastic
system becomes neutrally stable, at which motion of the structure is
possible without any external excitation?

Thus the response of the airplane structure to the atmospheric

turbulence and the flutter problem are linked together. When the
response of the structure to a finite disturbance is finite, the structure
is stable, and flutter will not occur. When the structure flutters at a
critical speed of flow, its response to a finite disturbance becomes
indefinite.
1.3 Literature review

The earliest study of flutter seems to have been made by

Lanchester [1], Bairstow and Fage [2] in 1916. In 1918, Blasius [3]
made some calculations after the failure of the lower wing of Albatross
D3 biplane. But the real development of the flutter analysis had to wait
for the development of Non-stationary airfoil theory by Kutta and
Joukowsky.

Glauret [4,5] published data on the force and moment acting on

a cylindrical body due to an arbitrary motion. In 1934, Theodorsens
[6] exact solution of a harmonically oscillating wing with a flap was
published.

The torsion flutter was first found by Glauret in 1929. It is

discussed in detail by Smilg [7]

Several types of single degree of freedom flutter involving

control surfaces at both subsonic and supersonic speeds have been
found [8,9], all requiring the fulfillment of certain special conditions
on the rotational axis locations, the reduced frequency and the mass
moment of inertia.

Pure bending flutter is possible for a cantilever swept wing if it

is heavy enough relative to the surrounding air and has a sufficiently
large sweep angle [10].
 The stability of more complicated motions can be determined by
calculating the energy input from the airstream. The bending torsion
case in an incompressible fluid has been calculated by J.H.Greidanus
and the energy coefficient in Bending-Torsion oscillations has been
given [11]

The use of Quasi-steady Aerodynamic theory for the flutter

analysis of the wings and excellent treatises in the field of
aeroelasticity are given by Y.C.Fung [12], E.H.Dowell [13,14],
L.Mirovitch [15] and others.

In the typical wing whose elastic axis (locus of shear centers)

and mass axis (locus of center of gravity) do not coincide, the nature of
oscillations is always coupled flexure-torsion. A vast literature exists
on the flexure-torsion problem of engineering structures. Evins [16]
has given comprehensive details about vibration fixture transducers
and instrumentation. Bisplinghoff and H.Ashley [17] has described the
elastic characteristics shape and inertial idealization.

A new method for determining mass and stiffness matrices from

modal test data is described by Alvin and Paterson [18]. This method
determines minimum order mass and stiffness matrices, which is used
to determine the optimum sensor location. Dugundji [19] examined
panel flutter and the rate of damping. The problem of two and three-
dimensional plate undergoing cyclic oscillations and aeroelastic
instability is investigated by Dowell [13,14].
 The behavior of tip loaded cantilever beam with an arbitrary
cross section using a power series solution technique for the out of
plane flexure and torsion case is discussed by Kosmataka [20]. This
includes a linear relation developed for locating shear center.
Abott [21] has suggested a technique for representing the shape of the
aerofoil through analytical relations.

The coupled flexure-torsion vibration response of beam under

deterministic and random load is investigated thoroughly by Eslimy
and Sobby [22] by use of normal mode method. The exact
determination of coupled flexure-torsion vibration characteristics of
uniform beam having single cross section symmetry is studied by
Dokumaci [23].

Literature in the context of excitation systems and design of

electromagnets is referred to by Chatterji [24] and Wilson [25]. Haisler
and Allen [26] and Kuhn [27] have described a procedure of
computing required centroidal and elastic properties including shear
flows. Mathematical formulations of flexure-torsion problem have
been broadly described by Mirovitch [15], Thomson [28], Hurty and
Rubinstein [29] and Y.C.Fung [12]. Talukedar, Kamle and Yadav [30]
have discussed an analytical method for flexure-torsion coupled
vibration of vehicles leading to aircraft application.

At present, subsonic flight is a daily event and supersonic and

hypersonic flights are a reality. Now aeroelastic analysis has become
an organic part of the design.
1.4 Summary of the present work

The present work involves the flutter analysis of a typical low

speed subsonic aircraft wing using a beam model as an idealized
structure representing the same. A Galerkin type of analysis using
normal mode superposition is adopted for the solution of the
differential equations of motion. Free vibration analysis of the wing is
carried out using the elementary beam model. Subsequently, the flutter
analysis of the aircraft wing is carried out using the elementary beam
model and quasi-steady aerodynamic theory, and the flutter speed is
compared with that obtained from NASTRAN. Provision is made to
account for dynamic coupling between the bending and torsional
degrees of freedom due to the fact that the shear center can be off the
centroidal axis for beams of unsymmetric sections.

As a benchmark problem, free vibration analysis of a cantilever

beam of a typical uniform rectangular section is chosen. The
NASTRAN results are generated using beam elements. Good
agreement between the results can be observed.

A typical subsonic (SARAS) aircraft wing (clean wing, with no

control surface effects) is then analysed using the program code. The
generated data for the elastic and inertial properties of the wing,
discretized as a collection of stepped beam elements, is employed for
the free vibration analysis and also the flutter analysis. These are then
compared with the results obtained from NASTRAN.
For the clean wing analysis for the SARAS aircraft, it is found that
the flutter speed is beyond the subsonic regime.i.e, it indicates the wing
does not flutter in the subsonic flow. Even for the case of a flutter
speed (determined by the present subsonic formulation and NASTRAN
code) that exceeds the limit of subsonic regime, agreement of the
results show that the computational procedure adopted here is
reasonably reliable. As a check, results are generated with reduced
stiffness parameters so that the flutter speed effectively falls in the
subsonic regime. Again results agree with those from NASTRAN,
showing the validity of the present method of analysis.
 CHAPTER -2

MATHEMATICAL FORMULATIONS

In this chapter, the mathematical formulation of subsonic flutter

analysis of a typical subsonic wing is presented. For low speed
subsonic aircrafts, the wings are usually unswept or the sweep angle
will usually be very small. A typical subsonic wing is shown in Fig 2.1
and a typical uniform rectangular wing is shown in Fig 2.2 For
aerodynamic reasons, a typical low speed subsonic wing is
characterized by high aspect ratio (span / mean chord) and a straight or
nearly straight configuration. This fact is advantageous for structural
analysis of the wing using a simple beam model, despite the complex
arrangement of the constituent structural elements. An airplane wing,
as an elastic body, has infinitely many degrees of freedom. But owing
to its particular construction, its elastic deformation in any chord wise
section can usually be described with sufficient accuracy by two
quantities: the deflection at a reference point, and the angle of
rotation about that point, i.e., the flexural and torsional deformations
respectively.

Each wing is assumed to behave like a cantilever, supported at

the axis of connectivity of the two wings, inside the fuselage. The wing
is visualized as a collection of stepped beam elements, each having its
respective elastic properties. A modal analysis method is used, using
the classical cantilever modes of beams. The present analysis is limited
to the clean wing, i.e., the ailerons are not involved in the analysis.
 ELASTIC AXIS INERTIA AXIS

Y Fig-2.1 (a)

y0 y

c
Fig-2.1 (b)

Fig-2.1 (a) A typical subsonic wing

Fig-2.1 (b) A typical aerofoil section

 U SHEAR CENTER AND CENTROID

. c
x

x
ELASTIC AXIS AND INERTIA AXIS
y

Fig-2.2 (a) Fig-2.2 (b)

w x

L
Fig-2.2(c)

Fig-2.2 (a) A typical uniform wing

Fig-2.2 (b) Section of the uniform wing

Fig-2.2 (c) Deflection of the uniform wing

 2.1 Formulation of equations of motion:

An unswept cantilever wing having a straight elastic axis

perpendicular to the fuselage, which is assumed to be fixed in space is
considered. The wing deformation can be measured by a deflection w
and a rotation about the elastic axis, w being positive downward
and is assumed positive if the leading edge up. The chordwise
displacement will be neglected. The frame of reference is chosen as
shown in Fig 2.1 (a). with the x-axis coinciding with the elastic axis.
Let y be the distance between the center of mass and the elastic axis at
any section, positive if the former lies behind the latter. Let c be the
chord length and yo be the distance of the elastic axis after the leading
edge. In a steady flow of speed U, the wing will have some elastic
deformation, which is however, of no concern to the problem of flutter.
In the following, the free motion of the wing following an initial
disturbance is considered.

Thus let w and be the deviations from the equilibrium state, and
let the inertia, elastic and aerodynamic forces correspond also to the
deviations from the steady-state values; then, for small disturbances,
the principle of superposition holds, and we have the following
equations of motion.
2 2w 2w 2
EI 2 + m 2 + my 2 + L = 0 --------- (2.1)
x 2 x t t

2 2w
GJ + I 2 + my 2 + M = 0 ---------- (2.2)
x x t t

where EI and GJ are the bending and torsional rigidity of the wing, m
and I are the mass and mass moment of inertia about the elastic axis of
the wing section at x, per unit length along the span, and L and M are
the aerodynamic lift and moment per unit span, respectively,

where,

U 2
L= cCL ----------- (2.3)
2

U 2 U 2 x
M= c CM =
2
c 2 (CM ) LE + 0 CL ----------- (2.4)
2 2 c

in which the lift and moment coefficients are given by the quasi-steady
subsonic aerodynamic theory as

dCL 1 dh 1 3 d
CL =
d + U dt + U 4 c x0 dt ----------- (2.5)

c d 1
( CM ) l . e . = CL ----------- (2.6)
8U dt 4
 The aerodynamic forces and moments are obtained by using the
quasi-steady strip theory, whereby local lift coefficient is
proportional to the instantaneous angle of attack . The derivative
dCL/d is considered to be constant, with a theoretical value of 2 for
incompressible flow and an experimental value of somewhat less than
2. Moreover the aerodynamic analysis is subject to the quasi-steady
assumption, which implies that only the instantaneous deformation is
important and the history of motion may be neglected.

Equations 2.5 and 2.6 give the lift and moment for variable w and .

From these equations we get,

2 2w 2w 2 U 2 dCL 1 w c 3 y0
EI 2 + m 2 + my 2 + c
U t + U 4 c t = 0
+
x 2 x t t 2 d
0< x<L

2w 2 U 2 2
GJ + my 2 + I 2 c
x x t t 2
c y0 1 dCL 1 w c 3 y0
+ + + = 0
8U t c 4 d U t U 4 c t
0< x< L

---------- (2.7) and (2.8)

The displacements w and are subject to the Boundary conditions

w
w= = = 0 at x=0
x

2 w 3 w
= = =0 at x= L ---------- (2.9)
x 2 x3 x

Should y and U be zero, Eqs (2.7) and (2.8) would be reduced

to two independent equations, one for w and one for . The terms
involving y and U indicate inertia and aerodynamic couplings.

Since Eqs (2.7) and (2.8) are linear equations with constant
coefficients, the solution can be written in the usual form

w( x, t ) = W ( x ) e t ( x, t ) = ( x ) e t

--------- (2.10)

where is generally complex. Introducing Eq. (2.10) into Eqs (2.7)

t
and (2.8) and dividing through out by e , we obtain the ordinary
differential equations
 U 2 dCL U dCL 3 y
( EIW ) '' ''
+
2
c
d
+
2
c
d
W + c 0 + 2 m(W + y ) = 0
4 c

---------(2.11a)

U 2 y 1 dC U
( GJ ) ' '
c2 0 L
2 c 4 d 2
y 1 dC y 1 3 y dC
c 2 0 L W + c 0 0 L +
c 4 d c 4 4 c d 8
2 ( my W + I ) = 0

-----------(2.11b)

d d2
( ) = ( ) ( ) = 2( )
' "
where and
dx dx

The boundary conditions retain the same form except that w

and are replaced by W and , respectively and partial derivatives of
w and with respect to x by total derivatives.

No closed form solution of equations (2.11) is possible; hence

an approximate solution is used. Before proceeding with the solution,
however, it will prove instructive to examine the effect of airflow
speed U on the parameter that gives the stability condition of the
system.
 2.2 Free vibration analysis of the wing

Free vibration analysis of the wing is performed using the Euler

beam model with the help of computer code as explained earlier. Good
agreement between the results confirms the validity of the present
beam model and the relevant computer code. The natural frequencies
of the equivalent beam model for one wing are determined by this code
and compared with NASTRAN.

In the beam formulation, the shear center offset from the

centroid brings about dynamic coupling between the bending and
torsional modes. This coupling is obvious from the bending and
torsional modes. No pure torsion and bending modes exist in the
structure.

For U = 0(free vibration), the above equations reduce to,

( EIW )" "

+ 2 m (W + y ) = 0
---------(2.12a)
0< x<L

( GJ ' ) + 2 ( my W + I ) = 0
'

--------- (2.12b)
0< x< L
The system can be shown to be self-adjoint and positive definite.

Using Galerkins method and assuming solution

n n+m
W = a j j = a j j ----------(2.13)
j =1 j = n +1

where j are the modal functions, which should satisfy the boundary
conditions as shown before.

The independent pure bending modes and pure torsional modes of the
uniform cantilever beam with symmetric sections (zero shear center
offset) are suitable modal functions.

Pure jth beam bending mode for the classical Euler beam with
cantilever boundary condition is given as [34]

j = (cosh an x cos an x) n (sinh an x sin an x) j=1,2n

------- (2.14a)

(cosh an L cos an L)
where, n =
(sinh an L sin an L )
Cantilever boundary conditions satisfied are

j ( x = 0) = 0 j ' ( x = 0) = 0 j = 1, 2,3.........n

The wave number an can be computed so as to satisfy the

characteristic equation,

cos( an L).cosh( an L) + 1 = 0

Pure j th torsional mode of the uniform beam is

2k 1
j = sin x j = n+k ; k = 1,2,3m
2L
----------(2.14b)

The boundary conditions satisfied here are

j ( x = 0) = 0 j = n + 1, n + 2.........n + m

'j ( x = L ) = 0

Introducing Eqs. (2.13) Into Eqs. (2.12)

n n n+m

a ( EI ) a m a my
"
j j
"
+ 2
j j + 2
j j =0 -----------(2.15a)
j =1 j =1 j = n +1

n+ m n n+m

a j ( GJ j ) + a my
'
' 2
j j + 2
a j I j = 0 -------(2.15b)
j = n +1 j =1 j = n +1
Multiplying equation (2.15a) by i (I=1,2n) and equation (2.15b) by
i (i=n+1,n+2,n+m), and integrating both results over the interval
0 x L we obtain algebraic Eigen value problem

n+m n+m

k a
j =1
ij j + 2
m a
j =1
ij j =0 (i=1,2n+m) -------(2.16)

where for the n bending modes alone

kij = i ( EI j ) dx =
L L
'' "
EI i" j " dx = kij
0 0

(i, j = 1,2n) .(2.17a)

Orthogonality of Normal modes yields the condition

kij = 0 for i j

and for the torsional modes alone,

kij = i ( GJ j ' ) dx = GJ i ' j ' dx = k ji

L ' L

0 0

i, j = n + 1, n + 2,...n + m
-------(2.17b)
Again kij = 0 for i j is the orthogonality condition
and
L
mij = m ji = mi j dx
0
i , j = 1,.n

--------(2.17c)
L
mij = m ji = my i j dx i = 1,2,n; j = n+1,n+2,...n+m
0

L
mij = m ji = I i j dx i , j = n + 1, n + 2,...n + m
0

---------(2.17d)

Orthogonality condition mij = 0 for i j

These are the symmetric stiffness and mass coefficients. Equation

(2.16) can be written in the matrix form

Ka = -2 Ma ----------(2.18)

where K and M are positive definite symmetric matrices. Hence, the

eigenvalue -2 must be real and positive, the square root of which

yields the natural circular frequency (rad/sec), {i.e., i = i2 for ith

mode}
 2.3 Flutter analysis of the wing:

Flutter analysis of the wing is also carried out using the same
elementary beam model. The quasi-steady aerodynamic theory is used
to obtain the aerodynamic forces interacting with the structure. First
the problem is solved taking one bending mode and one torsion mode
as a first estimate and then the result has been improved taking higher
modes. The same is solved in NASTRAN in pk-method. The results
obtained in NASTRAN matches the results obtained through the code.
The flutter speed obtained through the code has shown conservative
values.

The flutter speed obtained for the present configuration has

shown to be very high and also that the wing is very stiff. Hence the
stiffness of the wing is reduced by reducing the modulus of elasticity
and correspondingly the modulus of rigidity. The results obtained for
the wing with reduced stiffness parameters are typical for subsonic
flutter.

The eigenvalue is a continuous function of the air speed U.

When U is not zero, but infinitesimally small, the exponent is no
longer pure imaginary but complex, = + i. Of course, to
investigate this case, we must return to the non-self adjoint system. It
can be shown that for sufficiently small U and for (dCL / d) < 2 ,
the wing is losing energy to the surrounding air, so that the motion is
damped oscillatory, and hence asymptotically stable. The clear
implication is that is negative. As U increases, can become
positive, so that at the point at which changes sign, the motion
ceases to be damped oscillatory and becomes unstable. The air speed
corresponding to = 0 is known as critical speed and denoted by Ucr.
There are many critical values of U but, because in actual flight U
increases from an initially zero value, the lowest critical value is the
most important. One can distinguish between two critical cases,
depending on the value of . When = 0 and = 0 the wing is said
to be in critical divergent condition. When = 0 and 0 the
wing is said to be in critical flutter condition.

The above qualitative discussion can be substantiated by a

more quantitative analysis. To this end, we must derive and solve the
complete non-self adjoint eigenvalue problem. Introducing solution
(2.13) into Eqs.(2.11), multiplying Eq. (2.11a) by i (i = 1,2n) and

Eq. (2.11b) by i i (i = n+1,n+2..n+m) and integrating both results

over the interval 0 < x < L, we obtain the eigenvalue problem

[K + U2H + UL + 2M ] a = 0 ----------(2.19)
The matrices H and L are not symmetric. Their elements can be
shown to have the expressions

hij = 0 i, j = 1, 2,.....n

dCL L
hij =
2 d 0
ci j dx i = 1, 2,.....n; j = n + 1, n + 2,.......n + m

hij = 0 i = n + 1, n + 2,.......n + m; j = 1, 2,.....n

dCL L y 1
hij =
2 d 0
c 2 0 i j dx i, j = n + 1, n + 2,.......n + m
c 4

-----------.(2.20a)

dCL L
i, j = 1, 2,.....n
lij =
2 d 0
ci j dx

dCL L 3 y i = 1, 2,.....n
lij =
2 d
0
c 2 0 i j dx
4 c

j = n + 1, n + 2,.......n + m
dCL L y 1 i = n + 1, n + 2,.......n + m ;
lij =
2 d 0
c 2 0 i j dx
c 4

j = 1, 2,.....n
L y 1 3 y dC
2
lij = c 3 0 0 L i j dx
0
8 c 4 4 c d

i, j = n + 1, n + 2,.....n + m -------------(2.20b)
Using the same procedure, we can reduce the eigenvalue problem to
the standard form

k *a* = M *a* (2.21)

where
T T
a* = aT bT = aT aT (2.22)

is a 2 (n + m) vector and

0 1
K* =
( K + U H ) UL
2

1 0
M* = (2.23 a, b)
0 M

are 2 (n + m) x 2 (n + m) matrices.

The critical value Ucr of interest here is the lowest value of U for
which
= Re = 0

Flutter occurs if the real part of the eigenvalue i.e., = Re > 0

and the imaginary part of the eigenvalue i.e., = I m ( ) 0

For >0 and = 0 , divergence occurs.

The characteristic behavior of a typical mode that undergoes flutter
instability under varying airflow speeds U is shown in Fig.2.3

Stable
U<Ucr

Unstable (Flutter)
U>Ucr

Flutter boundary
U=Ucr

Fig 2.3 Behavior of typical mode amplitude when imaginary part 0

 To compute Ucr, one must solve the eigenvalue problem
repeatedly for increasing values of U. For small values of U, all the
eigen values r ( r = 1, 2,......, 2n + 2m ) have negative real parts. The first
value of u at which the real part of an eigenvalue reduces to zero is
Ucr.

0 U
Ucr

A first estimate of Ucr can be obtained by approximating W

and by means of a single term, n = m = 1 . Then, letting = i,
Ucr is the value of U, which permits the solution of the determinantal
equation

i 0 1 0
0 i 0 1
det =0
k11 U 2 h12 i m11 + Ul11 i m12 + Ul12

0 k22 + U 2 h22 i m12 + Ul21 i m22 + Ul22

------------ (2.24)
Equation (2.24) yields

4 ( m11m22 m12 2 ) i 3U m11l22 + m22l11 m12 ( l12 + l21 )

2 U 2 ( l11l22 l12l21 h12 m12 + h22 m11 ) + k22 m11 + k11m22
iU U 2 ( h12l21 + h22l11 ) + k22l11 k11l22 + k11 (U 2 h22 + k22 ) = 0

-----------(2.25)

Equating the imaginary part to zero, we obtain

U 2 ( h12l21 + h22l11 ) + k22l11 k11l22

=2

m12 ( l12 + l21 ) m11l22 + m22l11

so that substituting into the real part of Eq. 19, we can write the
quadratic equation in U2
AU 4 + BU 2 + C = 0
-----------(2.26)

where,
A = ( h12l21 + h22 l11 )
( h12l21 + h22l11 ) ( m11m22 m12 2 ) ( l11l22 l12l21 h12 m12 + h22 m11 )

m12 ( l12 + l21 ) ( m11l22 + m22l11 )
B = 2 ( h12 l21 + h22l11 )( k22 l11 k11l22 ) ( m11m22 m12 2 )
( h12l21 + h22 l11 )( k 22 m11 + k11m22 ) + ( k 22l11 k11l22 )( l11l22 l12l21 h12 m12 + h22 m11 )
m12 ( l12 + l21 ) ( m11l22 + m22l11 )

+ k11h22 m12 ( l12 + l21 ) ( m11l22 + m22l22 )

2

C = ( k22l11 k11l22 ) ( m11m22 m12 2 ) ( k22l11 k11l22 )

2

( k22 m11 + k11m22 ) m12 ( l12 + l21 ) ( m11l22 + m22l11 )

+ k11k22 m12 ( l12 + l21 ) ( m11l22 + m22l11 )
2

The solution of Eq. (2.26) is

B 1
U2 = B 2 4 AC -------------(2.27)
2A 2A

The first estimate of the critical value Ucr is the smallest positive
value of U that can be obtained.
2.4 Discretization and integration for tapered wing
(SARAS)

A typical aircraft wing (SARAS) is actually tapered along the

length. Thus, the elastic rigidity, inertia loading and aerodynamic
chord length distributions are not uniform.

To include the effects of varying section properties, the

following scheme of discretization and integration is adopted, but
using global modal functions i, satisfying cantilever boundary
conditions [equations (2.14 a, b)]

The entire wing of length L is discretized into, say, N elements.

For element r of length lr, the elastic rigidities (EI) r and (GJ) r, mass
per unit lengths mr and average chord cr are assumed to be constant
within the element.

Therefore integrals in equations (2.17) and (2.20) are

approximated by the following expressions

Bending:
L N xr +1

kij = EI i j dx ( EI )r
" "
i
"
j
"
dx
0 r =1 xr

L N xr +1

mij = mi j dx mr i, j = 1, 2, 3.........n
0 r =1
dx
xr
i j

--------------(2.28)
Torsion:

L N xr +1

kij = GJ i j dx ( GJ )r
' '
i
'
j
'
dx
0 r =1 xr

L N xr +1

mij = I i j dx ( I )r dx i j
0 r =1 xr

i , j = n + 1, n + 2........n + m
----------- (2.29)

Coupled inertia:

L N xr +1

mij = my i j dx mr y r dx
i j i = 1, 2.......n
0 r =1 xr

j = n + 1, n + 2.........n + m

---------(2.30)
The aerodynamic parameters of equations (2.20) can be likewise
obtained by piecewise integrations, where

L N xr +1

c dx c dx
0
i j
r =1
r
xr
i j

2 ( 0 )r
L
y0 1 N y 1 r+1
x

0 c c 4 i j dx i j dx
2
cr
r =1 c
r 4 x

and

L
y0 1 3 y0 dl L N
3 y0 r 1 3 y0 r r+1
x

0 c 8 c 4 4 c d i j dx cr i j dx
3

r =1 8 c r 4 4 c r xr

---------- (2.31)

These expressions are substituted in the eqn. (2.19) for flutter analysis
of tapered beam

The chord length of the wing is assumed to be varying linearly along

the length (i.e., from Root to Tip). The chord length (cr) for each
element is taken at the middle of each section (Table 3.3). These chord
lengths are substituted in the above expressions, which are in turn
substituted in Eqs (2.20).
 CHAPTER 3

NUMERICAL RESULTS

Based on the formulation given in chapter 2, a computer program is

written for the free vibration analysis and flutter analysis of aircraft
wing. The results have been validated using a standard package
NASTRAN. Further the results of some parametric studies have been
presented.

3.1 Free vibration analysis results

3.1.1 Uniform beam

In this section, to ascertain the correctness of the formulation, a bench
mark problem of an uniform cantilever beam is solved

Numerical data:
The following properties of the cantilever beam are used for the
analysis:
Length = 5m
Width = 2m
Thickness = 0.04m
Youngs Modulus of elasticity = E = 70 * 109 N/m2
Poissons ration = = 0.33
Shear Modulus of rigidity = G = 26.3 * 109 N/m2
Density of the material = s = 2700 kg/m3
Density of air = = 1.225 kg/m3
 U SHEAR CENTER AND CENTROID

x 2m

5m 0.04m

ELASTIC AXIS AND INERTIA AXIS

Fig-3.1 (a) Fig-3.1 (b)

Fig 3.1 (a) Planar view of uniform wing

Fig 3.1 (b) Sectional view of uniform wing

The natural frequencies of typical uniform beam with the properties
shown above are as shown in the following Table 3.1

Table 3.1

Natural frequencies of the uniform beam

Type and Program NASTRAN % Error b/w

Mode no. results Results Program &
in Hz in Hz NASTRAN
1 bending 1.316 1.330 1.05
2 bending 8.404 8.272 -1.59
3 bending 23.533 22.900 -2.76
4 bending 46.115 44.400 -3.86
1 torsion 6.360 6.340 -0.31
2 torsion 19.080 18.880 -1.06
3 torsion 31.798 30.960 -2.70
4 torsion 44.517 43.270 -2.81
3.1.2 The air craft wing

A Typical discretization of the aircraft wing (FE model and

Aerodynamic model) are shown in Figs 3.2 and 3.3

Numerical data:
The numerical data used for the actual wing and also for the wings
with reduced stiffness parameters are as shown below.

Case (1): ACTUAL WING

Youngs Modulus of elasticity = E = 72 * 109 N/m2

Poissons ratio = = 0.3
Shear Modulus of rigidity = G = 27.69 * 109 N/m2

Case (2a): WITH REDUCED STIFFNESS PARAMETERS

E* = 0.5E and G* = 0.5G

Youngs Modulus of elasticity = E* = 36 * 109 N/m2
Poissons ratio = = 0.3
Shear Modulus of rigidity = G* = 13.85 * 109 N/m2

Case (2b): WITH REDUCED STIFFNESS PARAMETERS

E* = 0.1E and G* = 0.1G

Youngs Modulus of elasticity = E* = 7.2 * 109 N/m2
Poissons ratio = = 0.3
Shear Modulus of rigidity = G* = 2.769 * 109 N/m2
Fig 3.2 An FE model of complete SARAS wing

Fig 3.3 Aerodynamic model of SARAS wing

 Table 3.2
Mass distribution and mass densities of the beam element
Ref [32] PD ST 0314

Sl. No. Ele L in C/S Area Density in Mass per

m (lr) in m2 Kg/m3 unit length
x 10-6(Ar) (s) r Kg/m (mr)
1 0.350 12871.0 9309.30 119.82
2 0.315 12037.0 7803.98 93.94
3 0.285 12318.0 9609.27 118.37
4 0.300 11400.0 9784.31 111.54
5 0.325 15438.0 7654.20 118.16
6 0.315 7348.1 7968.55 58.55
7 0.315 6826.9 49346.43 336.88
8 0.325 6470.2 41645.03 269.45
9 0.325 5314.3 44344.80 235.66
10 0.325 5110.9 40029.30 204.58
11 0.325 4947.1 34392.09 170.14
12 0.325 4817.5 29235.08 14.84
13 0.325 4153.8 26966.23 112.01
14 0.325 3939.4 22792.80 89.79
15 0.325 3455.0 17246.97 59.58
16 0.350 3259.6 9963.28 32.47
17 0.350 3132.6 6385.20 20.00
18 0.300 3153.5 6178.57 19.48
19 0.210 2803.0 9275.77 25.99
20 0.350 2557.7 6660.95 17.04
21 0.370 2009.6 7063.66 14.20
22 0.370 2119.7 6485.07 13.75

Total wing mass = 762.6 kg

 Table 3.3
Sectional properties and aerodynamic chord lengths
of elements of the wing

Ref [32] PD ST - 0314

Sl. Ele L in Izz in m4 Iyy in m4 J in m4 Chord

No. m x 10-4 x 10-4 x 10 -4
length(m)
1 0.350 4.0477 35.3980 70.3185 2.402
2 0.315 3.2949 29.8140 45.8763 2.326
3 0.285 3.1249 29.5620 25.7400 2.256
4 0.300 2.8388 23.3110 15.0045 2.191
5 0.325 3.7622 23.3330 9.0301 2.120
6 0.315 1.5680 11.6677 7.1700 2.047
7 0.315 1.4159 10.4640 5.6790 1.975
8 0.325 1.2581 9.3388 5.2815 1.902
9 0.325 0.9428 6.8012 4.3770 1.828
10 0.325 0.8215 6.4647 3.8130 1.754
11 0.325 0.7414 6.0293 3.3300 1.680
12 0.325 0.6706 5.4199 2.7525 1.606
13 0.325 0.5127 4.1906 2.3640 1.532
14 0.325 0.4474 3.6573 2.0250 1.458
15 0.325 0.3453 2.8600 1.6335 1.384
16 0.350 0.2958 2.3924 1.3110 1.307
17 0.350 0.2537 2.0252 1.0515 1.227
18 0.300 0.2295 1.8596 0.8730 1.153
19 0.210 0.1901 1.5036 0.7470 1.095
20 0.350 0.1525 1.2544 0.5865 1.031
21 0.370 0.1036 0.7833 0.4125 0.949
22 0.370 0.9409 0.6224 0.2940 0.865
 Table 3.4
Shear center position w.r.t. Centroidal axis

Ref [32] PD ST - 0314

Sl. Ele L in ZG in m YG in m
No. m x 10-3 x 10-3
1 0.350 -355.430 36.300
2 0.315 -432.248 -137.681
3 0.285 -379.560 -199.990
4 0.300 -434.200 -278.100
5 0.325 -161.000 -9.200
6 0.315 -1.648 39.906
7 0.315 -9.040 11.800
8 0.325 -6.333 7.034
9 0.325 -8.275 -3.881
10 0.325 -19.418 -0.664
11 0.325 -18.830 -3.497
12 0.325 -10.425 13.077
13 0.325 -16.588 -4.718
14 0.325 -14.829 -5.609
15 0.325 -6.369 -12.955
16 0.350 -0.028 -10.495
17 0.350 -2.468 8.142
18 0.300 2.517 -8.635
19 0.210 1.100 -7.300
20 0.350 -1.173 -6.608
21 0.370 -0.851 -5.770
22 0.370 -1.952 -5.314
 The above numerical data are used for the analysis of the
subsonic wing. The wing is visualized as a collection of stepped beam
elements, each having its respective properties as shown in the above
tables. The natural frequencies obtained for the wing for each case are
given below.

Table 3.5
Natural frequencies of the subsonic wing case (1)

Type and Present Stick model 3D model

Mode no. analysis results in Hz results in Hz
results Ref [32] Ref[32]
in Hz Table 4.3 Table 4.3
PD ST-0314 PD ST-0314
1 bending 7.331 7.119 7.087
2 bending 21.227 20.786 20.481
3 bending 49.800 48.538 47.781
4 bending 127.767 - -
1 torsion 57.168 56.500 56.338
2 torsion 123.302 - -
3 torsion 183.103 - -
4 torsion 388.522 - -

Type and %Error b/w %Error b/w

Mode no. program program and
and stick 3D model
model
1 bending 2.89 3.33
2 bending 2.07 3.51
3 bending 2.53 4.05
4 bending - -
1 torsion 1.17 1.45
2 torsion - -
3 torsion - -
4 torsion - -
 Table 3.6
Natural frequencies of the subsonic wing with E*=0.5E case (2a)

Type and Program NASTRAN % Error b/w

Mode no. results results Program &
in Hz in Hz NASTRAN
1 bending 5.180 5.03 -2.98
2 bending 15.009 14.60 -2.80
3 bending 35.930 34.10 -5.36
4 bending 103.49 102.20 -1.26
1 torsion 40.420 39.90 -1.30
2 torsion 90.132 92.35 2.40
3 torsion 128.800 129.67 0.67
4 torsion 274.580 266.00 -3.22

Table 3.7
Natural frequencies of the subsonic wing with E*=0.1E case (2b)

Type and Program NASTRAN % Error b/w

Mode no. results results Program &
in Hz in Hz NASTRAN
1 bending 2.31 2.24 -3.12
2 bending 6.71 6.50 -3.23
3 bending 16.07 15.20 -5.72
4 bending 44.56 43.32 -2.86
1 torsion 18.07 17.89 -1.00
2 torsion 39.13 38.80 -0.85
3 torsion 57.62 57.99 0.63
4 torsion 122.79 119.15 -3.05
3.2 Flutter analysis

3.2.1 Uniform beam:

The flutter analysis of the uniform rectangular beam is carried

out using the modal analysis method. The shear center in the case of
uniform section coincides with the centroid of the section. Hence there
will be no dynamic coupling in the case of uniform beam.
The Eigen values are obtained for increasing values of
velocities. The values obtained are plotted. The velocity v/s the real
part of the eigenvalue which is indicative of damping and the velocity
v/s the imaginary part of the eigen value which is the frequency in
rad/s are plotted from the complex eigen values obtained from the
present analysis for different modes. The same graphs are also plotted
for the values obtained from NASTRAN for different modes. The
graphs plotted are shown in the Figs 3.4 (a) and 3.4 (b)
The flutter speeds obtained from the present analysis and
through NASTRAN are shown in the following Table 3.8

Table 3.8
Uniform beam Flutter results

Present analysis NASTRAN

Flutter speed (m/s) 148 170
NASTRAN

UNIFORM BEAM

1.00E+00
REAL PART OF THE EIGEN VALUES ( )---->

0.00E+00
0.00E+00 5.00E+01 1.00E+02 1.50E+02 2.00E+02 2.50E+02

-1.00E+00
2nd bending
-2.00E+00 2nd torsion
1st torsion
3rd torsion
-3.00E+00
1st bending
3rd bending
-4.00E+00
4th bending

-5.00E+00

-6.00E+00

-7.00E+00
SPEED (U) (m/s)

PRESENT ANALYSIS

Fig 3.4 (a) Plot of Velocity (U) v/s real part () of the eigen value
NASTRAN

UNIFORM BEAM

IMAG PART OF THE EIGEN VALUE ()------> 300

250

1st bending
200 1st torsion

3rd bending

4th torsion
150
2nd bending

2nd torsion
100
3rd torsion

4th bending

50

0
0 20 40 60 80 100 120 140 160 180 200
SPEED (U) (m/s)

PRESENT ANALYSIS

Fig 3.4 (b) Plot of Velocity (U) v/s imaginary part () of the eigen value
3.2.2 Aircraft wing

The flutter analysis of the wing is also carried in the same way
with the inclusion of effect of shear center offset. The clean wing
modeled as having stepped beam elements is analysed as mentioned
before. The velocity v/s the real part of eigen values and velocity v/s
the imaginary part of eigen values are plotted from the complex eigen
values obtained from the present analysis and also from the
NASTRAN.

The velocity v/s damping curves i.e., v-g curves and the velocity
v/s frequency curves i.e., v-f curves can also be plotted. The relation
between the damping (g) values and eigen values and the relation
between the frequency (f) values and eigen values are as given below.

If the eigen value obtained is = +i, where is the real part

and is the imaginary part, then

g= x 2 and f =
2

Hence the v-g and v-f curves can be plotted from the eigen
values and the same can also be obtained from NASTRAN.
The velocity v/s Real part and velocity v/s Imaginary part are
also plotted for the wing with reduced stiffness parameters. All the
graphs obtained are shown below in the Figs 3.5, 3.6, and 3.7
NASTRAN

SARAS WING

REAL PART OF THE EIGEN VALUES ( )--- 10

0
0 100 200 300 400 500 600 700 800
2nd bending
-5

3rd bending
-10
1st torsion
>

-15
3rd torsion

-20 1st bending

2nd torsion
-25

-30

-35

SPEED (U) (m/s)

PRESENT ANALYSIS

Fig 3.5 (a) Plot of Velocity (U) v/s real part () of the eigen value
NASTRAN

SARAS WING

1200
IMAG PART OF THE EIGEN VALUE ( )------>

1000

800
1st bending
2nd bending
3rd bending
600 4th bending
1st torsion
2nd torsion
400 3rd torsion

200

0
0 100 200 300 400 500 600 700

SPEED (U) (m/s)

PRESENT ANALYSIS

Fig 3.5 (b) Plot of Velocity (U) v/s imaginary part () of the eigen value
NASTRAN

WING WITH REDUCED STIFFNESS PARAMETERS(0.5 E)

10
REAL PART OF THE EIGEN VALUES ( )---->

0
0 100 200 300 400 500 600 2nd bending
3rd bending
1st torsion
-5
1st bending
3rd torsion
2nd torsion
-10

-15

-20

SPEED (U) (m/s)

PRESENT ANALYSIS

Fig 3.6 (a) Plot of Velocity (U) v/s real part () of the eigen value
NASTRAN

WING WITH REDUCED STIFFNESS PARAMETERS(0.5E and 0.5G)

900
IMAG PART OF THE EIGEN VALUE ( )------>

800

700

1st bending
600
2nd bending
500 3rd bending
2nd torsion
400 1st torsion
4th bending
300
3rd torsion

200

100

0
0 100 200 300 400 500 600
SPEED (U) (m/s)

PRESENT ANALYSIS

Fig 3.6 (b) Plot of Velocity (U) v/s imaginary part () of the eigen value
NASTRAN

WING WITH REDUCED STIFFNESS PARAMETERS (0.1E)

2
REAL PART OF THE EIGEN VALUES ( )---->

0
0 50 100 150 200 250

-2 2nd bending
3rd bending
1st torsion
-4 2nd torsion
3rd torsion
1st bending
-6

-8

-10
SPEED (U) (m/s)

PRESENT ANALYSIS

Fig 3.7 (a) Plot of Velocity (U) v/s real part () of the eigen value
NASTRAN

WING WITH REDUCED STIFFNESS PARAMETERS (0.1E and 0.1G)

400
IMAG PART OF THE EIGEN VALUE ()------>

350

300

1st bending
250 2nd bending
3rd bending

200 1st torsion

2nd torsion
3rd torsion
150
4th bending

100

50

0
0 50 100 150 200 250

SPEED (U) (m/s)

PRESENT ANALYSIS

Fig 3.7 (b) Plot of Velocity (U) v/s imaginary part () of the eigen value
The flutter speeds obtained from the graphs for the actual wing and for
the wing with reduced stiffness parameters are shown in the following
Table 3.9

Table 3.9

Flutter results of the subsonic wing

Flutter speeds (m/s) Present analysis NASTRAN

Aircraft wing 650 680

Wing with reduced stiffness 450 470

parameters (0.5E and 0.5G)
Wing with reduced stiffness 208 223
parameters (0.1E and 0.1G)
 CHAPTER 4

CONCLUSIONS

4.1 DISCUSSION OF RESULTS

In the present analysis of subsonic aircraft wing, classical Euler

beam model is used. The wing is discretized as a collection of stepped
beam elements and is employed for free vibration analysis and the
flutter analysis.

The beam model used yields reasonably good results for natural
frequencies. The natural frequencies obtained are in good agreement
with previously obtained natural frequencies of the clean wing through
NASTRAN (Ref [32] PD ST-0314) showing the validity of the present
analysis.

For the flutter analysis of the clean wing, the aerodynamic forces
are obtained by the quasi-steady subsonic aerodynamic theory. The
flutter results obtained by the present method are compared with those
obtained from the NASTRAN. The present analysis gives a
conservative value of flutter speed as compared to NASTRAN values.
The flutter graphs obtained from present analysis and through
NASTRAN are shown in the previous chapter. The comparison shows
the validity of present method for flutter analysis.
 For the clean wing analysis of the SARAS aircraft, the flutter
speed is found to be beyond the subsonic regime, i.e., the wing does
not flutter in the subsonic flow. Even for the case of flutter speed that
exceeds the limit of subsonic regime, agreement of results show that
the method adopted is reasonably reliable. Analysis is carried out with
reduced stiffness parameters so that the flutter speed falls in the
subsonic regime. Again results agree with those from NASTRAN.

4.2 FURTHER SCOPE

1 The present work is limited to the clean wing analysis that

doesnt show flutter in the subsonic regime. However it is
necessary to check if the wing with control surfaces is prone
to subsonic flutter. The present method can be easily
extended to determine flutter boundaries of wing with control
surfaces

2 The flutter analysis of the T-tail is critical from the point of

design. The quasi-steady method can easily be extended to
the T-Tail assembly consisting of Horizontal tail, Vertical
tail, Rudder and Elevator.
REFERENCES:

1. Lanchester, F.W.: Torsional vibration of the Tail of an

Aeroplane. Aeronaut.Research Com.R & M.276, part i (July
1916)

2. Bairstow, L., and A.Fage: Oscillations of the Tail Plane and

Body of an Aeroplane in Flight. Aeronaut.Research Com.R
& M.276, part ii (July 1916)

3. Blasius,H: Umber Schwingungsercheiningen an Einholmigen

Unterflugeln.Z.Flugtech.u.Motorluftschif.16,39-42 (1925)

4. Glauret,H.: The Accelerated Motion of a Cylindrical Body

through a Fluid. Aeronaut.Research Com.R & M.1215 (1929)

5. Glauret,H.: The Force and Moment of an Oscillating

Aerofoil. Aeronaut.Research Com.R & M.1242 (1929)

6. Theodorsen, Th.: General Theory of Aerodynamic Instability

and the Mechanism of Flutter. NACA Rept.496 (1934)

7. Smilg, B.: The Instability of Pitching Oscillations of an

Airfoil in Subsonic Incompressible Potential Flow.
J.Aeronaut.Sci.16,691-696 (Nov.1949)
8. Cheilik,H., and H.Frissel: Theoretical Criteria for Single
Degree of Freedom Flutter at Supersonic Speeds. Cornell
Aeronaut.Lab.Rept.CAL-7A (May 1947)

9. Runyan,H.L.,H.J.Cunningham, and C.E.Watkins: Theoretical

Investigation of Several Types of Single-Degree-of-Freedom
Flutter.J.Aeronaut.Sci.19,101-110, 126 (1952) Comments by
K.P.Abichandani and R.M.Rosenberg and authors reply,
215-216; 503-504

10. Cunningham, H.J.: Analysis of Pure-Bending Flutter of a

Cantilever Swept Wing and Its Relation to Bending-Torsion
Flutter. NACA Tech.Note 2461 (1951)

11. Greidanus, J.H.: Low-Speed Flutter. J.Aeronaut.Sci.16,127-

128 (1949)

12. Y.C.Fung, An Introduction to the theory of aero elasticity,

John Wiley and sons, Inc

13. E.H.Dowell, Non linear oscillation of a fluttering plate,

AIAA J. Vol.4,No.7,July 1996.

14. E.H.Dowell,Theoretical and experimental panel flutter study,

AIAA J. Vol.3,No.12,Dec 1995.

15. L.Mirovitch, Elements of Vibration analysis, McGraw Hill

Book Co.,1975.
16. D.J.Evins, Modal Testing: Theory and Practice, Research
and Study Press Ltd. John Wiley and Sons Inc,1986.

17. R.L.Bisplinghoff and H.Ashley, Aero elasticity; Addison-

Wesley Publication 1957.

18. K.F.Alvin and L.D.Paterson, Method for determining

minimum order mass and stiffness Matrices from modal test
data: AIAA J. Vol.33,No.1,Jan 1995.

19. John Dugundji, Theoretical consideration of Panel flutter at

high supersonic Mach No., AIAA J. Vol.4,No.7,July 1966

20. J.B.Kosmatka, Flexure-Torsion behavior of prismatic beam,

AIAA J. Vol.31,No.1,Jan 1993.

21. Abott.I.H. and Von Doenhoff, Theory of Wing section,

McGraw Hill,Newyork,1949.

22. S.H.R.Eslimy-Islahany and A.J.Sobby, Response of bending

torsion coupled beam to deterministic and random loads, J. of
Sound and Vibration,195(2),1996.

23. E.Dokumaci, An Exact solution for coupled bending and

torsion vibration of uniform beam having single cross section
symmetry, J.of Sound and Vibration,119(3),1987.
24. Shyamlal Chatterji, Active stabilization technique of
dynamic machine, M.Tech.Thesis, IIT Kanpur, March 1991.

25. W.Wilson, Calculation and Design of Electrical Apparatus,

Chapman and Hall Ltd.,1955.

26. Haisler W.E. and Allen D.H., Introduction to Aerospace

structural analysis, John Wiley Publication., 1985.

27. Paul Kuhn, Stresses in Aircraft and Shell structure, McGraw

Hill Book Co.,Inc.

28. W.T.Thomson, Theory of Vibration with application,

George Allen and Unwin, Prantice Hall Publication, 1983.

29. Hurty W.C. and Rubinstein M.F., Dynamics of Structures,

Prantice Hall of India Pvt.Ltd.,1967.

30. S.Talukedar, S.Kamle and D.Yadav (1996), Track induced

heave pitch dynamics of vehicles with variable section
flexible attachment, Vehicle System Dynamics, (Paper to be
published)

31. Dowell E.H., D.A.Peters, R.H.Scanlan and F.Sisto, A

Modern course in Aeroelasticity, III edition, Kluwer
Academic Publishers.
32. Dr. S.Mukherjee, Manju, NAL Project Document, Dynamic
characterization SARAS wing and empennage. PD ST-0314

33. L.Mirovitch, Theory of Aero elasticity, McGraw Hill Book

Co., 1975.

34. Mario Paz, Structural Dynamics, Theory and Computation,II

edition,CBS Publishers and Distributors.

35. Ray W.Clough and Joseph Penzien, Dynamics of Structures,

II edition, McGraw-Hill International editions.

36 . Lecture notes by Dr.S.Mukherjee.

37. Dewey H.Hodges and G.Alvin Pierce, Introduction to

Structural Dynamics and Aeroelasticity, Cambridge University
Press.
 APPENDIX I A

FORTRAN code for generation of stiffness, mass and aerodynamic matrices

! AEROELASTIC FLUTTER ANALYSIS OF TAPERED WING(SARAS)

! QUASI STEADY AERODYNAMIC THEORY IS USED FOR AERODYNAMIC MATRICES
! NUMERICAL INTEGRATION (TRAPEZOIDAL RULE)
! IS USED FOR INTEGRATION OF THE MODAL FUNCTIONS

DIMENSION FUN(51000),SUMKT(10,10),SUMMT(10,10)
DIMENSION SUMKB(10,10),SUMMB(10,10),SUMMC(10,10)
DIMENSION STIFF(50,50),FMASS(50,50)
DIMENSION SUMH2(10,10),SUMH4(10,10)
DIMENSION SUML1(10,10),SUML2(10,10),SUML3(10,100),SUML4(10,10)
DIMENSION FHMAT(10,10),FLMAT(10,10)
DIMENSION TORSIONK(10,10),TORSIONM(10,10)
DIMENSION BENDINGK(10,10),BENDINGM(10,10)
DIMENSION COUPLEDM(10,10)
DIMENSION AEROH2(10,10),AEROH4(10,10)
DIMENSION AEROL1(10,10),AEROL2(10,10),AEROL3(10,10),AEROL4(10,10)
DIMENSION STIFF2(50,50),FMASS2(50,50)
DIMENSION WR(100),WI(100),FMASS1(50,50),STIFF1(50,50)
DIMENSION anl(10),a(10)

OPEN(1,FILE='wing.in')
OPEN(2,FILE='wing.out')

WRITE (*,*) 'ENTER THE TOTAL LENGTH OF THE BEAM'

READ (*,*) UL
WRITE (*,*) 'ENTER THE NUMBER OF DIVISIONS OF THE BEAM'
READ (*,*) NNN

DO 9999 N = 1,4
DO 9999 M = 1,4
TORSIONK (N, M) = 0.0
TORSIONM (N, M) = 0.0
BENDINGK (N, M) = 0.0
BENDINGM (N, M) = 0.0
COUPLEDM (N, M) = 0.0
AEROH2 (N, M) = 0.0
AEROH4 (N, M) = 0.0
AEROL1 (N, M) = 0.0
AEROL2 (N, M) = 0.0
AEROL3 (N, M) = 0.0
AEROL4 (N, M) = 0.0
9999 CONTINUE

DO 5555 III = 1,NNN

! Reading all the data of each section (material and inertia properties of each section)

WRITE (*,*) 'ENTER THE LOWER AND UPPER LIMITS OF SECTION',III

READ (1,*) FL,U
 WRITE (*,*) 'ENTER THE NUMBER OF DIVISIONS'
READ (1,*) ND
WRITE (*,*) 'ENTER THE NUMBER OF MODES'
READ (1,*) MO
WRITE (*,*) 'ENTER THE MODULUS OF ELASTICITY'
READ (1,*) E
WRITE (*,*) 'ENTER THE MODULUS OF RIGIDITY'
READ (1,*) G
WRITE (*,*) 'ENTER THE DENSITY OF THE AIR '
READ (1,*) RHOA
WRITE (*,*) 'ENTER THE ASPECT RATIO'
READ (1,*) AR
WRITE (*,*) 'ENTER IZZ'
READ (1,*) FIZZ
WRITE (*,*) 'ENTER THE CROSS SECTION AREA OF EACH DIVISION'
READ (1,*) AREA
WRITE (*,*) 'ENTER THE DENSITY OF THE MATERIAL'
READ (1,*) RHOM
WRITE (*,*) 'ENTER J'
READ (1,*) FJ
WRITE (*,*) 'ENTER IXX'
READ (1,*) FIXX
WRITE (*,*) 'ENTER THE SHEAR CENTER OFFSET'
READ (1,*) YTHETA
WRITE (*,*) 'ENTER THE CHORD LENGTH OF THE ELEMNET'
READ (1,*) CHORD

PI=4.0*ATAN (1.0)
!----------------------------------------------------------------------------------
DO 500 N=1,MO
DO 500 M=1,MO
H=(U-FL)/ND
N1=ND+1
S1=0.0
S2=0.0
S3=0.0
S4=0.0
S5=0.0
S6=0.0
S7=0.0
S8=0.0
S9=0.0
S10=0.0
S11=0.0

PI=4.0*ATAN(1.0)

! stiffness matrix for torsion -SUMKT(N,M)

! TORSIONK(N,M) - summation of integrals SUMKT(N,M) of each section
DO 50 I=1,N1
X=FL+(I-1)*H
CALL KT(N,M,X,UL,G,FJ,FKT)
FUN(I)=FKT
S1=2.0*FUN(I)+S1
50 CONTINUE
S1=S1-FUN(1)-FUN(N1)
SUMKT(N,M) = S1*H/2.0
TORSIONK(N,M) = TORSIONK(N,M) + SUMKT(N,M)
! mass matrix for torsion -SUMMT(N,M)
! TORSIONM(N,M) - summation of integrals SUMMT(N,M) of each section

DO 60 I=1,N1
X=FL+(I-1)*H
CALL MT(N,M,X,UL,RHOM,FIXX,FFMT)
FUN(I)=FFMT
S2=2.0*FUN(I)+S2
60 CONTINUE
S2=S2-FUN(1)-FUN(N1)
SUMMT(N,M) = S2*H/2.0
TORSIONM(N,M) = TORSIONM(N,M) + SUMMT(N,M)

! stiffness matrix for bending

DO 70 I=1,N1
X=FL+(I-1)*H
CALL KB(N,M,X,UL,FIZZ,E,FKB)
FUN(I)=FKB
S3=2.0*FUN(I)+S3
70 CONTINUE
S3=S3-FUN(1)-FUN(N1)
SUMKB(N,M) = S3*H/2.0
BENDINGK(N,M) = BENDINGK(N,M) + SUMKB(N,M)

! mass matrix for bending

DO 80 I=1,N1
X=FL+(I-1)*H
CALL MB(N,M,X,UL,RHOM,AREA,FMB)
FUN(I)=FMB
S4=2.0*FUN(I)+S4
80 CONTINUE
S4=S4-FUN(1)-FUN(N1)
SUMMB(N,M) = S4*H/2.0
BENDINGM(N,M) = BENDINGM(N,M) + SUMMB(N,M)
! coupled mass matrix
DO 90 I=1,N1
X=FL+(I-1)*H
CALL MC(N,M,X,UL,AREA,YTHETA,FMC)
FUN(I)=FMC
S5=2.0*FUN(I)+S5
90 CONTINUE
S5=S5-FUN(1)-FUN(N1)
SUMMC(N,M) = S5*H/2.0
COUPLEDM(N,M) = COUPLEDM(N,M) + SUMMC(N,M)

! aerodynamic matrices

DO 200 I=1,N1
X=FL+(I-1)*H
CALL HH2(N,M,X,UL,FH2)
FUN(I)=FH2
S6=2.0*FUN(I)+S6
200 CONTINUE
S6=S6-FUN(1)-FUN(N1)
SUMH2(N,M) = ((RHOA*PI)*S6*H/2.0)
AEROH2(N,M) = AEROH2(N,M) + SUMH2(N,M)
 DO 210 I=1,N1
X=FL+(I-1)*H
CALL HH4(N,M,X,UL,YTHETA,FH4)
FUN(I)=FH4
S7=2.0*FUN(I)+S7
210 CONTINUE
S7=S7-FUN(1)-FUN(N1)
SUMH4(N,M) = ((RHOA*PI)*S7*H/2.0)
AEROH4(N,M) = AEROH4(N,M) + SUMH4(N,M)

DO 220 I=1,N1
X=FL+(I-1)*H
CALL LL1(N,M,X,UL,FFL1)
FUN(I)=FFL1
S8=2.0*FUN(I)+S8
220 CONTINUE
S8=S8-FUN(1)-FUN(N1)
SUML1(N,M) = RHOA*PI * S8*H/2.0
AEROL1(N,M) = AEROL1(N,M) + SUML1(N,M)

DO 230 I=1,N1
X=FL+(I-1)*H
CALL LL2(N,M,X,UL,YTHETA,FFL2)
FUN(I)=FFL2
S9=2.0*FUN(I)+S9
230 CONTINUE
S9=S9-FUN(1)-FUN(N1)
SUML2(N,M) = RHOA * PI * S9*H/2.0
AEROL2(N,M) = AEROL2(N,M) + SUML2(N,M)

DO 240 I=1,N1
X=FL+(I-1)*H
CALL LL3(N,M,X,UL,YTHETA,FFL3)
FUN(I)=FFL3
S10=2.0*FUN(I)+S10
240 CONTINUE
S10=S10-FUN(1)-FUN(N1)
SUML3(N,M) = -( (RHOA * PI ) * S10*H/2.0)
AEROL3(N,M) = AEROL3(N,M) + SUML3(N,M)

DO 250 I=1,N1
X=FL+(I-1)*H
CALL LL4(N,M,X,UL,DCL,YTHETA,FFL4)
FUN(I)=FFL4
S11=2.0*FUN(I)+S11
250 CONTINUE
S11=S11-FUN(1)-FUN(N1)
SUML4(N,M) = (0.5*RHOA)*S11*H/2.0
AEROL4(N,M) = AEROL4(N,M) + SUML4(N,M)

500 CONTINUE
 WRITE (2,*) '"STIFFNESS MATRIX FOR PURE TORSION"'
DO 111 N=1,MO
DO 110 M=1,MO
WRITE (2,15) TORSIONK(N,M)
110 CONTINUE
WRITE (2,*)
111 CONTINUE

WRITE (2,*) '"MASS MATRIX FOR PURE TORSION"'

DO 121 N=1,MO
DO 120 M=1,MO
WRITE (2,15) TORSIONM(N,M)
120 CONTINUE
WRITE (2,*)
121 CONTINUE

WRITE (2,*) '"STIFFNESS MATRIX FOR PURE BENDING"'

DO 131 N=1,MO
DO 130 M=1,MO
WRITE (2,15) BENDINGK(N,M)
130 CONTINUE
WRITE (2,*)
131 CONTINUE

WRITE (2,*) '"MASS MATRIX FOR PURE BENDING"'

DO 141 N=1,MO
DO 140 M=1,MO
WRITE (2,15) BENDINGM(N,M)
140 CONTINUE
WRITE (2,*)
141 CONTINUE

WRITE (2,*) '"COUPLED MASS MATRIX"'

DO 151 N=1,MO
DO 150 M=1,MO
WRITE (2,15) COUPLEDM(N,M)
150 CONTINUE
WRITE (2,*)
151 CONTINUE

WRITE (2,*) '"H2 MATRIX"'

DO 301 N=1,MO
DO 300 M=1,MO
WRITE (2,16) AEROH2(N,M)
300 CONTINUE
WRITE (2,*)
301 CONTINUE

WRITE (2,*) '"H4 MATRIX"'

DO 311 N=1,MO
DO 310 M=1,MO
WRITE (2,16) AEROH4(N,M)
310 CONTINUE
311 CONTINUE
 WRITE (2,*) '"L1 MATRIX"'
DO 321 N=1,MO
DO 320 M=1,MO
WRITE (2,16) AEROL1(N,M)
320 CONTINUE
WRITE (2,*)
321 CONTINUE

WRITE (2,*) '"L2 MATRIX"'

DO 331 N=1,MO
DO 330 M=1,MO
WRITE (2,16) AEROL2(N,M)
330 CONTINUE
WRITE (2,*)
331 CONTINUE

WRITE (2,*) '"L3 MATRIX"'

DO 341 N=1,MO
DO 340 M=1,MO
WRITE (2,16) AEROL3(N,M)
340 CONTINUE
WRITE (2,*)
341 CONTINUE

WRITE (2,*) '"L4 MATRIX"'

DO 351 N=1,MO
DO 350 M=1,MO
WRITE (2,16) AEROL4(N,M)
350 CONTINUE
WRITE (2,*)
351 CONTINUE

15 FORMAT (F15.4,$)
16 FORMAT (F19.15,$)

5555 CONTINUE

STOP
END
!--------------------------------------------------------------------------------------------------------------
! SUBROUTINES FOR THE COMPUTATION OF MODAL FUNCTIONS
!-------------------------------------------------------------------------------------------------------------

! STIFFNESS MATRIX FOR BENDING

SUBROUTINE KB(N,M,X,U,FIZZ,E,FKB)

DIMENSION anl(10),a(10)

PI=4.0*ATAN(1.0)

anl(1) = 1.8753
anl(2) = 4.6941
anl(3) = 7.8547
anl(4) = 10.9953

SIG1 = ( COS(anl(N)) + COSH(anl(N)) ) / ( SIN(anl(N)) +SINH(anl(N)) )

SIG2 = ( COS(anl(M)) + COSH(anl(M)) ) / ( SIN(anl(M)) +SINH(anl(M)) )

a(1) = 1.8753/U
a(2) = 4.6941/U
a(3) = 7.8547/U
a(4) = 10.9953/U

FK1 = (COSH(a(N)*X)+COS(a(N)*X)) - (SIG1 * (SINH(a(N)*X) + SIN(a(N)*X)))

FK2 = (COSH(a(M)*X)+COS(a(M)*X)) - (SIG2 * (SINH(a(M)*X) + SIN(a(M)*X)))

AA = (a(N))**2
BB = (a(M))**2

FKB = E * FIZZ * (FK1*FK2) * (AA) * (BB)

RETURN
END
!-----------------------------------------------------------------------
! MASS MATRIX FOR BENDING
SUBROUTINE MB(N,M,X,U,RHOM,AREA,FMB)

DIMENSION anl(10),a(10)

PI=4.0*ATAN(1.0)

anl(1) = 1.8753
anl(2) = 4.6941
anl(3) = 7.8547
anl(4) = 10.9953

SIG1 = ( COS(anl(N)) + COSH(anl(N)) ) / ( SIN(anl(N)) +SINH(anl(N)) )

SIG2 = ( COS(anl(M)) + COSH(anl(M)) ) / ( SIN(anl(M)) +SINH(anl(M)) )

a(1) = 1.8753/U
a(2) = 4.6941/U
a(3) = 7.8547/U
a(4) = 10.9953/U

FM1 = (COSH(a(N)*X)-COS(a(N)*X)) - (SIG1 * (SINH(a(N)*X) - SIN(a(N)*X)))

FM2 = (COSH(a(M)*X)-COS(a(M)*X)) - (SIG2 * (SINH(a(M)*X) - SIN(a(M)*X)))
 FMB = (RHOM*AREA) * FM1 * FM2

RETURN
END
!------------------------------------------------------------------------
! STIFFNESS MATRIX FOR TORSION
SUBROUTINE KT(N,M,X,U,G,FJ,FKT)

PI=4.0*ATAN(1.0)

A = ((2.0*N-1)*PI*X)/(2.0*U)
B = ((2.0*M-1)*PI*X)/(2.0*U)

C = (2.0*N-1)*PI/(2.0*U)
D = (2.0*M-1)*PI/(2.0*U)

FKT = G * FJ * (COS(A)*COS(B)) * C * D

RETURN
END
!------------------------------------------------------------------------
! MASS MATRIX FOR TORSION

SUBROUTINE MT(N,M,X,U,RHOM,FIXX,FFMT)

PI=4.0*ATAN(1.0)

A = ((2*N-1)*PI*X)/(2*U)
B = ((2*M-1)*PI*X)/(2*U)

FFMT = (RHOM*FIXX) * (SIN(A) * SIN(B))

RETURN
END
!--------------------------------------------------------------------
! COUPLED MASS MATRIX (IF THERE IS SHEAR CENTER OFFSET)
SUBROUTINE MC(N,M,X,U,AREA,YTHETA,FMC)
DIMENSION anl(10),a(10)

PI=4.0*ATAN(1.0)

anl(1) = 1.8753
anl(2) = 4.6941
anl(3) = 7.8547
anl(4) = 10.9953

SIG1 = ( COS(anl(N)) + COSH(anl(N)) ) / ( SIN(anl(N)) +SINH(anl(N)) )

a(1) = 1.8753/U
a(2) = 4.6941/U
a(3) = 7.8547/U
a(4) = 10.9953/U

B = ((2*M-1)*PI*X)/(2*U)

FM1 = (COSH(a(N)*X)-COS(a(N)*X)) - (SIG1 * (SINH(a(N)*X) - SIN(a(N)*X)))

FM2 = SIN(B)

FMC = (RHOM*AREA) * YTHETA * FM1*FM2

 RETURN
END

!-------------------------------------------------------------------------
! AERO DYNAMIC MATRICES

SUBROUTINE HH2(N,M,X,U,CHORD,FH2)
DIMENSION anl(10),a(10)

PI=4.0*ATAN(1.0)

anl(1) = 1.8753
anl(2) = 4.6941
anl(3) = 7.8547
anl(4) = 10.9953

SIG1 = ( COS(anl(N))+COSH(anl(N)) ) / ( SIN(anl(N))+SINH(anl(N)) )

a(1) = 1.8753/U
a(2) = 4.6941/U
a(3) = 7.8547/U
a(4) = 10.9953/U

B = ((2*M-1)*PI*X)/(2*U)

FH21=(COSH(a(N)*X)-COS(a(N)*X)) - (SIG1 * (SINH(a(N)*X) - SIN(a(N)*X)))

FH22=SIN(B)

FH2 = CHORD*FH21*FH22

RETURN
END

!------------------------------------------------------------------------

SUBROUTINE HH4(N,M,X,U,YTHETA,CHORD,FH4)

PI=4.0*ATAN(1.0)

A = ((2*N-1)*PI*X)/(2*U)
B = ((2*M-1)*PI*X)/(2*U)

Y0 = (CHORD/2.0) - YTHETA

FH4 = ( (CHORD**2.0) * (Y0/CHORD-0.25) * (SIN(A) * SIN(B)) )

RETURN
END
 SUBROUTINE LL1(N,M,X,U,CHORD,FFL1)
DIMENSION anl(10),a(10)

PI=4.0*ATAN(1.0)

anl(1) = 1.8753
anl(2) = 4.6941
anl(3) = 7.8547
anl(4) = 10.9953

SIG1 = ( COS(anl(N)) + COSH(anl(N)) ) / ( SIN(anl(N)) +SINH(anl(N)) )

SIG2 = ( COS(anl(M)) + COSH(anl(M)) ) / ( SIN(anl(M)) +SINH(anl(M)) )

a(1) = 1.8753/U
a(2) = 4.6941/U
a(3) = 7.8547/U
a(4) = 10.9953/U

FL11= (COSH(a(N)*X)-COS(a(N)*X)) - (SIG1 * (SINH(a(N)*X) - SIN(a(N)*X)))

FL12= (COSH(a(M)*X)-COS(a(M)*X)) - (SIG2 * (SINH(a(M)*X) - SIN(a(M)*X)))

FFL1 = CHORD*(FL11*FL12)

RETURN
END
!------------------------------------------------------------------------
SUBROUTINE LL2(N,M,X,U,YTHETA,CHORD,FFL2)
DIMENSION anl(10),a(10)

PI=4.0*ATAN(1.0)

anl(1) = 1.8753
anl(2) = 4.6941
anl(3) = 7.8547
anl(4) = 10.9953

SIG1 = ( COS(anl(N)) + COSH(anl(N)) ) / ( SIN(anl(N)) +SINH(anl(N)) )

a(1) = 1.8753/U
a(2) = 4.6941/U
a(3) = 7.8547/U
a(4) = 10.9953/U

B = ((2*M-1)*PI*X)/(2*U)

FL21= (COSH(a(N)*X)-COS(a(N)*X)) - (SIG1 * (SINH(a(N)*X) - SIN(a(N)*X)))

FL22= SIN(B)

Y0 = (CHORD/2.0) - YTHETA

FFL2 = ( (CHORD**2) * (0.75-Y0/CHORD) * FL21*FL22 )

RETURN
END
 SUBROUTINE LL3(N,M,X,U,YTHETA,CHORD,FFL3)
DIMENSION anl(10),a(10)

PI=4.0*ATAN(1.0)

anl(1) = 1.8753
anl(2) = 4.6941
anl(3) = 7.8547
anl(4) = 10.9953

SIG2 = ( COS(anl(M)) + COSH(anl(M)) ) / ( SIN(anl(M)) +SINH(anl(M)) )

a(1) = 1.8753/U
a(2) = 4.6941/U
a(3) = 7.8547/U
a(4) = 10.9953/U

B = ((2*N-1)*PI*X)/(2*U)

FL31= SIN(B)
FL32= (COSH(a(M)*X)-COS(a(M)*X)) - (SIG2 * (SINH(a(M)*X) - SIN(a(M)*X)))

Y0 = (CHORD/2.0) - YTHETA

FFL3 = ( (CHORD**2) * (Y0/CHORD-0.25) * FL31*FL32 )

RETURN
END
!------------------------------------------------------------------------
SUBROUTINE LL4(N,M,X,U,DCL,YTHETA,CHORD,FFL4)

PI=4.0*ATAN(1.0)

A = ((2*N-1)*PI*X)/(2*U)
B = ((2*M-1)*PI*X)/(2*U)

Y0 = (CHORD/2.0) - YTHETA

FFL4 = ( (CHORD**3) * ((PI/8.0) - ( (Y0/CHORD-0.25) * (0.75-Y0/CHORD)*(2.0*PI) ))

* (SIN(A) * SIN(B)) )

RETURN
END
!------------------------------------------------------------------------
 APPENDIX I B

MATLAB code for the eigen value analysis

% kb = bending stiffness matrix

% mb = bending mass matrix
% kt = torsional stiffness matrix
% mt = torsional mass matrix
% mc coupled mass matrix
% h1,h2,h3,h4,l1,l2,l3,l4 = aerodynamic matrices
% the above matrices have to be input here for the eigen value analysis

kc = zeros(4,4)

% "H1 MATRIX"
h1 = zeros(4,4);

% "H3 MATRIX"
h3 = zeros(4,4);

% stiffness matrix
stiff1 = [kb;kc];
stiff2 = [kc;kt];
stiff = [stiff1,stiff2];
ss = stiff

% mass matrix
mass1 = [mb;mc];
mass2 = [mc;mt];
mass = [mass1,mass2];
mm = mass

% aerodynamic matrices

aeroh1 = [h1;h3];
aeroh2 = [h2;h4];
hh = [aeroh1,aeroh2];

aerol1 = [l1;l3];
aerol2 = [l2;l4];
ll = [aerol1,aerol2];
% n = number of modes
n=4

% finding eigenvalues for different speeds

u=-50.0;
for nn=1:17
u = u + 50.0
uu(nn) = u ;

a1 = zeros(2*n);
a2 = eye(2*n);
a3 = -( ss + (uu(nn)^2 * (hh)) );
a4 = -uu(nn)*(ll);

aa1 = [a1,a2];
aa2 = [a3,a4];
aa = [aa1;aa2];

b1 = eye(2*n);
b2 = zeros(2*n);
b3 = zeros(2*n);
b4 = mm;

bb1 = [b1,b2];
bb2 = [b3,b4];
bb = [bb1;bb2];

[v2,d2] = eig(aa,bb,'qz');

dd = (sort(diag(d2)))

end

% the eigen values obtained are complex .

% the eigen values obtained will not be in order which has to be sorted and the
graphs can be plotted using the plot command