FINITE ELEMENT-BASED FLUTTER ANALYSIS OF

CABLE-SUSPENDED BRIDGES
By Ahmad Namini, 1 Associate Member, ASCE,
Pedro Albrecht,2 Member, ASCE, and Harold Bosch3
Downloaded from ascelibrary.org by University of Tasmania Library on 03/22/19. Copyright ASCE. For personal use only; all rights reserved.

ABSTRACT: In this paper, an approach to computational flutter analysis is pre-

sented, which permits the determination of the critical wind velocity that initiates
damping- and stiffness-driven flutter of cable-suspended bridges. Dynamic response
may be coupled elastically and/or aerodynamically, with the developed algorithm,
the so-called pK-F method determining the preflutter and postflutter responses by
solution to the modal equations of motion. The entire method is presented in matrix
form, so as to be easily implemented into finite element systems. As an example,
the pK-F method is applied to the Luling cable-stayed bridge, which was experi-
mentally tested for flutter in the wind tunnel. Also, examples are provided for
flutter determination with different deck sections. The pK-F method has proved
reliable in its methodology and efficient in its use. The generality of the method
permits more flutter scenarios to be examined, including construction stages and
different flutter derivatives for different portions of the bridge structure.

INTRODUCTION

Since the famous Tacoma Narrows Bridge failure of 1940, the bridge
engineering community has been faced with the design consideration of
aerodynamic flutter. Flutter is the dynamic instability that develops from
the mutual interaction of elastic, inertial, damping, and self-excited aero-
dynamic forces, whereby at some critical wind velocity the bridge oscillates
in a divergent, destructive manner. Flutter design consists of analyzing a
proposed bridge deck configuration for determination of the lowest wind
velocity that initiates instability. The wind velocity at flutter should be higher
than meteorological possible wind velocities at the bridge site.
At present, flutter analysis combines both experimental and analytical
procedures. The pioneering free-oscillation method (Scanlan and Tomko
1971) represented self-excited lift^L,,, and drag, Dp, forces as well as the
aerodynamic moment, Ma, by so-called flutter derivatives. Flutter deriva-
tives relate forces, per unit span, to the vertical bending or heave, h, lateral
or sway, p, and torsional, a, deformations and associated velocities as

L„ = \ 9U\2B) KHf(K) | + KHUK) ^ + K2Hl(K)a (1)

Dp = ± PU2(2B) KP*(K) | + KP$(K) JJ- + K2P;(K)a (2)

'Asst. Prof., Dept. of Civ. and Arch. Engrg., Univ. of Miami, Coral Gables, FL
33124.
2
Prof., Dept. of Civ. Engrg., Univ. of Maryland, College Park, M D 20742.
3
Res. Struct. Engr., F H W A Turner-Fairbank Highway Res. Ctr., 6300 George-
town Pike, McLean, V A 22101.
Note. Discussion open until November 1, 1992. To extend the closing date one
month, a written request must be filed with the ASCE Manager of Journals. The
manuscript for this paper was submitted for review and possible publication on June
16, 1990. This paper is part of the Journal of Structural Engineering, Vol. 118, No.
6, June, 1992. © A S C E , ISSN 0733-9445/92/0006-1509/$1.00 + $.15 per page. Paper
No. 26629.

1509

J. Struct. Eng., 1992, 118(6): 1509-1526

 2
Ma = \ PU*(2B ) KAf(K) ^+ KA$(K)~ + K2A${K)K (3)

where H*, P*, A* = nondimensional flutter derivatives, i = 1, 2, 3,

dependent upon K\ p = air mass density; U = wind velocity; K = reduced
Downloaded from ascelibrary.org by University of Tasmania Library on 03/22/19. Copyright ASCE. For personal use only; all rights reserved.

frequency of response (K = Bu>IU); B = deck width; and u> = response

