16TH INTERNATIONAL CONFERENCE ON COMPOSITE MATERIALS

AEROELASTIC DIVERGENCE AND FREE VIBRATION OF

TAPERED COMPOSITE WINGS
Mengchun Yu*, Chyanbin Hwu*
* National Cheng Kung University, Department of Aeronautics and Astronautics

Keywords: aeroelastic, divergence, free vibration, wing, tapered, eigenvalue

Abstract Numerous different analytical model of the

composite wing structures such as the classical beam
In aeroelasticity, divergence of wing structures
model [2,3], the coupled bending-torsion model [4-
is an important static phenomenon. On the other
6], and the refined models taking warping restraint
hand, the free vibration analysis is significant to
[7-10], transverse shear deformation [11], shell
investigate the dynamic characteristics of the
bending strain [12,13], and cross-sectional materials
structure. Recently, a comprehensive wing model
and geometries [14,15] into account have been
considering the airfoil shape was established.
developed in the past several decades. However,
However, the analysis of tapered wing was still not
they just treated the wing as plane beam without
taken into account so far. In this paper, this model is
considering the airfoil shape due to the complexity.
extended to discuss the tapered composite wings. By
Recently, a comprehensive model considering the
employing the Hamilton’s principle and following
shape of airfoil for composite wing structures was
the standard procedure of finite element formulation,
therefore developed by our co-workers [16,17]. This
an elementwise comprehensive model was developed.
model can treat the more realistic wing structures
Because this model can be applied to both static and
such as stiffened multicell wings that are composed
dynamic analyses, divergence and free vibration will
of the laminated skins, stringers, ribs, and spars. No
be studied. In this model, aerodynamic force vector
matter how comprehensive the model is, it is still
containing the lift and moment of the wings was
restricted to the cases of the wings with uniform
approximated by the strip theory, which will then
airfoil shape along span direction.
lead to a standard eigen-relation for solving the
To study the tapered effect of the wing
divergence dynamic pressure by neglecting the
structure and to improve the computational
inertial terms. Besides, this model will also provide
efficiency of the comprehensive model, in this
another eigen-relation for solving the natural
article the elementwise comprehensive model is
frequencies when all the external forces are
derived from the Hamilton principle with the
eliminated.
concept of finite element method embedded in. Both
the accuracy and computational efficiency of the
elementwise comprehensive model are compared
1 Introduction with the results that are provided by the existing
The fundamental work concerning the papers or the commercial finite element software
divergence instability of swept metallic wings was packages (ANSYS) in the aeroelastic divergence and
done about fifty years ago. It was shown that low free vibration analyses.
static aeroelastic divergence speeds were associated 2 Elementwise Comprehensive Model of
with the swept-forward wings unless they were Composite Wing
stiffened enough. Because the metallic wing is
limited to its material properties, the swept-forward 2.1 Assumptions and the Comprehensive Model
wing aircraft was considered as an impossible task If the cover skin of the wing is made of
for a time. Until the aeroelastic tailoring concept for composite laminates, the entire wing structure may
the composite wing structures was raised and studied be simulated as a composite sandwich plate [16,17]
by Krone [1], relevant research about swept in which the wing skins, stringers, and spar flanges
composite wing sprouts once again. are simulated as the faces resisting the in-plane force
and bending moment; furthermore, the spar webs
1
 Mengchun Yu, Chyanbin Hwu

and ribs are simulated as the core resisting the boundary conditions can be expressed in matrix
transverse shear and the force normal to the face. form as follows:
With the transverse stiffening components, it is

F ′ − F0 + p = I 0 d (2)
usually assumed that the wing chordwise section is
strong enough to avoid the deformation. So that
some aerodynamic parameters are constants as the F = K 1d + K 2 d ′ , F0 = K 0 d + K 1T d ′ (3)
airfoil shape of the wing is fixed. Besides, because
the wing cross section must have a streamline shape
commonly referred to as an airfoil section, the K 1d + K 2d′ = Fˆ on y = y f and y = y0 (4)
thickness of the sandwich will be a function of the
where F, p and d are, respectively, the section force,
airfoil and the thickness is not too small to neglect
surface force and displacement vectors; K0, K1 and
the transverse shear deformation. Based upon these
K2 are stiffness matrices related to the extensional,
considerations shown above, a comprehensive
coupling and bending stiffness of the composites; I0
model was developed [17]. In this model, the
is a matrix related to the mass, gravity center and
displacement field can be expressed as follows:
moment of inertia. The overdot and the prime denote,
z respectively, the differentiation with respect to time
and space. The overhat stands for the prescribed
value on the boundaries. Please refer to Hwu and
Gai [17] for detailed explanation.
2.2 Deduction of the Elementwise Comprehensive
x Model
In order to refine the comprehensive model, the
model is rederived from Hamilton’s principle [18] as
shown by the following equation:
y t2
(π − κ )dt = 0
Fig. 1. Geometry of the composite wing [16]
δ ∫t1
(5a)

