- Email: [email protected]
Dynamic instability of suspension bridges
Cocnpulcrs & Strvcfwes Vol. 41, No. 6, pp. 1321-1328, 1991 Printed in Great Britain.
DYNAMIC
004s7949/91 s3.lxJ + 0.00 Q 1991 Civil-Camp Ltd and Pergamon Prw plc
INSTABILITY
OF SUSPENSION
BRIDGES
T. J. A. AGAR Department of Civil Engineering, Rankine. Building, University of Glasgow, Glasgow G12 8LT, U.K. Abstract-Suspension bridges are long, slender flexible structures which have the potential to be susceptible to a variety of types of wind-induced instabilities, the most serious of which are divergence (due to stationary wind forces) and flutter (due to aerodynamic forces). Flutter occurs at certain wind speeds where aerodynamic forces acting on the deck feed energy into an oscillating structure, so increasing the vibration amplitudes. If this situation is approached the basic safety of the bridge is threatened. This paper describes a computational method for predicting flutter speed based on a modal technique. A selection of the lowest vertical and torsional natural mode shapes is included with the aerodynamic forces in an interaction analysis which yields an unsymmetric matrix eigenvalue problem, the roots of which indicate if flutter is possible. The paper addresses the question of how the degree of refinement of the basic structural model and the number of natural modes included in the interaction analysis affect the flutter speed predictions.
INTRODUCI’ION
In the early part of this century, developments in the aesthetics of bridge building and improvements in materials led to the construction of progressively longer, structurally more efficient and slender bridges. It was only after the Tacoma Narrows suspension bridge collapsed shortly after its completion in 1940 that the potentially unstable behaviour of this type of bridge under wind action began to be investigated. It subsequently became appreciated that such bridges could be liable to violent aerodynamic oscillations that could not be diagnosed from static analyses taking account of steady wind forces, even those for maximum design wind speeds. One of the essential requirements of modern suspension bridge design is to avoid significant levels of wind excited oscillations of which there are two main types: limited amplitude (non-divergent) oscillations produced by vortex shedding; and divergent oscillations produced by both galloping and classical flutter types of instability. The former class of behaviour may, in limit state terminology, be considered primarily as a ‘serviceability’ problem responsible mainly for excessive levels of vibration and having a potential for serious fatigue damage in the long term. In contrast, the latter class, in particular flutter, may be considered to be an ‘ultimate’ condition. Flutter is a self-excited oscillatory instability of a body suspended in an air stream. It involves the interaction of aerodynamic, inertial and elastic structural forces such that, at certain critical wind speeds, the aerodynamic forces act to feed energy into the oscillating structures and increase the magnitude of vibration, sometimes to catastrophic levels. The critical flutter speed of a bridge depends on its vertical and torsional natural vibration characteristics and also on the cross-sectional shape of its deck (since this us
4,/6-N
affects the aerodynamic forces that act during pitch and heave oscillations). Current U.K. Code guidance on stability is available for structures with individual spans not exceeding 200 m[l]. However, for spans greater than this, advice is that the stability should be verified by wind-tunnel tests, which are both time-consuming and expensive. As it was useful at the design stage to be able to examine the effect of different structural cor@uration of a major suspension bridge being designed by Mott, Hay & Anderson, London for a Far East client in relation to, among other things, the flutter problem, a computational package ANSUSP was developed by the author [2] to analyse suspension bridge three-dimensional dynamic behaviour. Initially, the program predicted flutter behaviour in a module using a time integration approach and which included the effects of geometrical nonlinearity. The aerodynamic deck excitation forces used were the analytic Theodorsen functions appropriate to an oscillating flat plate. This paper describes some features of a modal technique for predicting flutter (developed at Glasgow University) which has been incorporated into the ANSUSP package.
DYNAMIC ANALYSIS
The equation in terms of a nodal displacement vector {V} for a structural system can be expressed in conventional matrix notation as
[Ml{lf}+[Cl{il}+[KlIU}={P}
(1)
where the matrices [Ml, [Cl, [K] and {P} have the usual meaning. For a linear elastic structural system, a commonly used technique for solving eqn (1) is that of modal synthesis. This has the computational
1321
1322
T. J. A. AGAR
advantage that the response of a structure discretized as a large order n degrees of freedom configuration can be effectively reduced to one with m degrees of freedom (m Q n), where m is the number of modes included. In general, there is a nonlinear relationship between forces and displacements due to significant changes in the geometry of a suspension bridge as it deflects. Hence the stiffness matrix [K] is displacement dependent and the modal synthesis technique is not strictly appropriate. An alternative for solving equation (1) is numerically to integrate the equation of motion directly to give a time-history response, and various methods have been available for perform this [3-61. One of these methods, the Newmark fl method, was initially used in modules of ANSUSP to predict general structural response including fiutter. Details of this approach have heen given elsewhere [7].
Classical flutter is an aeroelastic phenomenon in which two degrees of freedom of a structure, a rotation and a translation, couple together in a flowdriven unstable oscillation. Coupling of the two degrees of freedom, indispensable for flat-plate thinaerofoil flutter has come to be the identifying sign for classical flutter. Single degree of freedom flutter can be associated with bluff, unfaired bodies where the Sow is strongly separated. Some of the older suspension bridge deck sections could give rise to single degree of freedom flutter, but modern designs tend to make the aerodynamics of the deck closely resemble that of a flat plate. Consequently the aerodynamic lift L, and moment Ma acting on an oscillating deck cross-section (Fig. 1) have been assumed to be reasonably approximated by the theoretically derived Theodorsen functions [S] which depend on deck vertical displacement h and twist ~1, their time derivatives and the reduced frequency of vibration of the system {orb/l/). It has not been possible to develop expressions for aerodynamic coefficients associated with bluff bodies
b
B _I=
b
---=-$_ _Ic
2
---
H h
k-f
h
Fig. 1. Aerodynamic forces on oscillating flat plate.
from basic fluid-flow principles. However, Scanlan and Tomko have shown experimentally that for small displacements the lift and moment may be treated as linear functions of h, a and their first two derivatives [9]. The coeEcient m~tipliers on these derivatives have also been determined experimentally in Japan and France [lo, 11, 121. structural model The program ANSUSP idealizes a suspension bridge as a three-dimensional framework as described by Agar[7& It is a two-cable idealization and comprises the types of structural elements shown in Fig. 2. Since the flutter phenomonen is normally considered to involve motion in the xy plane, unaccompanied by a transverse motion, vertical and torsional bridge vibrations are considered to consist essentially of (i) x, y, 8, displacements of cable nodes; (ii) y, 8,, d, displacements of deck nodes; (iii) x, 6, displa~men~ of tower top nodes. The global stiffness matrix is assembled from contributions from the deck, hanger, tower and cable element stiffness matrices, the last of these including a geometric component that accounts for the gravity stiffness of the structure [13]. The inertial character-
Cableelements
Rigid arm links
Fig. 2. Suspension bridge idealization in 3-D space.
l3ynamic ins~bi~t~ of
istics of the structure are modelled by lumping of the member element mass at the nodes. Having formed the stiffness and mass matrices, the resulting free vibration eigenvalue problem is solved by a sim~tan~us iteration algorithm for a specified number of lower vibration modes [x] [14]. These natural frequency and mode shape characteristics are permanently stored in files for use in subsequent modal flutter analyses. Modal $utter analysis For a suspension bridge subject to a transverse wind flow the right-hand side of eqn (1) becomes V,,}
=
Ml{UJ + PIW
where
1m-J%
w*1= (P I- WlT‘f1[xl) and [Xx)is a matrix containing m natural modes and [II] is a diagonal matrix cont~ning the ~~spon~ng circular frequencies squared. Looking for a solution to eqn (3) of the form {4) = {@}e” leads to the standard matrix eigenvalue form of order 2m
where ]q=[-F*
et*]
and
{xI=~:}.
Solution of eqn (4) yields a set of t and (01. Omitting vector and matrix brackets hereafter, for a complex conjugate pair of eigenvalues 1 = p + io and 3, = p - io and corresponding eigenvalues Cp= p + iq and Q, = p - iq the characteristic motion may be expressed as r$ =epf[(p +q)sinwt
+(q -p)coswt].
(5)
To ensure dynamic stability in a system, the real parts p of all the eigenvalues should be negative. The onset of flutter instability occurs at the lowest wind speed for which an eigenvalue is complex and has a zero real part. Divergence instability is indicated by a real ~genvalue becoming zero. Having the natural modes available, the calculation of matrix [a of order 2m in eqn (4) is straightforward, only the two upper sub-matrices C* and K* requiring significant computation. In calculating C*, normally C, the structural damping
1323
matrix, is conservatively taken as being zero. The eigensolution of eqn (4) is performed as follows: (a) The real unsymmetric matrix K is reduced to a matrix ii of upper Hessenberg form by similarity tr~sfo~ations. The ~~nvalues of 2 are the same as those of 2. (b) The eigenvalues Izi of a (and hence A, are evaluated using the Double QR method, which is a computationally more efficient extension of Francis’s QR method [ 151. (c) If required, any eigenvector {zj> resending to eigenvalue li is then computed by a back substitution process.
(2)
where the coefficients of [A] and [B] arise from the Theodorsen equations or expe~mentally measured flutter derivatives described above. Applying the modal synthesis technique leads to an m th order system in terms of generalized coordinates
]C*l= Kl’(tC] -
suspension bridges
SOLUTION CHARACTERISTICS
Since the method predicts eigenvalues and vectors defining the response for the given wind speeds, investigation of flutter stability involves a sweep through a range of increasing windspeeds until instability is indicated. Furthermore, the solution of the eigenproblem as in (a)-(c) above at any given wind speed must be iterative in itself. This is because of the aerodynamic terms are functions of response frequency, and re-analysis has to be performed until o,~~, the circular frequency used to evaluate the aerodynamic terms, agrees with the response of interest. Initially, with very low wind speeds, the response frequencies are almost equal to the natural frequencies and the real parts of the eigenvalues are negative, indicating decaying motions. If the response frequencies are plotted on an Argand diagram they all he just to the left of the imaginary axis. When analyses are carried out at successively higher wind speeds, all eigenvalues start to migrate along paths further away from the ima~nary axis, co~s~nding to increasing damping, while oscillation frequencies may increase or decrease from the natural frequency values. Eventually at some higher wind speed, one of these complex eigenvalue paths turn and heads back to intersect the axis. This corresponds to an oscillation with zero ~mping and the onset of flutter instability. Figure 3 shows an example of the response frequencies plotted on an Argand diagram where the flutter speed V, = 76.9 m/s. The characteristics and problems associated with the modal flutter solution have been described in detail by the author elsewhere [16]. NUMERICAL ACCURACY OF FLlJTl’JIR SPEED PREDICTIONS Efect
of basic idealisation
Figure 4 shows a typical structural idealization of the Sevem Bridge where the left, main and right side span decking is represented by 5, 16 and 5 beam elements, respectively. With the properties as shown in Fig. 4, natural frequency analyses in ANSUSP for various idealization refinements yield the natural frequency predictions as indicated in Fig. 5. Table 1
T. J. A. AGAR
1324
Mode 12
0-0-o
m-O-*-
laa++~-a--e_.
!
Mode 10 3.0
12% (““i 2.5
HI-o-e-o-o-o
2.0
1.0
o-e-.-~-
C.
I
I
-0.3
-0.2
Mode 2
c
0
-0.1 IJ
Fig. 3. Argand plot for 5-16-5 idealization of Sevem Bridge with first seven symmetric modes from Fig. 5. The (0) represent plotted points for wind speeds V (m/s): 10, 20, 30, 40, 50, 60, 65, 70, 72.5, 75, 77.5.
shows the CPU times required to compute the lowest 13 and 30 frequencies and mode shapes for the range of idealizations considered. All computation has been performed on a Dee uVAX II computer. Efect of natural modes inclua’ed If a modal eigenvalue analysis is performed, for the 5-16-5 model first including only modes 2 and 8, i.e. the flutter response is considered to be a combination of the fundamental symmetric flexural and torsional modes, then the Argand diagram of Fig. 6 is produced, where the vibration frequency at flutter (1 = 0) follows a path starting at the fundamental torsional frequency while the path of the bending frequency is directly away from the o axis. Including only these two modes, a flutter windspeed V, of -71.0 m/s and flutter frequency W, of 1.545 rad/s is predicted.
Table 1. CPU time for natural frequency analysis for various structural idealizations of the Scvem Bridge
Idealization
Number of degrees of freedom
2-6-2 3-8-3 5-16-5 8-26-8 10-32-10 12-36-12 17-54-17
71 107 215 359 449 561 773
CPU time (set) 13 modes 30 modes 36 55 112 183 240 276 414
140 184 344 627 763 887 1321
An analysis including the additional mode 6 leads to Fig. 7 where flutter is predicted at V, 76.9 m/s and at an oscillation frequency o,= 1.49 rad/s. Now consider an analysis including modes 2, 3, 6, 7, 8, 10, 12, i.e. the lowest 7 symmetric modes from Fig. 4. For the (symmetric) Severn Bridge, the lowest flutter speed corresponds to a symmetric displacement pattern with antisymmetric natural modes not contributing to the response. The Argand diagram is shown in Fig. 3 where the flutter prediction is again at V, = 76.9 m/s and or = 1.495 rad/s. Table 2 shows the flutter speed predictions and modal involvement in the flutter response obtained by including different combinations of modes in an analysis. The conclusion to be drawn from a range of the above types of analysis is that, in this case, the flutter response consists of a combination of modes 2,6 and 8, i.e. the first and third symmetric flexural and fundamental torsional modes. Figure 8 summarizes the flutter speed predictions obtained from the range of idealizations for various mode shapes included in the modal flutter analysis. It can be seen that the flutter speed predictions are (apart from the relatively coarse 2-6-2 and 3-8-3 idealizations) very insensitive to the degree of refinement of the starting structural model. However, two distinct values of flutter speed of N 71 and N 77 m/s are predicted depending on which mode shapes are included in the modal flutter analysis. Figure 9 shows the lowest three symmetric flexural and the lowest symmetric torsional mode shapes. Distance between cables = 22.87
g1’5rI
r[
1 10.2m
5 elements
I-
8 elements 488m
I:
305 m
-I
Structural properties
Cables
: CSA
Hangers
:
Towers
:
Mass CSA Mass Itlong.) W&y)
= = = = = = =
0.158m’ 1339 kg/m 0.0146 m2 124 ka/m 2.13& 8.52 m4 207 x lo6 kN/m2
Deck
: Itvert) J (torsion) Mass Mass moment of inertia Overall width
= = = =
1.079m4 3.995m4 10350 kg/m 696ooo kg-m*/m
=
31.85 m
Fig. 4. Sevem Bridge 5-165 idealization and properties.
Dynamic instability of suspension bridges
1325
Table 2. Flutter speed predictions and modal involvement in response for analyses with different combinations of symmetric mode shapes for 5-16-5 idealization Mode shape
Modal iavolvcmcnt as % flutter response
Flexural 2 3 6 7 10 12
100.0 -
85.4 14.6 -
65.7 34.3 -
57.3 9.8 30.7 2.1 UO 0.1
59.3 9.1 31.6 -
Torsional 8 18
100.0 -
100.0 -
100.0 -
100.0 -
98.6 1.4
q(rad/s) Vhmls) Note -,
1.545 71.0
1.551 70.9
1.495 76.9
1.495 76.9
1.495 76.9
indicates mode not included in analysis.
DISCUSSION
and flexural (mode 2) modes with V,= 77 m/s obtained by including at least modes 2, 6 and 8, there may appear to be some inconsistency. However, on examination the inference from the former analysis is that if the bridge were artificially constrained to allow
Initially, on comparing the flutter prediction of V,= 71 m/s obtained by including only the similarly shaped fundamental symmetric torsional (mode 8) Natural frequency o (radk)
A
_---
Mode
Designation+
-029
3-S-T(C)
*___*_-__*_____--
,*fl-
/
/
1’
8-
+ Mode designation: e.g. 3-S-F(C,S) implies 3rd symmetricflexural mode involving displacements in the centreand side spans.
1’ f
I
I 7-
__~_--.---_-~--- -_-_-
023 2-A-T(C)
,/’ /
/ 6-
.---
./ /
.’
/ _.____--.----.----
.C
.--_-__-_-.21
l-S-T(S) l-A-T(S)
.----_
2-S-T(C,S)
----_.18
l____~_____-~--__~----_.--__-- -413 l-A-T(C) l______._-_.---_.---
&i w l/.--,*-/
- - - -412
/‘-
/
--
,.______ /-----
.---._.--_-_-_-_ _.---.----.--------.g
---_.
4-S-F(C) 11 P-A-F(S) 10 l-S-F(S) *A-F(C)
.@ _*-----He*-----
$ -__ ---
$====$==I==-_==:;
1;-J(y
1 .~.----__-__.__-.-_-_.
.,c
-a_.
------.---.----.--------43
--------a
4
l-A-F(S) 2_S_F(CS)
n
Fig. 5. Natural frequencies of Swem Bridge for various idealizations.
1326
T. J. A. AGAR
t
!2.5
V, = 76.9 mls q = t .49 radls (Wtrial following
jO.5
I
-0.3
CL
wind speeds V (m/s): IO, 20, 30, 40, SO,60, 65, 70, 72.5.
flexural motion
proportional to the fundamental symmetric mode only, flutter could appear as predicted at 71 m/s, However, in reality the response induced by the aerodynamic forces will be the combination of all mode-shapes that the aerodynamics drive. The flexural response at flutter shown in Fig. 10 is consistent with that of Bleich, who reported the effect of incomplete matching between flexural and torsional modes as increasing flutter speed as the initially half-sine wave fundamental flexural mode develops a reverse curvature 1171.
80-
Predicted Flutter Speed V, (m/s)
I
-0.1
0
W
CL
.
Fig. 6. Argand pIot for 5-16-5 idealization of Sevem Bridge with modes 2 and 8. The (0) represent plotted points for
t
I
-0.2
Fig. 7. Argand plot for 5-46-5 idealization of Sevem Bridge with modes 2,6 and 8. The (e) represenl plotted points for wind speeds V(m/s): 10,20,30,40,50,60, 65,70,72.5,75.0, 77.5.
The conclusion is that in order to obtain the correct flutter prediction then modes 2, 6 and 8 at least, must be included. However, as can be seen from Fig, 8 there is no significant change in the flutter speed prediction if additional mode shapes are included. Comparison of the flutter wind speed V, = 76.9 m/s obtained above with predictions by other methods are shown in Table 3. Use of the numerical integration time-history module in ANSUSP (which uses exactly the same formulation for the aerodynamic Modes included + 2,8 (B2,3,8 A 2,6,8 o 2,3,6,8 v 2,3,6,8,18 02,3,6,7,8,10,12,18
Fig. 8. Flutter speed dependence on basic structural idealization and natural modes included.
Dynamic ~s~bitity of suspension bridges Table 3. Comparison of flutter windspeed predictions for the Sevem Bridge by various methods Method
Reference
Multi degree-of-freedom, modal Multi degree-of-freedom, time&story Two degrees-of-freedom, time-history (by author) Two degrees-of-freedom, modal (by author) Two degrees-of-freedom, experimental and calculation (Smith-NPL) Selberg semi-empirical
This method
Flutter speed (m/s) 76.9 77.5
Agar VI -
65.6
-
65.6
Smith [18]
65.7
Selberg [19]
68
force coefficients) produces a value in the region of 77.5 m/s with a flexural response showing a reverse cu~atu~ very similar to that shown in Fig. 10. The other predictions in Table 3 are based on modelling the behaviour with 2 degrees of freedom only and give consistently lower results. The two predictions performed by the author use the same Theodorsen aerodynamic forces as the multi-degree of freedom modal method described here and the time integration procedure. Hence there is a substantiation for the difference in flutter predictions obtained from multi-degree of freedom and two-degree of freedom models.
Mode2
1327
In practice, the static divergence instability is predicted to occur at 74.7 m/s, so that the theoretical increase in flutter speed (from approximately 66 m/s for two-degree of freedom to approximately 77 m/s for multi-degree of freedom methods) is to some extent academic in this case. From the computational efficiency viewpoint, comparison of the modal flutter method described here and the time integration approach shows the former to be much superior. With the modal approach the majority of the computation is involved in the initial natural frequency analysis (e.g. see Table 1: 334 CPU seconds to obtain 30 modes of the 5-16-5 model), while the computation involved in the subsequent modal eigenvalue analysis is modest. It can be shown that an approximate semi-empirical expression for the number of CPU seconds needed for convergence of the modal flutter solution for each windspeed is given by 0.03(171)2~5 where fi is the number of natural modes included in the eigenvalue flutter analysis. By comparison, the multi-degree-of-freedom timehistory integration method with the same S-16-5 structural idealization required N 15OOCPUs per wind speed considered. Con~uently, use of the modal flutter analysis method described in the paper represents a saving in computation by a factor of between 20 and 50 over the time integration approach for comparable accuracy.
(I-S-F)
E
1
~odr3(2-S-F
Moda 6 (3-S-F)
c _---
---_ -I, --.
Mod. 8 (I-S-T
1
-z
-----b--e
Fig. 9. Natural mode shapes 2, 3, 6 and 8 for Sevem Bridge analysis. Tower displacements shown at X5 magnification.
T. J. A. AGAR
1328
0.125
0.25
0.375
0.5
Fig. 10. Deck vertical displacements at flutter from mode shapes 2, 3 and 6. Curves represent displacements at times indicated:
time
Curve
1 2 3 4 T =
CUN’Z
time
t, t + T
5
t + (4T/8)
t + O-/8) t + (2T/8) t + (3T/8)
6 7 8
t + (ST/E) t + (6T/8) 6 + (7T/8)
Period of oscillation = 2x/w,.
CONCLUSIONS A method of computing the flutter speeds of suspension bridges has been presented. The method uses a multi-degree-of-freedom modal technique and includes the aerodynamic force interaction in the governing equations of motion to produce an unsymmetric eigenvalue problem. The characteristic roots of the matrix involved are indicative of the nature of possible motions that can exist at a given wind speed. A study carried out to consider effects of the degree of refinement of the basic structural model and the number of modes included in the eigenvalue flutter analysis has shown that flutter speed predictions are in-sensitive to the basic suspension bridge model, but that it may be important to include more than the fundamental mode shapes on which the simpler 2-degree of freedom methods are based. The modal method described is demonstrated to be computationally much more efficient than an alternative time-history integration method, but yields comparable flutter speed prediction accuracy. REFERENCES 1. Bridge Aerodynamics, Proposed British Design Rules. Thomas Telford, London (1981). 2. T. J. A. Agar, Analysis of suspension bridges-programme ANSUSP user guide. Mott Hay & Anderson Computer Department (1980). 3. E. L. Wilson and R. W. Clough, Dynamic response by a step-by-step analysis. Proc. Symp. on the Use of Computers in Civil Engineering, Lisbon. (1982). 4. K. J. Bathe and E. L. Wilson, Stability and accuracy analysis of direct integration methods. Earthquake Eng. Struct. Dynumic 1, 283-291 (1973). 5. H. M. Hilber, T. J. R. Hughes and R. L. Taylor, Improved numerical dissipation for time integration algorithms in structural dynamics. Earthquake Eng. Struct. Dynamic 5, 283-292 (1977).
6. N. M. Newmark, A method of computation for structural dynamics. Proc. ASCE J. 66, 67-94 (1959). 7. T. J. A. Agar, The analysis of aerodynamics flutter of suspension bridges, Comput. Strut. 30, 593-600 (1988). 8. T. Theodorsen, General theory of aerodynamic instability and the mechanism of flutter, NACA Technical Report, No. 496 (1935). 9. R. H. Scanlan and J. J. Tomko, Airfoil and bridge deck flutter derivatives. Proc. ASCE J. 97, 1717-1737 (1971). 10. T. Okubo and N. Narita, A comparative study of aerodynamic forces acting on cable-stayed bridge girders. Proc. 2nd US-Japan Seminar on Wind Effects on Structures, Kyoto, University of Tokyo Press, pp. 271-283
1
(1976). T. Okubo and K. Yokoyama, Some approaches for improving wind stability of cable-stayed bridges. Proc. 4th Int. Conf. on Wind Efects on Buildings and Structures, London, 1985. Cambridge University Press,
1
pp. 241-249, 1986. H. Loiseau and E. Szechenyi, Etude du comportement aeroelastique du tablier d’un pont a haubans, T.P.
1975-75, Office National d’Etudes et de Reeherches Aerospatiales, Chatillon, France. 13. J. S. Przemieniecki, Theory oj Matrix Structural Analysis. McGraw-Hill, New York (1968). 14. R. B. Corr and A. Jennings, A simultaneous iteration algorithm for symmetric eigenvalue problems. Int. J. Num. Meth Engng 10, 647-663 (1976). 15. J. G. F. Francis, The QR transformation, parts I and II. Computer J. 4, 265-271, 332-345 (1961). 16. T. J. A. Agar, Aerodynamic flutter analysis of suspension bridges by a modal technique. Eng. Structures 11, 75582 (1989). 17. F. Bleich, The flutter theory. The Mathematical Theory of Vibration in Suspension Bridges (Edited by Bleich, F., McCullough, C. B., Rosecrans, R. and Vincent, G. S.), U.S. Government Printing Office, Washing, DC, Ch. 7, pp. 241-281 (1950). 18. I. P. Smith, The aeroelastic stability of the Sevem suspension bridge. NPL Aero Report 1105 (1964). 19. A. Selberg, Oscillation and aerodynamic stability of suspension bridges. Acta Polytech. Stand., Civil Eng. No. 13 (1961).
DYNAMIC
004s7949/91 s3.lxJ + 0.00 Q 1991 Civil-Camp Ltd and Pergamon Prw plc
INSTABILITY
OF SUSPENSION
BRIDGES
T. J. A. AGAR Department of Civil Engineering, Rankine. Building, University of Glasgow, Glasgow G12 8LT, U.K. Abstract-Suspension bridges are long, slender flexible structures which have the potential to be susceptible to a variety of types of wind-induced instabilities, the most serious of which are divergence (due to stationary wind forces) and flutter (due to aerodynamic forces). Flutter occurs at certain wind speeds where aerodynamic forces acting on the deck feed energy into an oscillating structure, so increasing the vibration amplitudes. If this situation is approached the basic safety of the bridge is threatened. This paper describes a computational method for predicting flutter speed based on a modal technique. A selection of the lowest vertical and torsional natural mode shapes is included with the aerodynamic forces in an interaction analysis which yields an unsymmetric matrix eigenvalue problem, the roots of which indicate if flutter is possible. The paper addresses the question of how the degree of refinement of the basic structural model and the number of natural modes included in the interaction analysis affect the flutter speed predictions.
INTRODUCI’ION
In the early part of this century, developments in the aesthetics of bridge building and improvements in materials led to the construction of progressively longer, structurally more efficient and slender bridges. It was only after the Tacoma Narrows suspension bridge collapsed shortly after its completion in 1940 that the potentially unstable behaviour of this type of bridge under wind action began to be investigated. It subsequently became appreciated that such bridges could be liable to violent aerodynamic oscillations that could not be diagnosed from static analyses taking account of steady wind forces, even those for maximum design wind speeds. One of the essential requirements of modern suspension bridge design is to avoid significant levels of wind excited oscillations of which there are two main types: limited amplitude (non-divergent) oscillations produced by vortex shedding; and divergent oscillations produced by both galloping and classical flutter types of instability. The former class of behaviour may, in limit state terminology, be considered primarily as a ‘serviceability’ problem responsible mainly for excessive levels of vibration and having a potential for serious fatigue damage in the long term. In contrast, the latter class, in particular flutter, may be considered to be an ‘ultimate’ condition. Flutter is a self-excited oscillatory instability of a body suspended in an air stream. It involves the interaction of aerodynamic, inertial and elastic structural forces such that, at certain critical wind speeds, the aerodynamic forces act to feed energy into the oscillating structures and increase the magnitude of vibration, sometimes to catastrophic levels. The critical flutter speed of a bridge depends on its vertical and torsional natural vibration characteristics and also on the cross-sectional shape of its deck (since this us
4,/6-N
affects the aerodynamic forces that act during pitch and heave oscillations). Current U.K. Code guidance on stability is available for structures with individual spans not exceeding 200 m[l]. However, for spans greater than this, advice is that the stability should be verified by wind-tunnel tests, which are both time-consuming and expensive. As it was useful at the design stage to be able to examine the effect of different structural cor@uration of a major suspension bridge being designed by Mott, Hay & Anderson, London for a Far East client in relation to, among other things, the flutter problem, a computational package ANSUSP was developed by the author [2] to analyse suspension bridge three-dimensional dynamic behaviour. Initially, the program predicted flutter behaviour in a module using a time integration approach and which included the effects of geometrical nonlinearity. The aerodynamic deck excitation forces used were the analytic Theodorsen functions appropriate to an oscillating flat plate. This paper describes some features of a modal technique for predicting flutter (developed at Glasgow University) which has been incorporated into the ANSUSP package.
DYNAMIC ANALYSIS
The equation in terms of a nodal displacement vector {V} for a structural system can be expressed in conventional matrix notation as
[Ml{lf}+[Cl{il}+[KlIU}={P}
(1)
where the matrices [Ml, [Cl, [K] and {P} have the usual meaning. For a linear elastic structural system, a commonly used technique for solving eqn (1) is that of modal synthesis. This has the computational
1321
1322
T. J. A. AGAR
advantage that the response of a structure discretized as a large order n degrees of freedom configuration can be effectively reduced to one with m degrees of freedom (m Q n), where m is the number of modes included. In general, there is a nonlinear relationship between forces and displacements due to significant changes in the geometry of a suspension bridge as it deflects. Hence the stiffness matrix [K] is displacement dependent and the modal synthesis technique is not strictly appropriate. An alternative for solving equation (1) is numerically to integrate the equation of motion directly to give a time-history response, and various methods have been available for perform this [3-61. One of these methods, the Newmark fl method, was initially used in modules of ANSUSP to predict general structural response including fiutter. Details of this approach have heen given elsewhere [7].
Classical flutter is an aeroelastic phenomenon in which two degrees of freedom of a structure, a rotation and a translation, couple together in a flowdriven unstable oscillation. Coupling of the two degrees of freedom, indispensable for flat-plate thinaerofoil flutter has come to be the identifying sign for classical flutter. Single degree of freedom flutter can be associated with bluff, unfaired bodies where the Sow is strongly separated. Some of the older suspension bridge deck sections could give rise to single degree of freedom flutter, but modern designs tend to make the aerodynamics of the deck closely resemble that of a flat plate. Consequently the aerodynamic lift L, and moment Ma acting on an oscillating deck cross-section (Fig. 1) have been assumed to be reasonably approximated by the theoretically derived Theodorsen functions [S] which depend on deck vertical displacement h and twist ~1, their time derivatives and the reduced frequency of vibration of the system {orb/l/). It has not been possible to develop expressions for aerodynamic coefficients associated with bluff bodies
b
B _I=
b
---=-$_ _Ic
2
---
H h
k-f
h
Fig. 1. Aerodynamic forces on oscillating flat plate.
from basic fluid-flow principles. However, Scanlan and Tomko have shown experimentally that for small displacements the lift and moment may be treated as linear functions of h, a and their first two derivatives [9]. The coeEcient m~tipliers on these derivatives have also been determined experimentally in Japan and France [lo, 11, 121. structural model The program ANSUSP idealizes a suspension bridge as a three-dimensional framework as described by Agar[7& It is a two-cable idealization and comprises the types of structural elements shown in Fig. 2. Since the flutter phenomonen is normally considered to involve motion in the xy plane, unaccompanied by a transverse motion, vertical and torsional bridge vibrations are considered to consist essentially of (i) x, y, 8, displacements of cable nodes; (ii) y, 8,, d, displacements of deck nodes; (iii) x, 6, displa~men~ of tower top nodes. The global stiffness matrix is assembled from contributions from the deck, hanger, tower and cable element stiffness matrices, the last of these including a geometric component that accounts for the gravity stiffness of the structure [13]. The inertial character-
Cableelements
Rigid arm links
Fig. 2. Suspension bridge idealization in 3-D space.
l3ynamic ins~bi~t~ of
istics of the structure are modelled by lumping of the member element mass at the nodes. Having formed the stiffness and mass matrices, the resulting free vibration eigenvalue problem is solved by a sim~tan~us iteration algorithm for a specified number of lower vibration modes [x] [14]. These natural frequency and mode shape characteristics are permanently stored in files for use in subsequent modal flutter analyses. Modal $utter analysis For a suspension bridge subject to a transverse wind flow the right-hand side of eqn (1) becomes V,,}
=
Ml{UJ + PIW
where
1m-J%
w*1= (P I- WlT‘f1[xl) and [Xx)is a matrix containing m natural modes and [II] is a diagonal matrix cont~ning the ~~spon~ng circular frequencies squared. Looking for a solution to eqn (3) of the form {4) = {@}e” leads to the standard matrix eigenvalue form of order 2m
where ]q=[-F*
et*]
and
{xI=~:}.
Solution of eqn (4) yields a set of t and (01. Omitting vector and matrix brackets hereafter, for a complex conjugate pair of eigenvalues 1 = p + io and 3, = p - io and corresponding eigenvalues Cp= p + iq and Q, = p - iq the characteristic motion may be expressed as r$ =epf[(p +q)sinwt
+(q -p)coswt].
(5)
To ensure dynamic stability in a system, the real parts p of all the eigenvalues should be negative. The onset of flutter instability occurs at the lowest wind speed for which an eigenvalue is complex and has a zero real part. Divergence instability is indicated by a real ~genvalue becoming zero. Having the natural modes available, the calculation of matrix [a of order 2m in eqn (4) is straightforward, only the two upper sub-matrices C* and K* requiring significant computation. In calculating C*, normally C, the structural damping
1323
matrix, is conservatively taken as being zero. The eigensolution of eqn (4) is performed as follows: (a) The real unsymmetric matrix K is reduced to a matrix ii of upper Hessenberg form by similarity tr~sfo~ations. The ~~nvalues of 2 are the same as those of 2. (b) The eigenvalues Izi of a (and hence A, are evaluated using the Double QR method, which is a computationally more efficient extension of Francis’s QR method [ 151. (c) If required, any eigenvector {zj> resending to eigenvalue li is then computed by a back substitution process.
(2)
where the coefficients of [A] and [B] arise from the Theodorsen equations or expe~mentally measured flutter derivatives described above. Applying the modal synthesis technique leads to an m th order system in terms of generalized coordinates
]C*l= Kl’(tC] -
suspension bridges
SOLUTION CHARACTERISTICS
Since the method predicts eigenvalues and vectors defining the response for the given wind speeds, investigation of flutter stability involves a sweep through a range of increasing windspeeds until instability is indicated. Furthermore, the solution of the eigenproblem as in (a)-(c) above at any given wind speed must be iterative in itself. This is because of the aerodynamic terms are functions of response frequency, and re-analysis has to be performed until o,~~, the circular frequency used to evaluate the aerodynamic terms, agrees with the response of interest. Initially, with very low wind speeds, the response frequencies are almost equal to the natural frequencies and the real parts of the eigenvalues are negative, indicating decaying motions. If the response frequencies are plotted on an Argand diagram they all he just to the left of the imaginary axis. When analyses are carried out at successively higher wind speeds, all eigenvalues start to migrate along paths further away from the ima~nary axis, co~s~nding to increasing damping, while oscillation frequencies may increase or decrease from the natural frequency values. Eventually at some higher wind speed, one of these complex eigenvalue paths turn and heads back to intersect the axis. This corresponds to an oscillation with zero ~mping and the onset of flutter instability. Figure 3 shows an example of the response frequencies plotted on an Argand diagram where the flutter speed V, = 76.9 m/s. The characteristics and problems associated with the modal flutter solution have been described in detail by the author elsewhere [16]. NUMERICAL ACCURACY OF FLlJTl’JIR SPEED PREDICTIONS Efect
of basic idealisation
Figure 4 shows a typical structural idealization of the Sevem Bridge where the left, main and right side span decking is represented by 5, 16 and 5 beam elements, respectively. With the properties as shown in Fig. 4, natural frequency analyses in ANSUSP for various idealization refinements yield the natural frequency predictions as indicated in Fig. 5. Table 1
T. J. A. AGAR
1324
Mode 12
0-0-o
m-O-*-
laa++~-a--e_.
!
Mode 10 3.0
12% (““i 2.5
HI-o-e-o-o-o
2.0
1.0
o-e-.-~-
C.
I
I
-0.3
-0.2
Mode 2
c
0
-0.1 IJ
Fig. 3. Argand plot for 5-16-5 idealization of Sevem Bridge with first seven symmetric modes from Fig. 5. The (0) represent plotted points for wind speeds V (m/s): 10, 20, 30, 40, 50, 60, 65, 70, 72.5, 75, 77.5.
shows the CPU times required to compute the lowest 13 and 30 frequencies and mode shapes for the range of idealizations considered. All computation has been performed on a Dee uVAX II computer. Efect of natural modes inclua’ed If a modal eigenvalue analysis is performed, for the 5-16-5 model first including only modes 2 and 8, i.e. the flutter response is considered to be a combination of the fundamental symmetric flexural and torsional modes, then the Argand diagram of Fig. 6 is produced, where the vibration frequency at flutter (1 = 0) follows a path starting at the fundamental torsional frequency while the path of the bending frequency is directly away from the o axis. Including only these two modes, a flutter windspeed V, of -71.0 m/s and flutter frequency W, of 1.545 rad/s is predicted.
Table 1. CPU time for natural frequency analysis for various structural idealizations of the Scvem Bridge
Idealization
Number of degrees of freedom
2-6-2 3-8-3 5-16-5 8-26-8 10-32-10 12-36-12 17-54-17
71 107 215 359 449 561 773
CPU time (set) 13 modes 30 modes 36 55 112 183 240 276 414
140 184 344 627 763 887 1321
An analysis including the additional mode 6 leads to Fig. 7 where flutter is predicted at V, 76.9 m/s and at an oscillation frequency o,= 1.49 rad/s. Now consider an analysis including modes 2, 3, 6, 7, 8, 10, 12, i.e. the lowest 7 symmetric modes from Fig. 4. For the (symmetric) Severn Bridge, the lowest flutter speed corresponds to a symmetric displacement pattern with antisymmetric natural modes not contributing to the response. The Argand diagram is shown in Fig. 3 where the flutter prediction is again at V, = 76.9 m/s and or = 1.495 rad/s. Table 2 shows the flutter speed predictions and modal involvement in the flutter response obtained by including different combinations of modes in an analysis. The conclusion to be drawn from a range of the above types of analysis is that, in this case, the flutter response consists of a combination of modes 2,6 and 8, i.e. the first and third symmetric flexural and fundamental torsional modes. Figure 8 summarizes the flutter speed predictions obtained from the range of idealizations for various mode shapes included in the modal flutter analysis. It can be seen that the flutter speed predictions are (apart from the relatively coarse 2-6-2 and 3-8-3 idealizations) very insensitive to the degree of refinement of the starting structural model. However, two distinct values of flutter speed of N 71 and N 77 m/s are predicted depending on which mode shapes are included in the modal flutter analysis. Figure 9 shows the lowest three symmetric flexural and the lowest symmetric torsional mode shapes. Distance between cables = 22.87
g1’5rI
r[
1 10.2m
5 elements
I-
8 elements 488m
I:
305 m
-I
Structural properties
Cables
: CSA
Hangers
:
Towers
:
Mass CSA Mass Itlong.) W&y)
= = = = = = =
0.158m’ 1339 kg/m 0.0146 m2 124 ka/m 2.13& 8.52 m4 207 x lo6 kN/m2
Deck
: Itvert) J (torsion) Mass Mass moment of inertia Overall width
= = = =
1.079m4 3.995m4 10350 kg/m 696ooo kg-m*/m
=
31.85 m
Fig. 4. Sevem Bridge 5-165 idealization and properties.
Dynamic instability of suspension bridges
1325
Table 2. Flutter speed predictions and modal involvement in response for analyses with different combinations of symmetric mode shapes for 5-16-5 idealization Mode shape
Modal iavolvcmcnt as % flutter response
Flexural 2 3 6 7 10 12
100.0 -
85.4 14.6 -
65.7 34.3 -
57.3 9.8 30.7 2.1 UO 0.1
59.3 9.1 31.6 -
Torsional 8 18
100.0 -
100.0 -
100.0 -
100.0 -
98.6 1.4
q(rad/s) Vhmls) Note -,
1.545 71.0
1.551 70.9
1.495 76.9
1.495 76.9
1.495 76.9
indicates mode not included in analysis.
DISCUSSION
and flexural (mode 2) modes with V,= 77 m/s obtained by including at least modes 2, 6 and 8, there may appear to be some inconsistency. However, on examination the inference from the former analysis is that if the bridge were artificially constrained to allow
Initially, on comparing the flutter prediction of V,= 71 m/s obtained by including only the similarly shaped fundamental symmetric torsional (mode 8) Natural frequency o (radk)
A
_---
Mode
Designation+
-029
3-S-T(C)
*___*_-__*_____--
,*fl-
/
/
1’
8-
+ Mode designation: e.g. 3-S-F(C,S) implies 3rd symmetricflexural mode involving displacements in the centreand side spans.
1’ f
I
I 7-
__~_--.---_-~--- -_-_-
023 2-A-T(C)
,/’ /
/ 6-
.---
./ /
.’
/ _.____--.----.----
.C
.--_-__-_-.21
l-S-T(S) l-A-T(S)
.----_
2-S-T(C,S)
----_.18
l____~_____-~--__~----_.--__-- -413 l-A-T(C) l______._-_.---_.---
&i w l/.--,*-/
- - - -412
/‘-
/
--
,.______ /-----
.---._.--_-_-_-_ _.---.----.--------.g
---_.
4-S-F(C) 11 P-A-F(S) 10 l-S-F(S) *A-F(C)
.@ _*-----He*-----
$ -__ ---
$====$==I==-_==:;
1;-J(y
1 .~.----__-__.__-.-_-_.
.,c
-a_.
------.---.----.--------43
--------a
4
l-A-F(S) 2_S_F(CS)
n
Fig. 5. Natural frequencies of Swem Bridge for various idealizations.
1326
T. J. A. AGAR
t
!2.5
V, = 76.9 mls q = t .49 radls (Wtrial following
jO.5
I
-0.3
CL
wind speeds V (m/s): IO, 20, 30, 40, SO,60, 65, 70, 72.5.
flexural motion
proportional to the fundamental symmetric mode only, flutter could appear as predicted at 71 m/s, However, in reality the response induced by the aerodynamic forces will be the combination of all mode-shapes that the aerodynamics drive. The flexural response at flutter shown in Fig. 10 is consistent with that of Bleich, who reported the effect of incomplete matching between flexural and torsional modes as increasing flutter speed as the initially half-sine wave fundamental flexural mode develops a reverse curvature 1171.
80-
Predicted Flutter Speed V, (m/s)
I
-0.1
0
W
CL
.
Fig. 6. Argand pIot for 5-16-5 idealization of Sevem Bridge with modes 2 and 8. The (0) represent plotted points for
t
I
-0.2
Fig. 7. Argand plot for 5-46-5 idealization of Sevem Bridge with modes 2,6 and 8. The (e) represenl plotted points for wind speeds V(m/s): 10,20,30,40,50,60, 65,70,72.5,75.0, 77.5.
The conclusion is that in order to obtain the correct flutter prediction then modes 2, 6 and 8 at least, must be included. However, as can be seen from Fig, 8 there is no significant change in the flutter speed prediction if additional mode shapes are included. Comparison of the flutter wind speed V, = 76.9 m/s obtained above with predictions by other methods are shown in Table 3. Use of the numerical integration time-history module in ANSUSP (which uses exactly the same formulation for the aerodynamic Modes included + 2,8 (B2,3,8 A 2,6,8 o 2,3,6,8 v 2,3,6,8,18 02,3,6,7,8,10,12,18
Fig. 8. Flutter speed dependence on basic structural idealization and natural modes included.
Dynamic ~s~bitity of suspension bridges Table 3. Comparison of flutter windspeed predictions for the Sevem Bridge by various methods Method
Reference
Multi degree-of-freedom, modal Multi degree-of-freedom, time&story Two degrees-of-freedom, time-history (by author) Two degrees-of-freedom, modal (by author) Two degrees-of-freedom, experimental and calculation (Smith-NPL) Selberg semi-empirical
This method
Flutter speed (m/s) 76.9 77.5
Agar VI -
65.6
-
65.6
Smith [18]
65.7
Selberg [19]
68
force coefficients) produces a value in the region of 77.5 m/s with a flexural response showing a reverse cu~atu~ very similar to that shown in Fig. 10. The other predictions in Table 3 are based on modelling the behaviour with 2 degrees of freedom only and give consistently lower results. The two predictions performed by the author use the same Theodorsen aerodynamic forces as the multi-degree of freedom modal method described here and the time integration procedure. Hence there is a substantiation for the difference in flutter predictions obtained from multi-degree of freedom and two-degree of freedom models.
Mode2
1327
In practice, the static divergence instability is predicted to occur at 74.7 m/s, so that the theoretical increase in flutter speed (from approximately 66 m/s for two-degree of freedom to approximately 77 m/s for multi-degree of freedom methods) is to some extent academic in this case. From the computational efficiency viewpoint, comparison of the modal flutter method described here and the time integration approach shows the former to be much superior. With the modal approach the majority of the computation is involved in the initial natural frequency analysis (e.g. see Table 1: 334 CPU seconds to obtain 30 modes of the 5-16-5 model), while the computation involved in the subsequent modal eigenvalue analysis is modest. It can be shown that an approximate semi-empirical expression for the number of CPU seconds needed for convergence of the modal flutter solution for each windspeed is given by 0.03(171)2~5 where fi is the number of natural modes included in the eigenvalue flutter analysis. By comparison, the multi-degree-of-freedom timehistory integration method with the same S-16-5 structural idealization required N 15OOCPUs per wind speed considered. Con~uently, use of the modal flutter analysis method described in the paper represents a saving in computation by a factor of between 20 and 50 over the time integration approach for comparable accuracy.
(I-S-F)
E
1
~odr3(2-S-F
Moda 6 (3-S-F)
c _---
---_ -I, --.
Mod. 8 (I-S-T
1
-z
-----b--e
Fig. 9. Natural mode shapes 2, 3, 6 and 8 for Sevem Bridge analysis. Tower displacements shown at X5 magnification.
T. J. A. AGAR
1328
0.125
0.25
0.375
0.5
Fig. 10. Deck vertical displacements at flutter from mode shapes 2, 3 and 6. Curves represent displacements at times indicated:
time
Curve
1 2 3 4 T =
CUN’Z
time
t, t + T
5
t + (4T/8)
t + O-/8) t + (2T/8) t + (3T/8)
6 7 8
t + (ST/E) t + (6T/8) 6 + (7T/8)
Period of oscillation = 2x/w,.
CONCLUSIONS A method of computing the flutter speeds of suspension bridges has been presented. The method uses a multi-degree-of-freedom modal technique and includes the aerodynamic force interaction in the governing equations of motion to produce an unsymmetric eigenvalue problem. The characteristic roots of the matrix involved are indicative of the nature of possible motions that can exist at a given wind speed. A study carried out to consider effects of the degree of refinement of the basic structural model and the number of modes included in the eigenvalue flutter analysis has shown that flutter speed predictions are in-sensitive to the basic suspension bridge model, but that it may be important to include more than the fundamental mode shapes on which the simpler 2-degree of freedom methods are based. The modal method described is demonstrated to be computationally much more efficient than an alternative time-history integration method, but yields comparable flutter speed prediction accuracy. REFERENCES 1. Bridge Aerodynamics, Proposed British Design Rules. Thomas Telford, London (1981). 2. T. J. A. Agar, Analysis of suspension bridges-programme ANSUSP user guide. Mott Hay & Anderson Computer Department (1980). 3. E. L. Wilson and R. W. Clough, Dynamic response by a step-by-step analysis. Proc. Symp. on the Use of Computers in Civil Engineering, Lisbon. (1982). 4. K. J. Bathe and E. L. Wilson, Stability and accuracy analysis of direct integration methods. Earthquake Eng. Struct. Dynumic 1, 283-291 (1973). 5. H. M. Hilber, T. J. R. Hughes and R. L. Taylor, Improved numerical dissipation for time integration algorithms in structural dynamics. Earthquake Eng. Struct. Dynamic 5, 283-292 (1977).
6. N. M. Newmark, A method of computation for structural dynamics. Proc. ASCE J. 66, 67-94 (1959). 7. T. J. A. Agar, The analysis of aerodynamics flutter of suspension bridges, Comput. Strut. 30, 593-600 (1988). 8. T. Theodorsen, General theory of aerodynamic instability and the mechanism of flutter, NACA Technical Report, No. 496 (1935). 9. R. H. Scanlan and J. J. Tomko, Airfoil and bridge deck flutter derivatives. Proc. ASCE J. 97, 1717-1737 (1971). 10. T. Okubo and N. Narita, A comparative study of aerodynamic forces acting on cable-stayed bridge girders. Proc. 2nd US-Japan Seminar on Wind Effects on Structures, Kyoto, University of Tokyo Press, pp. 271-283
1
(1976). T. Okubo and K. Yokoyama, Some approaches for improving wind stability of cable-stayed bridges. Proc. 4th Int. Conf. on Wind Efects on Buildings and Structures, London, 1985. Cambridge University Press,
1
pp. 241-249, 1986. H. Loiseau and E. Szechenyi, Etude du comportement aeroelastique du tablier d’un pont a haubans, T.P.
1975-75, Office National d’Etudes et de Reeherches Aerospatiales, Chatillon, France. 13. J. S. Przemieniecki, Theory oj Matrix Structural Analysis. McGraw-Hill, New York (1968). 14. R. B. Corr and A. Jennings, A simultaneous iteration algorithm for symmetric eigenvalue problems. Int. J. Num. Meth Engng 10, 647-663 (1976). 15. J. G. F. Francis, The QR transformation, parts I and II. Computer J. 4, 265-271, 332-345 (1961). 16. T. J. A. Agar, Aerodynamic flutter analysis of suspension bridges by a modal technique. Eng. Structures 11, 75582 (1989). 17. F. Bleich, The flutter theory. The Mathematical Theory of Vibration in Suspension Bridges (Edited by Bleich, F., McCullough, C. B., Rosecrans, R. and Vincent, G. S.), U.S. Government Printing Office, Washing, DC, Ch. 7, pp. 241-281 (1950). 18. I. P. Smith, The aeroelastic stability of the Sevem suspension bridge. NPL Aero Report 1105 (1964). 19. A. Selberg, Oscillation and aerodynamic stability of suspension bridges. Acta Polytech. Stand., Civil Eng. No. 13 (1961).