circular frequency. These flutter derivatives are determined in a controlled
wind tunnel environment for a representative portion of the bridge deck
over a wide range of wind velocities with the experimental techniques for
their extraction described thoroughly by Scanlan (1975). Flutter derivatives
are typically plotted as a function of the reduced velocity, defined as
U/(NB), which is 2^/K, where N = response frequency in Hz (cps). Note
that in (l)-(3), the lift and drag forces are coupled aerodynamically with
the moment but not to each other.
In order to determine the critical wind velocity that initiates flutter, modal
analysis is typically employed, with natural modes associated with predom-
inantly torsional motion being the most critical. For cable-suspended bridges,
torsional and sway motions are typically coupled due to changes in bridge
elevation and/or an appreciable vertical offset between the deck cross sec-
tion's center of gravity and elastic axis. Computationally, certain simplifi-
cations are employed in a modal analysis that represents a worst-case sce-
nario in flutter generation. Scanlan (1987) describes a scenario where all,
or most, of the wind energy feeds into a single mode, or a very restricted
pair of modes, that is unstable. Therefore, only modes suspected of being
flutter prone need be accounted for in the modal representation of overall
bridge response. Another simplification is that structural frequencies are
usually not strongly modified by aerodynamic forces. Therefore, classical
aircraft-type flutter, referred to also as stiffness-driven flutter, in which two
modes (usually bending and torsional) coalesce to a single flutter frequency,
will not occur, and is justified since its occurrence is rare. However, Scanlan
(1978, 1981) presented a modal energy and root finding approach to the
stiffness-driven flutter solution.
In damping-driven flutter, all flutter derivatives associated with aerody-
namic coupling (i.e., H*, H*, P*, P*,A\) are usually neglected, which
aerodynamically, but not elastically, uncouples the equations of motion.
Instead, single-degree-of-freedom flutter, referred to also as damping-driven
flutter, is determined. This form physically interprets the onset of flutter as
the wind velocity that produces enough negative aerodynamic damping to
offset the bridge's own mechanical damping. Computationally, but not phys-
ically, damping-driven flutter may be viewed as a special case of stiffness-
driven flutter, in which all aerodynamic coupling terms are neglected.
This paper presents a general numerical procedure, the so-called pK-F
method, which is modeled after the aircraft industry's p-K method (Hassig
1971) for the solution of preflutter and postflutter response of airfoils. Aero-
dynamic coupling is employed through the use of all flutter derivatives of
(l)-(3). To facilitate a general-purpose algorithm suitable for computer
implementation, the numerical procedure is presented in matrix form uti-
lizing a space frame idealization, which can be incorporated in finite element
code for cable-suspended bridge analyses (see Fig. 1 for the positive direc-
tions of the 12 degrees of freedom). Finally, as an example of the pK-F
method's usefulness, it is applied to the Luling cable-stayed bridge, which
was experimentally tested for flutter (Bosch 1977).

1510

J. Struct. Eng., 1992, 118(6): 1509-1526

Downloaded from ascelibrary.org by University of Tasmania Library on 03/22/19. Copyright ASCE. For personal use only; all rights reserved.

FIG. 1. Positive Directions of 12-Degree-of-Freedom Space Frame Member

MEMBER AERODYNAMIC MATRICES

Eqs. (l)-(3) represent aerodynamic forces per unit of deck span. To

convert these uniformly distributed forces into member end effects, a simple
lumping procedure is adopted whereby one-half of the load is assumed to
act at each member end [See Figs. 2(a-c)]. Using a space frame element
discretization of the overall bridge, the ith member's end aerodynamic joint
loads, A,-, is subdivided into stiffness, A„-, and damping, Adh components
as

A,- = ^ PU2 (A*/q< + JjKAt (4)

with q, and q, being space frame deformation and velocity nodal vectors,
respectively. With L being the member length, the nonzero terms of the
aerodynamic stiffness matrix of order 12 are
4,(2,4) = 4,(8,10) = -BLK2H% (5a)
4,(3,4) = 4,(9,10) = BLK^n (56)
4,(4,4) = 4,(10,10) = B2LK2A$ (5c)
and the nonzero terms of the aerodynamic damping matrix of order 12 are
4 d (2,2) = 4,(8,8) = BLKHJ ; . . . (6a)
4,(2,4) = 4,(8,10) = -B^LKHi (6b)
4,(3,3) = 4,(9,9) = BLKPt (6c)
4,(3,4) = 4,(9,10) = B^LKFt (6d)
2
4,(4,2) = 4,(10,8) = -B LKA? (6e)
Ad^A) = ^(10,10) = B*LKAi - (6/)
Since aerodynamic forces are nonconservative, the stiffness and damping
1511

J. Struct. Eng., 1992, 118(6): 1509-1526

Downloaded from ascelibrary.org by University of Tasmania Library on 03/22/19. Copyright ASCE. For personal use only; all rights reserved.

LIFTING FORCE PER UMIT LENGTH LUMPED LIFTING FORCE

(a)

MOMENT PER UNIT LENGTH LUMPED MOMENT

(c)
FIG. 2. (a) Lumped Aerodynamic Lift Force; (b) Lumped Aerodynamic Drag Force;
and (c) Lumped Aerodynamic Moment

components of (4) need not be symmetrical as is the case. Also, the member
aerodynamic matrix A,-, continually changes over a range of wind velocities
and response frequencies. In this implementation of the pK-F method,
member aerodynamic matrices, A,-, are determined at up to eight unique
reduced velocities, [U/(NB)]h 1, 2, 3, . . . , 8. The eight values of [U/(NB)]
should represent reduced velocities for which flutter derivatives are known
from experimental investigations.
1512

J. Struct. Eng., 1992, 118(6): 1509-1526

 FLUTTER EQUATIONS OF MOTION

Using either a lumped- or consistent-mass formulation and corresponding