where π is potential energy consisting of strain

u ( x , y , z , t ) = zθ ( y , t ) (1a) energy and the work done by the body forces and
surface tractions and κ is the kinetic energy as
follows,
{ }
v ( x, y, z , t ) = v0 ( y, t ) + z β f ( y, t ) + x β r ( y, t ) (1b)

∫ (Tˆ u )dS
⎛1 ⎞
w(x, y, z , t ) = w f ( y, t ) − xθ ( y, t )
π= ∫ ⎜⎝ 2 σ
V
ij ε ij − f i ui ⎟dV −
⎠ Sσ
i i (5b)
(1c)
where u, v, and w are the displacement components
∫
1
κ= ρu i u i dV (5c)
in the direction x (chordwise), y (spanwise), z 2 V
(thicknesswise), respectively (Fig. 1); t is the time
variable, v0 the midplane displacement in the y By introducing Eq. (3) into Eq. (5b) and (5c) and
direction, wf the displacement (positive upward) in subtracting Eq. (5c) from Eq. (5b), we may obtain
the z direction measured along the reference axis, π −κ

∫ [d′ K d′ + 2d K d′ + d K d − d I d − 2p d]dy
and θ the rotation angle with respect to x axis due to =
1 T T T T T
2 1 0 0
twist around the reference axis (positive nose up), 2 y (6)
that is, βx=θ. Furthermore, βf is the rotation angle of −F d[ ]
ˆT yf
y0
x-z plane with respect to the y axis measured at the
reference axis and βr is the variation rate of βf along The governing equations and boundary conditions
the x axis, hence βx=βf+xβr. Through the above derived from Eq. (6) agree with Eq. (2) and Eq. (4)
interpretation, it is obvious that there are five basic derived in [17].
deformation functions for the stiffened composite
wing structures as follows, v0, wf, θ, βf, and βr. The
equations of motion, constitutive equation, and

2
 AEROELASTIC DIVERGENCE AND FREE VIBRATION OF TAPERED
COMPOSITE WINGS

ye Le 1
qn ≡ ρV n2 = q cos 2 Λ (10)
2
Node 1 Node 2 Node 3 where ρ , V n , q , and Λ are the density of the
airflow, airflow velocity normal to the leading edge,
Fig. 2. Local coordinate of each element dynamic pressure, and the angle of sweep
respectively. The lift curve slope coefficient a is a
Because y-axis is the primary coordinate aerodynamic parameter defined as
variable in the wing structure. Thus the wing
structures are simulated as a beam divided in n dC L AR
a≡ = a0 (11)
elements connected along the y axis. Each element dα AR + 4 cos Λ
has 3 nodes (Fig. 2) and the shape function of each
where a0 is the corresponding 2D lift curve slope and
node is expressed by 2nd order polynomials. By
AR is the aspect ratio.
following the procedures stated between eqns.(5)-(6),
With the aerodynamic functions listed above,
the energy expression of each element can be written
the external aerodynamic force vector fe can be
in terms of the nodal displacement vector ue as
simplified further as
π −κ = ∫
⎧1 T 2
[ ]T ⎫
⎨ u e g K e − g J e u e − g (t )u e f e ⎬dy f e = N T ( y )p( y )dy = q n (p 0 + K a u e )
∫
2
(7)
ye ⎩2 ⎭ (12)
in which g is a function of time, and Ke, Je and fe are,
where the value of p 0 is constant, but the value of
respectively, element stiffness matrix, element
inertia matrix and element external force vector the aerodynamic stiffness matrix K a varies with the
determined through the following relations: interaction between the aerodynamics (L, T) and the
displacement function (ue).
∫ {N′ K N′ + 2N K }
T T T
Ke = 2 1 N′ + NT K 0 N dy (8a) Because divergence is a static phenomenon in
ye
aeroelasticity, we can neglect the time effect by
eliminating the inertial terms Je. Equation (7) can

∫
1 therefore be rewritten as follows
Je = NT I 0 Ndy (8b)
2 ye
π −κ =
1 T
2
(
u e K e u e − u Te q n p 0 + K a u e ) (13)

fe =
1
2 ye∫
NT pdy + NT Fˆ [ ] yf
y0 (8c) After substituting Eq. (13) into Eq. (5a), the
equilibrium equation of the divergence analysis is
where N(y) is the shape function matrix. expressed as
2.2.1 Divergence
In this article, we will use the aerodynamic
⎡1
( T
) ( T
⎢ 2 K e + K e − qn K a + K a )⎤⎥u e = qnp 0 (14)
⎣ ⎦
strip theory and the known results for two-
dimensional flow to approximate the lift and the The determinant of the coefficient matrix of ue
pitching moment [19]. The relation between the should be zero to avoid the trivial solutions, which
structure deformation and aerodynamic forces (lift will then lead to the result of divergence,
force L and pitching moment T) can be expressed as
[ (
L( y ) = q n acθ 0 + ac θ ( y ) − w′f tan Λ )] (9a)
1
2
( ) (
K e + K Te − qn K a + K Ta = 0 ) (15)

[( ) ( )]
Obviously, it is an eigen problem and the eigen-
T ( y ) = qn aceθ0 + c 2Cm, ac + ac eθ ( y ) − ew′f tan Λ (9b) value and eigen-vector indicate the normal
divergence dynamic pressure and the nodal
where c is the chord length, θ 0 the initial angle of displacement vector respectively. After the
attack, Cm,ac the pitching moment coefficient about superposition of all the elements accomplishing, the
the aerodynamic center, and e the distance between eigen problem can be solved.
the lines of aerodynamic and flexural centers. The
normal dynamic pressure qn is

3
 Mengchun Yu, Chyanbin Hwu

2.2.2 Free Vibration As shown in Fig. 3, even the wing is just

To determine the natural frequency of the stiffened divided into few elements, the divergence dynamic
wing structures, the values of all components in the pressures (of different aspect ratio) calculated by the
body force distribution vector p~ and the prescribed elementwise comprehensive model still converge
very rapidly . To further examine the computational
surface traction Tˆi are both set as zero in the
efficiency, the example of composite wing structure
equation. So that the element external force vector with NACA 2412 airfoil used by Hwu and Tsai [16]
can be expressed as fe=0. Furthermore, g(t) is is adopted here. The present method only spends
assumed as a harmonic motion with the natural 1/10 of the time consumed in Ref. [16].
frequency ω in the free vibration analysis. By using Next, the accuracy of present model is checked
Eq. (7) and Hamilton’s principle again, the equation by comparing with the results of the uniform
of motion is derived as metallic wing structure (without tapered) evaluated
1
2
[( ) (
K e + K Te + ω 2 J e + J Te u e = 0 )] (16)
by Weisshaar [20] and Librescu [21]. As shown in
Fig. 4, the present data also agree with them
especially when AR = 4. Besides, Fig. 4 also reveals
The nontrivial solutions of Eq. (16) exist only when that the divergence dynamic pressure will decrease
the determinant of the coefficient matrix of ue as the values of aspect ratio AR increase.
becomes zero, which will then provide us the
following equation for solving natural frequencies:
1
2
(K e )
+ K eT + ω 2 1
2
(M e )
+ M eT = 0 (17)

As shown above, it is also an eigen problem

and the eigen-value and eigen-vector indicate the
square of natural frequency and the nodal
displacement vector respectively. After the
superposition of all the elements accomplishing, the
eigen problem can be solved. It is similar to the
divergence analysis.
3 Numerical Results

3.1 Divergence
In this section, the calculating efficiency of the Fig. 4. Comparison with the existing data
elementwise comprehensive model will be shown
first. And the accuracy of the present model will be
examined by comparing with the existing numerical
solutions provided by other papers.