elastic stiffness formulation (which may or may not account for elastic cou-
pling), mass, M,, and stiffness, Kf, matrices are written for each element.
Though proportional viscous damping will shortly be assumed, for com-
Downloaded from ascelibrary.org by University of Tasmania Library on 03/22/19. Copyright ASCE. For personal use only; all rights reserved.

pleteness, assume that member damping, D,, matrices exist. Thus, by map-
ping all matrices onto the global axes, the assembled equilibrium equations
are given as

Mq + Dq + Kq = ^ pC/2 U , q + ^ A d qj (7)

The modal analysis approach for dynamic motion states that the response
at all degrees of freedom is separated into spatial (natural modes) and time-
dependent (generalized coordinates) components as

q=n (8)
where <& = an rc-row by m-column matrix of natural mode shapes (r'th
column corresponding to r'th natural mode); £ = an m-row vector of gen-
eralized coordinates (ith entry corresponding to ith generalized coordinate);
and n and m = the total number of degrees of freedom and total number
of modes permitted to participate in the response, respectively. Note that
natural modes may be any desired modes thought to be susceptible to flutter,
either coupled or uncoupled elastically.
Eq. (8) is substituted into (7), and premultiplied by the transpose of the
natural modes matrix to obtain, as a result of the orthogonality between
modes m-equations of flutter response,

M«g + D«g + K*| = \ pt/2 f Af| + i ASIJ (9)

The generalized mass, Ms, damping (assuming proportional damping be-

tween modes), Dg, and stiffness, Kg, are diagonal matrices, and, thus, they
are uncoupled between modes. However, the generalized aerodynamic ma-
trices, Af and Af, are nondiagonal, and flutter response is thereby coupled.
Since flutter response is typically harmonic, it is assumed that the gen-
eralized coordinate is of a damped, harmonic form
g = R exp(jjrf) (10)
where R = response amplitude; JJL = (8 + «)co; 8 = logarithmic decrement;
and co = circular frequency. Eq. (10) can be rewritten in terms of the reduced
frequency, K = asBIU, as

£ = R exp f | t f ( 8 + i)t\ (11)

and a complex variable, p, is defined to be

p = K{h + i) (12)
Substituting (12) and (11) into (9), yields equations of motion as
1513

J. Struct. Eng., 1992, 118(6): 1509-1526

 MM§) P2 + DMf)p + *«

IpirUi + |ASP Rexp(|/rt) = 0 (13)

Downloaded from ascelibrary.org by University of Tasmania Library on 03/22/19. Copyright ASCE. For personal use only; all rights reserved.

The aerodynamic forces as reflected by flutter derivatives are experimentally

derived. Damping forces are accounted for by the generalized damping
matrix, D g . Therefore, in (13), the complex variable, p, acting on the gen-
eralized aerodynamic damping matrix, is assumed to only account for the
undamped response, i.e., the imaginary part of p that corresponds to the
term iK of (12). Therefore, (13) is rewritten to its final form as

M* ( ! ) p2 + ,D* ( | ) p + K*

\plpUl + \k#K Rexp(!^) =0 (14)

For a nontrivial solution to exist, the response amplitude in (14) cannot be

zero, and the exponential is never zero. Therefore, the determinant of the
remaining expression is zero and is written as

F = M* ( ! ) p2 + V* ( | ) p + K* - \ PU2 (Af + ! A|iX

(15)
where F = the flutter determinant, which is a function of complex variable
p. Therefore, for a given wind velocity, U, and reduced frequency, K, the
value, noting real and imaginary parts, of p that makes the flutter deter-
minant equal zero gives the actual oscillatory response. The logarithmic
decrement can then be computed as

5 = ^ (16)
Im(p)
and the response frequency as

<o = f Im(P) (17)

where Re and Im = the real and imaginary parts of a complex variable,

respectively. Depending on the sign of the logarithmic decrement, the re-
sponse is defined to be
8 < 0, stable (18a)
8 = 0, neutrally stable (18b)
8 > 0, unstable (18c)
The wind velocity that produces the neutrally stable condition is termed the
flutter velocity, Uf, with the corresponding flutter frequency, u/.
1514

J. Struct. Eng., 1992, 118(6): 1509-1526

 pK-F METHOD

Computationally, the pK-F method provides an algorithm to solve for

the flutter velocity, Uf, of (15). First, one must choose a range of wind
velocities, from t/low to f/high, in increments of Umc for which to check for
flutter susceptibility. The steps in the solution technique, the pK-F method,
Downloaded from ascelibrary.org by University of Tasmania Library on 03/22/19. Copyright ASCE. For personal use only; all rights reserved.

are as follows:

1. Compute generalized mass, damping, and stiffness matrices.

2. Compute a series of generalized aerodynamics matrices for a set of pre-
defined reduced frequencies.
3. Compute two starting values, pu and p,,>2, for each natural mode in the
flutter equations, with subscript i corresponding to the /th natural mode.
4. Compute the current wind velocity, Ucm, from within the range of wind
velocities t/,ow to I/high, starting with t/,ow, ending with Uhigh, and incrementing
uinc.
5. Loop over the /th natural mode, / = 1, 2, 3, . . . , m where m is the
total number of natural modes.
6. Determine the reduced frequencies corresponding to the current wind
velocity and response frequency, Kx and K2, where Kt = lm(pu) and K2 =
7. Determine the generalized aerodynamic matrices corresponding to Kx
and K2; interpolate or extrapolate if necessary.
8. Form the flutter equation, (15), for (p,,, K{) and (p, 2 , K2).
9. Determine the flutter determinants corresponding to (p,tl, Kx) and (p,a,
K2), as Ft and F2, respectively.
10. Using a linear prediction scheme, choose the next value of p as

P,, = ' " % : p;f> d9)

11. Repeat steps 6-11 until F2 < tolerance, where the tolerance value is
nearly zero. [Note that the next Kx corresponds to the previous Im(p,2) and
the next K2 corresponds to the previous lm(p,}).]
12. When the flutter determinant is below the tolerance, the true oscillatory
response, pic, of a given natural mode and wind velocity is found. Therefore,
the logarithmic decrement, (16), and frequency, (17), are then computed.
13. Repeat steps 5-12 for all natural modes.
14. Establish starting values of plx and pj2 for the next wind velocity and
each natural mode by prorating the true response pic for / = 1, 2, 3, . . . , m
as

(20)
p"=^irnr
p,,2 = 1.01 p , ^ (21)

15. Repeat steps 4-14 for all wind velocities in the predefined range.

The values of C/,ow, Uhigh, and Uinc, must be selected with care. With regard
to the lowest wind velocity, C/,ow, for the algorithm to be stable, starting
values of p in step 3 should be determined from a nearly undamped (8 ~
1515

J. Struct. Eng., 1992, 118(6): 1509-1526

 0) case. Physically, an undamped case corresponds to aerodynamic forces
generated at extremely low wind velocities. Therefore, U,ow should be nearly,
but not exactly, zero since a zero velocity produces no force, and thereby
would not start the iterative process of steps 6-12. Good starting values of
p are
Downloaded from ascelibrary.org by University of Tasmania Library on 03/22/19. Copyright ASCE. For personal use only; all rights reserved.

p,,x = -0.015,. + iS, (22)

Pi.2 = iS, (23)
where St = (co,5)/[/low is the reduced frequency corresponding to the ith
natural mode, which, in fact, is a natural mode's Strouhal Number. Also,
(22) uses Theodorsen's (1934) circulation constants at a reduced frequency
of zero for the other starting value of p.
The upper bounds on the wind velocity can be as large as desirable, but
should be restricted to meteorologically possible values. The increment of
wind velocity must be small so that the prorating in step 14 of the response
from one wind velocity to the next should not corrupt preflutter results.
The physical interpretation of (20) and (21) is that for small increments of
wind velocity, proportional aerodynamic forces exist from one wind velocity
to the next. This is justified only when the incremental wind velocity is
small.
In step 7, the generalized aerodynamic matrices are continually changing
for a given reduced frequency. Therefore, care must be taken in forming a
representative set of flutter derivatives over a wide range of reduced fre-
quencies. As stated earlier, the pK-F method, as implemented into a pro-
duction program, has a maximum of eight unique entries, which encompass
the entire range of experimental wind tunnel reduced frequencies, with
linear interpolation and extrapolation.
In step 10, two known approximations of roots and corresponding flutter
determinant values are shown. The next approximation of the root is linearly
extrapolated, (19) (see Fig. 3). The linear prediction scheme is numerically
stable for small-order matrices (Bathe 1982). Flutter equations are typically
formed for only a few natural modes.

LULING CABLE-STAYED BRIDGE

The 837.6-m long, four-lane Luling cable-stayed bridge on Interstate 310

spans the Mississippi River between the towns of Luling and Destrehan in
F(p) k

FIG. 3. Linear Root Predictor

1516

J. Struct. Eng., 1992, 118(6): 1509-1526

 St. Charles Parish, about 20 km west of New Orleans, Louisiana, which
opened to traffic in October 1983. The bridge has five spans: a 372.5-m
center span between the towers, two 154.8-m and 150.9-m anchor spans,
and two 79.2-m approach spans. The cables are arranged in a double-plane
fan configuration with 12 cables in each plane.
Downloaded from ascelibrary.org by University of Tasmania Library on 03/22/19. Copyright ASCE. For personal use only; all rights reserved.

The design engineers evaluated 15 deck sections and several railing types
for aerodynamic stability. Many were paired sections, meaning that sections
of same configuration had either an asphalt or a concrete roadway. The
sections consisted of single-box, double-box, and open-truss configurations.
The railings were of the open-slotted perpendicular, tapered closed-slotted,
and tapered with a welded fascia plate types.
Models of the 15 proposed deck sections were built to a 1/60 scale and
then tested in the wind tunnel at various angles of incidence (Bosch 1977).
The models were 1.52 m long, corresponding to a 91.5-m length of the
bridge deck, and 0.43 and 0.51 m wide. After all models had been tested
in the wind tunnel over a wide range of wind velocities, the final section,
shown in Fig. 4, was chosen. It exhibited efficient static and dynamic be-
havior, while maintaining flutter stability below the design wind speed of
241 km/hr. The good flutter performance resulted mainly from the smooth
wind flow over the tapered railing and welded fascia plate, and a high
torsional natural frequency.
To illustrate the pK-F method's effectiveness, a flutter analysis was per-
formed on the Luling bridge with the final design cross section. The first
data needed for flutter analysis are natural modes thought to be susceptible
to flutter. The first bending, first sway, and first torsional modes are assumed
in the response analysis. The subspace iteration method with Lanczos start-
ing iteration vectors was used to extract bridge natural modes (Namini 1990).
Only pure modes are used in analysis, in order to illustrate aerodynamic
coupling effects only. Table 1 depicts the natural frequencies with Figs. 5,

,
C i"' /u'"U"ij&-cnjrpU'v \ \ \

28.0 m

FIG. 4. Luling Bridge Deck Section Configuration

TABLE 1. Luling Cable-Stayed Bridge Natural Frequencies and Mode Shapes

Mode designation Natural frequency (Hz) Mode shape
(1) (2) (3)
Bending 0.4282 see Fig. 5
Sway 0.7092 see Fig. 6
Torsion 1.2240 see Fig. 7

1517

J. Struct. Eng., 1992, 118(6): 1509-1526

 FIG. 5. Luling Bridge Fundamental Bending Mode
Downloaded from ascelibrary.org by University of Tasmania Library on 03/22/19. Copyright ASCE. For personal use only; all rights reserved.

FIG. 6. Luling Bridge Fundamental Sway Mode

"%.

FIG. 7. Luling Bridge Fundamental Torsion Mode

6, and 7 displaying the three-dimensional mode shapes for the bending,

sway, and torsional modes, respectively.
Flutter derivatives H\, A*, and A* were determined, while coupling
derivatives H%, H%, and A* were determined to be negligible. Flutter
derivatives associated with drag force, (2), were not determined, but Scanlan
(1987) states quasi-static approximations as

Pt = - | CD (24)

n=~f* : (25)
K daw
P t - -K"2da
— (26)
w
where CD = the static drag coefficient; and aw = the wind angle of inci-
dence. Calculations determined that flutter derivatives, P% and P* were
negligible, and, therefore, were neglected. Figs. 8, 9, 10, and 11 illustrate
the H*, Pt, A*, and A* flutter derivatives at various wind angles of
incidence, respectively.
With the given flutter derivatives and modal information, the pK-F method
described earlier produced the flutter profile, which is tabulated in Table
2. For each analysis, the torsional mode was the critical flutter condition
and is intuitively seen by the positive nature of the A * flutter derivative,
hence, negative damping contribution. Also, the response had no coupling
1518

J. Struct. Eng., 1992, 118(6): 1509-1526

 FLUTTER DERIVATIVE H1 STAR VS REDUCED VELOCITY
FOR LULING CABLE-STAYED BRIDGE

Anqle of Incidence (deq)

1.0- O — 0-4.00
Downloaded from ascelibrary.org by University of Tasmania Library on 03/22/19. Copyright ASCE. For personal use only; all rights reserved.

x 0.0 1 "••'•iSigs^J;

-1.0
0)
Q

% -2.0-^

-3.0
-0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5
Reduced Velocity

FIG. 8. Luling Bridge Flutter Derivative H* versus Reduced Velocity

FLUTTER DERIVATIVE P1 STAR VS REDUCED VELOCITY

FOR LULING CABLE-STAYED BRIDGE

0.2
Anqle of Incidence (deq)
O — O-4.CT0
0.1 @ — ©-1.75
£ A
A
A+0.00
A+2.00
0.0 • D+4.00
• — P+6.00
-0.1 :
H;
0)
Q
-0.2-
* -ft
6;
-0.3-
-8
-0.4
-0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
Reduced Velocity
FIG. 9. Luling Bridge Flutter Derivative P? versus Reduced Velocity