Fig. 5. Tapered effect in aeroelastic divergence of

the wing structure

Fig. 3. Convergence study of the present model

4
 AEROELASTIC DIVERGENCE AND FREE VIBRATION OF TAPERED
COMPOSITE WINGS

Before investigating the tapered effect in wing Table 1 shows the natural frequencies got from
structures, the taper ratio rt is denoted as the present model in this article well agree with the
solution of Hwu & Gai and ANSYS. It deserves to
be mentioned that the calculating efficiency is also
croot − ctip improved as shown in Table 2.
rt = (18)
croot
Table 1. Natural frequencies comparison for the
where croot and ctip are the chord lengths of the wing wing structure with uniform span (Unit:Hz)
root and wing tip, respectively. Mode No. Present
Hwu & Gai
ANSYS
The examples in a series of tapered ratios are [17]
solved for AR = 4 and 5. The numerical results are I 16.89 16.15 15.75
II 99.12 96.45 93.84
shown in Fig. 5. It shows that once the value of III 108.69 110.34 114.94
tapered ratio increases, the value of divergence IV 256.83 252.45 245.18
pressure will increase, too. V 325.23 333.85 343.82
3.2 Free Vibration
Table 2. Time consumed comparison
In this section, the present model is first Hwu & Gai
Present ANSYS
examined by comparison with the existing numerical [17]
data of the wing structure with uniform span Time consumed
2 25 30
(minute)
evaluated by Hwu and Gai [17]. The geometry
properties of the wing are chordwise length c=0.1m
and spanwise length L=0.4m shown in Fig. 6 (with Table 3. Effects of the taper ratio on the basic
croot = ctip ) where croot and ctip means the chord natural frequencies (Mode I) (Unit:Hz)
Taper ratio Mode No. Present ANSYS
length at the wing root and wing tip respectively. I 16.89 15.25
The detailed material properties of wing are listed in II 99.12 90.35
[17]. 1.0 III 108.69 113.16
IV 256.83 234.64
ctip V 325.23 275.61
I 17.38 14.37
II 99.88 85.11
0.95 III 113.85 117.66
IV 256.83 224.07
V 330.82 280.47
Wing Tip I 19.08 16.47
II 101.88 88.43
0.80 III 132.73 125.47
IV 256.89 224.35
V 350.35 277.64
y I 20.49 17.13
L=400mm II 103.39 87.02
0.70 III 148.41 133.02
IV 257.09 217.85
V 366.40 271.96
I 22.21 17.91
Wing Root II 105.23 85.39
0.60 III 167.08 155.67
IV 257.58 210.67
V 385.92 268.10
I 35.50 22.90
II 122.80 80.28
0.20 III 267.76 120.55
x IV 282.12 180.53
croot V 465.98 213.95

Fig. 6. Geometry of the tapered composite wing

5
 Mengchun Yu, Chyanbin Hwu

[2] Bisplinghoff, R. L., Asheley, H., and Halfman, R. L.,