between modal responses, since coupling flutter derivatives were neglected,

and only pure natural modes were used. Fig. 12 illustrates the flutter profile
for the Luling cable-stayed bridge along with the original designer's wind
angle boundary profile. The chosen cross section is acceptable for flutter
stability since the onset of flutter occurs at wind speed larger than the wind
angle boundary profile.
To view the preflutter and postflutter behavior, Figs. 13 and 14 depict
logarithmic decrement-versus-wind velocity and frequency-versus-wind ve-
locity diagrams, respectively at a 6.0° wind angle of incidence. From the
path of the torsional mode shown in Fig. 13, it is obvious that flutter occurs
1519

J. Struct. Eng., 1992, 118(6): 1509-1526

 FLUTTER DERIVATIVE A2 STAR VS REDUCED VELOCITY
FOR LULING CABLE-STAYED BRIDGE

0.35
Anqle of Incidence (deq)
O — 0-4.00
m—• —1.75
Downloaded from ascelibrary.org by University of Tasmania Library on 03/22/19. Copyright ASCE. For personal use only; all rights reserved.

3 0.25- A A+0.00
A A+2.00
3! • n+4.00
S> 0 . 1 5 - • — B+6.00

5> 0.05
o
a> -•-
iSS=8
1-0.05

-0.15
-0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
Reduced Velocity

FIG. 10. Luling Bridge Flutter Derivative A* versus Reduced Velocity

FLUTTER DERIVATIVE A3 STAR VS REDUCED VELOCITY

FOR LULING CABLE-STAYED BRIDGE

0.5
Anqle of Incidence (deq)
O-— 0 - 4 . 0 0
O
0.4- • - -1.75
A- -A+0.00
A- - A + 2.00
< 0.3- •- -D+4.00
-^+6.00
0.2

Q
0.1 n—, $
o
0.0 =|;

-0.1
-0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
Reduced Velocity

FIG. 11. Luling Bridge Flutter Derivative A * versus Reduced Velocity

without any regard to the bending or sway mode types. The instant of flutter
formation, which occurs when the logarithmic decrement becomes zero, is
at a wind velocity of 136 km/hr.
Before the onset of flutter, the assumed mechanical damping of 0.07%,
which was computed experimentally for the rigid model and thus adopted
here, provided enough damping to offset aerodynamic effects. From 136
km/hr and greater, the aerodynamic moment produced more negative damp-
ing than the structure's own mechanical damping. The response frequency,
as shown in Fig. 14, displays a reduction in the torsional flutter response
frequency, which is mathematically explained by the positive values of flutter
1520

J. Struct. Eng., 1992, 118(6): 1509-1526

 TABLE 2. Luling Cable-Stayed Bridge Flutter Profile
Flutter Criterion Based on Flutter Criterion Based on
Current Study Scanlan (1987)
Angle of
incidence Velocity Frequency Velocity Frequency
Downloaded from ascelibrary.org by University of Tasmania Library on 03/22/19. Copyright ASCE. For personal use only; all rights reserved.

(degrees) (km/hr) (Hz) (km/hr) (Hz)

(D (2) (3) (4) (5)
-4.00 391 1.082 394 1.076
-1.75 362 1.053 356 1.044
0.00 345 1.075 346 1.069
2.00 315 1.132 317 1.122
4.00 233 1.157 233 1.144
6.00 136 1.230 135 1.161

FLUTTER VELOCITY PROFILE

FOR LULING CABLE-STAYED BRIDGE
500
- Profile
Boundary
400

300-

o
CD 200
>

,- 100-

- 4 - 2 0 2 4 6
Angle of Incidence (deg)

FIG. 12. Luling Bridge Flutter Profile

derivative A*. The bending and sway modes are stable, since they both
possess negative values of logarithmic decrement, thus convergent oscilla-
tory behavior. The same general preflutter and postflutter tendencies oc-
curred for the other angles of incidence.
With Ld being the deck length, Scanlan (1987) and Scanlan and Jones
(1990) provide a flutter stability criterion for j'th mode as
4 M co,.
H*{K)Ghi + Pt{K)GPi + At(K)Ga (27)
pB4Ld co,
where

1
(28)
l + ^AHK)Ga

and
1521

J. Struct. Eng., 1992, 118(6): 1509-1526

 LOGARITHMIC DECREMENT VS WIND VELOCITY
ANGLE OF INCIDENCE = 6.00 DEG
LULING CABLE-STAYED BRIDGE

0.025
Mode Type
• - Bending
Downloaded from ascelibrary.org by University of Tasmania Library on 03/22/19. Copyright ASCE. For personal use only; all rights reserved.

— Torsion
—•Sway
2 0.012-

0.000

I v
->-\ X- • .
gi - 0 . 0 1 2 -
\
i
I
1
-0.025'
o 50 100 150 200 250 300
Wind Velocity (km/hr)
FIG. 13. Logarithmic Decrement versus Wind Velocity (Damping-Driven Flutter)

FREQUENCY VS WIND VELOCITY

ANGLE OF INCIDENCE = 6.00 DEG
LULING CABLE-STAYED BRIDGE