“Aeroelasticity”. Addison-Wesley, Cambridge, MA,
Chaps. 2 and 3, 1955.
[3] Megson, T. H. G., “Aircraft structures-for
engineering students”. 2nd ed., Edward Arnold,
London, Chap. 8, 1990.
[4] Weisshaar, T. A. , and Foist, B. L. , “Vibration
tailoring of advanced composite lifting surfaces,”
Journal of Aircraft, Vol. 22, pp.141-147, 1985.
[5] Chandra, R., Stemple, A. D., and Chopra, I., “Thin-
walled composite beams under bending, torsional,
and extensional loads,” Journal of Aircraft, Vol. 27,
pp.619-626, 1990.
[6] Banerjee, J. R., and Williams, F. W., “Free vibration
of composite beams: an exact method using symbolic
computation,” Journal of Aircraft, Vol. 32, No. 3,
Fig. 7. The trend of the natural frequency (Mode I) pp.636-642, 1995.
[7] Crawley, E. F., and Dugundji, J., “Frequency
Then we consider a tapered composite wing determination and nondimensionalization for
structure with NACA 2412 airfoil and the same composite cantilever plates,” Journal of Sound and
material properties as the preceding example. But Vibration, Vol. 72, No. 1, pp. 1-10, 1980.
the values of taper ratio are set as a series of values [8] Lottati, I., “Flutter and divergence aeroelastic
1.00, 0.95, 0.80, 0.70, 0.60, and 0.20. In order to characteristics for composite forward swept
cantilevered wing,” Journal of Aircraft, Vol. 22,
check the accuracy, the natural frequencies are
No.11, pp. 1-10, 1985.
evaluated in two ways, the present model and
[9] Librescu, L., and Khdeir, A.A., “Aeroelastic
ANSYS. The results of the present model do not
divergence of swept-forward composite wings
agree well with ANSYS shown in Table 3. The including warping restraint effect,” AIAA Journal,
difference may result from that the mesh-process Vol. 26, No. 11, pp. 1373-1377, 1988.
encounter something trouble for the geometry [10] Oyibo, G. A. and Bentson, J., “Exact solutions to
variation (tapered wing). However, we observe a oscillations of composite aircraft wings with warping
trend that the basic natural frequency will increase constraint and elastic coupling” AIAA, Vol.28, No. 6,
when the taper ratio rt decrease, as shown in Fig. 7. pp.1075-1081, 1990.
[11] Librescu, L., and Song, O., “On the static aeroelastic
4 Conclusion
tailoring of composite aircraft swept wings modelled
By embedding the finite element formulation as thin-walled beams structures,” Composite
in the comprehensive model [17], the elementwise Engineering, Special Issue, Use of Composite in
model in this article can deal generally with the Rotor Craft and Smart Structures, Vol. 2, No. 5-7, pp.
uniform and tapered composite wing structures. 497-512, 1992.
Besides, the excellent calculating efficiency saves a [12] Volovoi, V. V., and Hodges, D. H., “Single- and
lot of computational time that is consumed before. multicelled composite thin-walled beams,” AIAA
Journal, Vol. 40, No. 5, pp. 960-965, 2002.
The divergence analysis and the free vibration
analysis both are eigen-value problems. With the [13] Volovoi, V. V., and Hodges, D. H., “Theory of
anisotropic thin-walled beams,” Journal of Applied
usual procedure for finite element formulation, the
Mechanics, Vol. 67, No. 3, pp. 453-459, 2000.
eigenvalue such as divergence dynamic pressure and
[14] Yu, W., Volovoi, V. V., Hodges, D. H., and Hong, X.,
the natural frequencies could be solved similarly. It
“Validation of the variational asymptotic beam
can also be observed that the divergence dynamic sectional analysis,” AIAA Journal, Vol. 40, No. 10,
pressure will increase but the fundamental natural pp. 2105-2112, 2002.
frequency will decrease when the taper ratio rt [15] Yu, W. , Hodges, D. H., Volovoi, V. V., and Cesnik,
increases. C.E.S., “On Timoshenko-like modeling of initially
References curved and twisted composite beams,” International
Journal of Solids and Structures, Vol. 39, pp. 5101-
[1] Krone, Jr. N. J. “Divergence elimination with 5121, 2002.
advanced composites". AIAA Paper 75-1009,
L.A.,CA. , 1975

6
 AEROELASTIC DIVERGENCE AND FREE VIBRATION OF TAPERED
COMPOSITE WINGS

[16] Hwu, C. and Tsai, “Aeroelastic divergence of

stiffened composite multicell wing structures,”
Journal of Aircraft, Vol.39, No. 2, pp.242-251, 2002.
[17] Hwu, C. and H.S. Gai, “Vibration analysis of
composite wing structures by a matrix form
comprehensive model,” AIAA Journal, Vol.41, No.11,
pp.2261-2273, 1997.
[18] Clough, R.W. and Penzien, J., “Dynamics of
structures”, McGraw-Hill, New York, 1993.
[19] Oyibo, G. A., “Generic approach to determine
optimum aeroelastic characteristics for composite
forward-swept-wing aircraft,” AIAA Journal, Vol.22,
No.1, pp.117-123, 1984.
[20] Weisshaar, T. A., “Divergence of forward swept
composite wings,” Journal of Aircraft, Vol.17, No. 6,
pp.442-448, 1980.
[21] Librescu, L., and Simovich J., “General formulation
for the aeroelastic divergence of composite swept-
forward wing structures,” Journal of Aircraft, Vol.25,
No. 4, pp.364-371, 1988.