1.80
Mode Type
• - Bending
— Torsion
1.50 — • Sway

1.20

0.90

0.60--

0.30-
0 50 100 150 200 250 300
Wind Velocity (km/hr)

FIG. 14. Frequency versus Wind Velocity (Damping-Driven Flutter)

Bib,
#/ = U (29)

The variables Gu, Gp,,and Gai = the modal integrals in the heave, sway,
and torsional directions, respectively; while /, and £,- = generalized inertia
and damping ratio, respectively. For the lowest wind velocity that validates
the equality of (27), flutter has been initiated with the corresponding flutter
frequency being computed from (28). Table 2 details the flutter velocity,
1522

J. Struct. Eng., 1992, 118(6): 1509-1526

 TABLE 3. Stiffness-Driven Flutter Derivatives (Scanlan 1981)
Reduced
velocity
UI(NB) A*, A\ At H\ H*2 H%
d) (4) (5)
Downloaded from ascelibrary.org by University of Tasmania Library on 03/22/19. Copyright ASCE. For personal use only; all rights reserved.

(2) (3) (6) (7)

2 0 0 0 -0.67 0 0
4 0 -0.03 0 -1.50 0 -0.05
6 0.75 -0.05 0.50 -2.05 0.70 -1.25
8 0.70 -0.10 1.00 -3.25 2.25 -3.35
10 0.68 -0.14 1.46 -4.25 4.25 -4.00
12 0.70 -0.16 1.69 -5.50 8.90 -5.00

LOGARITHMIC DECREMENT VS WIND VELOCITY

STIFFNESS-DRIVEN FLUTTER

Mode Type
• - Bending
*- 0.200- — Torsion
reme

r
o
Q)
0.100-
I
J
o f
0.000-
hmi

\
•r -0.100- ^s
D
q ** ^
-0.200-

-0.300- , * , 1

200 400 600 800

Wind Velocity (km/hr)

FIG. 15. Logarithmic Decrement versus Wind Velocity (Stiffness-Driven Flutter)

FREQUENCY VS WIND VELOCITY

STIFFNESS-DRIVEN FLUTTER

1.800
Mode Type
• - -Bending
1.500 •—Torsion

1.200

0.900

0.600

0.300
0 200 400 600 800
Wind Velocity (km/hr)

FIG. 16. Frequency versus Wind Velocity (Stiffness-Driven Flutter)

1523

J. Struct. Eng., 1992, 118(6): 1509-1526

 Uf, with the corresponding flutter frequency, oif as determined by (27) and
(28). Good agreement with the pK-F is noted.
To test the pK-F method for stiffness-driven flutter, certain assumptions
were made. In Scanlan's (1981) example demonstrating the root-finding
approach to stiffness-driven flutter, flutter derivatives for a closed trape-
Downloaded from ascelibrary.org by University of Tasmania Library on 03/22/19. Copyright ASCE. For personal use only; all rights reserved.

zoidal box were detailed. In this study, those flutter derivatives, listed in
Table 3, are used to analyze the Luling cable-stayed bridge again. For
assumed mechanical damping of 1% for the fundamental bending and tor-
sion natural modes, the response logarithmic decrement and frequency ver-
sus wind velocity is plotted in Figs. 15 and 16, respectively. Flutter occurs
at 740 km/hr, noted by zero logarithmic decrement. The flutter type is of
the stiffness-driven form since the response frequencies coalesce. The post-
flutter response is not readily available because the pK-F method cannot
distinguish between modes in its solution algorithm if all response frequen-
cies are the same, as is the case here.
Though the stiffness-driven flutter example is fictitious, Scanlan's root-
finding approach was used in a hand computation with a computed flutter
velocity of 783 km/hr being found. This compares well with the pK-F meth-
od's computed velocity.

COMPUTATIONAL CONSIDERATIONS

The pK-F method has been coded into a software system named LEAF
(Library of Engineering Software for Aerodynamic Flutter), currently in
production use by the staff of the Federal Highway Administration's G.S.
Vincent Memorial Aerodynamic Laboratory. The LEAF system executes
on a microcomputer with a maximum of 12 unique natural modes permitted
to participate in the flutter response equations. Also, a maximum of eight
unique representations of flutter derivatives is permitted. The LEAF system,
running on an IBM PS/2 Model 70, typically executes the pK-F method for
a given mode and wind velocity in less than five s, with the Luling cable-
stayed bridge's run times not exceeding three min.

CONCLUSIONS

The developed pK-F method for flutter analysis of cable-suspended bridges

has proved reliable in its solution methodology and execution efficiency.
The method permits bridge designers to possess a single algorithm for damp-
ing- and stiffness-driven flutter. The depicted matrix form of the pK-F
method permits its insertion into finite element systems. The generality of
the method permits more flutter scenarios to be examined, including con-
struction stages and different flutter derivatives for different portions of the
bridge structure.

ACKNOWLEDGMENTS

The writers wish to express their gratitude to all personnel at the Federal
Highway Administration's Turner-Fairbank Research Center for their as-
, sistance and expert criticism. Also, the first writer would like to recognize
the National Highway Institute for the Research Fellowship Grant, which
provided the support funds to complete this research.
1524

J. Struct. Eng., 1992, 118(6): 1509-1526

 APPENDIX I. REFERENCES

Bathe, K. J. (1982). Finite element procedures in engineering analysis. Prentice-Hall,

Inc., Englewood Cliffs, N.J.
Bosch, H. R. (1977). "Aerodynamic investigations of the Luling, Louisiana cable-
stayed bridge." Report No. FHWA-RD-77-161, Federal Highway Administration,
Downloaded from ascelibrary.org by University of Tasmania Library on 03/22/19. Copyright ASCE. For personal use only; all rights reserved.

Offices of Res. and Development, Washington, D.C.

Hassig, H. L. (1971). "An approximate true damping solution of the flutter equation
by determinant iteration." /. Aircr., (8)11, 876-880.
Namini, A. H. (1990). "Microcomputer-based extraction of cable-stayed bridge nat-
ural modes." Proc, Forth Rail Bridge Centenary Conf., 1, 499-510.
Scanlan, R. H. (1975). "Recent methods in the application of test results to the wind
design of long, suspended-span bridges." Report No. FHWA-RD-75-115, Federal
Highway Administration, Offices of Res. and Development, Washington, D.C.
Scanlan, R. H. (1978). "The action of flexible bridges under wind, Part I." /. Sound
Vib., 60(2), 187-199.
Scanlan, R. H. (1981). "State-of-the-art methods for calculating flutter, vortex-in-
duced, and buffeting response of bridge structures." Report No. FHWA/RD-80/
50, Federal Highway Administration Offices of Res. and Development, Washing-
ton, D.C.
Scanlan, R. H. (1987). "Interpreting aeroelastic models of cable-stayed bridges." /.
Engrg. Mech., ASCE, 113(4), 555-575.
Scanlan, R. H., and Jones, N. P. (1990). "Aeroelastic analysis of cable-stayed bridges."
J. Struct. Engrg., ASCE, 116(2), 279-297.
Scanlan, R. H., and Tomko, J. J. (1971). "Airfoil and bridge deck flutter deriva-
tives." J. Engrg. Mech. Div., ASCE, 97(6), 1717-1737.
Theodorsen, T. (1934). "General theory of aerodynamic instability and the mech-
anism of flutter." NACA Report No. 496, U.S. Government Printing Office, Wash-
ington, D.C.

APPENDIX II. NOTATION

The following symbols are used in this paper:

A = aerodynamic joint load vector;

Ad = aerodynamic damping component matrix;
Af = generalized aerodynamic damping component matrix;
A* = torsional flutter derivative, i = 1, 2, 3;
K = aerodynamic stiffness component matrix;
Af = generalized aerodynamic stiffness component matrix;
B = deck width;
CD = static drag coefficient;
D = damping matrix;
D" = generalized damping matrix;
DP = aerodynamic drag force;
F = flutter determinant;
Gh = heave modal integral;
Gp = sway modal integral;
Ga = torsional modal integral;
H* = heave flutter derivative, i = 1, 2, 3;
h = heave deformation;
I = generalized inertia;
i = unit value along imaginary axis of complex plane;
K = stiffness matrix;
K* = generalized stiffness matrix;
1525

J. Struct. Eng., 1992, 118(6): 1509-1526

 K = reduced frequency;
L = member length;
Ld = deck length;
Lh = aerodynamic lift force;
M = mass matrix;
Downloaded from ascelibrary.org by University of Tasmania Library on 03/22/19. Copyright ASCE. For personal use only; all rights reserved.

Mg = generalized mass matrix;

Ma = aerodynamic moment;
m = number of modes permitted to participate in response;
N = response frequency in Hertz;
n = number of distinct degrees of freedom;
P* = lateral flutter derivative, i = 1, 2, 3;
p = sway deformation;
q = joint deformation vector;
R = response amplitude vector;
S = Strouhal Number;
t = time;
U = wind velocity;
C/cur = current wind velocity;
Uf = wind velocity at flutter;
Uhigh — ending or highest wind velocity of range;
Uinc = increment of wind velocity;
£Aow = starting or lowest wind velocity of range;
a = torsional deformation;
aw = wind angle of incidence;
8 = response logarithmic decrement;
I = damping ratio;
u, = damped harmonic frequency;
£, = generalized coordinate vector;
p = air mass density;
<& = natural mode shape matrix;
co = response circular frequency; and
oif = response circular frequency at flutter.

1526

J. Struct. Eng., 1992, 118(6): 1509-1526