Integrated finite strip flutter analysis of bridges

Computers and Structures 212 (2019) 145–161 Contents lists available at ScienceDirect Computers and Structures journal homepage: www.elsevier.com/locate/compstruc Integrated finite strip flutter analysis of bridges Hamidreza Naderian a,b,⇑, Moe M.S. Cheung b,c, Majid Mohammadian d, Elena Dragomirescu d a Department of Civil and Mineral Engineering, University of Toronto, Canada Department of Civil and Environmental Engineering, Hong Kong University of Science and Technology, Hong Kong c Western China Earthquake and Hazards Mitigation Research Centre, Sichuan University, Chengdu, China d Department of Civil and Environmental Engineering, University of Ottawa, Canada b a r t i c l e i n f o Article history: Received 14 April 2018 Accepted 10 October 2018 Keywords: Flutter Long-span cable-stayed bridge Aerodynamic Aeroelasticity B3 spline Integrated finite strip method a b s t r a c t The recently developed integrated finite strip method (IFSM) is extended to the aerodynamic flutter analysis of bridges under wind effects. The methodology is capable of performing a three-dimensional (3D) flutter analysis in the IFSM environment, and the solution falls into the category of both multi-mode and full-mode flutter analysis. Aerodynamic stiffness and damping matrices, as well as structural property matrices, are derived by IFSM. In addition, an optimal scheme is proposed for solving the flutter eigenvalue problem. Furthermore, a simple technique has been developed in order to handle different end boundary conditions. The proposed finite strip solution is very straightforward in terms of amount of input data, boundary conditions, modeling, and the flutter analysis process. Moreover, the convergence rate of the method is very high due to the semi-analytical and localization nature of the IFSM. Benchmark numerical investigations are presented, including the study of a simply supported thin flat shell and a model of the Kap Shui Mun Bridge, an existing long-span cable-stayed bridge. The numerical results show that IFSM significantly improves the convergence of the critical flutter frequencies, and therefore leads to smaller storage requirements and faster flutter eigenvalue extraction. Ó 2018 Elsevier Ltd. All rights reserved. 1. Introduction One of the main challenges in structures surrounded by wind flow is aerodynamic instability due to the effects of self-excited forces. If the wind-structure interaction leads to an increase in the magnitude of oscillatory motions, then aeroelastic instability occurs. This is caused by aerodynamic forces induced on the structure as a consequence of its motion. Such aerodynamic forces are called self-excited forces, and this self-excited oscillatory instability phenomenon is called flutter. Some examples of structures critical in countering flutter are airfoils in the aerospace industry like aircraft wings, and long-span cable-stayed bridges in civil engineering. In flutter, at a critical wind speed, the self-excited aerodynamic forces acting on an oscillatory system start feeding the energy of the system instead of dissipating it, which leads to destructive divergent vibrations in the structure. The aerodynamic behavior of a structure depends on the interaction between the aerodynamic forces of inertia and damping, and the structure’s stiffness properties. Therefore, a multi-disciplinary study between ⇑ Corresponding author at: 35 Saint George Street, Dept. of Civil and Mineral Engineering University of Toronto, Toronto, ON, M5S 1A4, Canada. E-mail address: [email protected] (H. Naderian). https://doi.org/10.1016/j.compstruc.2018.10.003 0045-7949/Ó 2018 Elsevier Ltd. All rights reserved. solid-mechanics, dynamics, and fluid-dynamics can investigate the aeroelasticity of an airfoil or a long-span bridge. Modern long-span bridges are usually constructed with very slender deck structures, which results in an increased risk of instability under extreme dynamic loads like wind and earthquakes. With the advent of new advanced materials such as composite fiber-reinforced polymers (FRP), longer-span bridges are being built throughout the world. On the other hand, these types of materials are very lightweight, and consequently, the stability of bridge structures against wind effects has become more critical. The nature of wind load is very complex and turbulent, which means it is not easy to understand the structural behavior of bridges under wind flow. Flutter is one of the most important criteria in the analysis and design of suspension and long-span cable-stayed bridges. The determination of the dynamic and aerodynamic characteristics of a particular bridge is a crucial step in solving the flutter problem. Since the 1960s, a number of publications have dealt with both analytical and experimental studies on flutter in long-span bridges. Tanaka et al. [1], and Ge et al. [2] investigated flutter in plates and bridges by multi-mode and fullmode approaches. Huang et al. [3] studied the aerodynamic behavior of bridges under skew wind by performing a series of section model tests. Dowell et al. [4] investigated the flutter of rectangular plates under a three-dimensional axial flow through numerical and 146 H. Naderian et al. / Computers and Structures 212 (2019) 145–161 experimental tests, and Wu et al. [5] presented a framework for linear and nonlinear aeroelastic analyses of cable-supported bridges. The topic of bridge flutter is still highly important, and researchers have tried to develop different types of experimental and analytical techniques for a better understanding of the flutter phenomenon [6–12]. The finite strip method (FSM), pioneered by Y. K. Cheung and M. S. Cheung [13,14], is one of the most efficient and accurate numerical techniques for modelling long-span structures with highly harmonic behaviour along the longitudinal direction. Among these, one can mention buckling and free vibration of thin-walled structures, as well as the structural investigation of long-span bridges. Despite the merits of the finite strip method, its application to bridge engineering has been limited to modeling only the bridge deck, while other structural components of the structure, including piers, towers, cables, etc., have been modeled as special boundary conditions of the deck [13]. For that reason, and with the advancement of other numerical techniques, the popularity of the finite strip method among bridge engineers has declined in recent decades. In order to overcome this obstacle, Naderian et al. [15,16] developed the so-called Integrated Finite Strip Method (IFSM), which is able to model an entire cable-stayed bridge system. IFSM is, in fact, an evolved version of the Spline Finite Strip Method (SFSM) [17], which is itself a more comprehensive format of the conventional Finite Strip Method (FSM), at the expense of incorporating more degrees of freedom. However, as a tool, this is more powerful than the ordinary finite strip method in terms of dealing with structures with complex geometry, boundary, and loading conditions. Moreover, SFSM is also capable of modeling multispan bridges and structures subjected to not only uniformly distributed forces but also concentrated forces. In IFSM, the bridge deck is divided into a number of flat shell spline finite strips in the longitudinal direction, and B3 spline functions define the longitudinal shape function of the displacement function. Using B3 splines, the continuity over the second derivative of the displacement function is achieved. Furthermore, the splines have localization properties through which the bandwidth of the global property matrices becomes very narrow. As a result, the computational time can be dramatically decreased. When it comes to complicated analyses, for instance, vibration-based real-time structural health monitoring of huge bridge structures under earthquake, typhoon, or high traffic loads, where numerous iterative dynamic analyses must be performed consistently, even a small percentage of computational savings is highly appreciated. Therefore, IFSM can be a suitable platform for the vibration investigation of modern bridge structures. Naderian et al. [16] recently extended the finite strip method to the seismic analysis of long-span cable-stayed bridges. In this paper, the IFSM is further extended to aerodynamic flutter analysis. It is worth mentioning that M. S. Cheung et al. [18] used a combination of the spline finite strip method and the finite element method for flutter analysis of long-span cable-stayed bridges. Despite the merits of their method, there are still considerable concerns about the efficiency and accuracy of the solution, which are addressed in the following: (I) the aerodynamic stiffness and damping matrices are not derived by the spline finite strip method but are based on lumped aerodynamic forces. The distribution of the aerodynamic forces along the nodal lines of the strip is based on trapezoidal discretization, while in the spline finite strip method it must be done by B3 spline discretization; (II) Only the deck is modeled by spline finite strips, while other components like piers, towers, and cables are simulated by equivalent beam finite elements. In other words, their numerical solution is a partial finite strip solution; and (III) The aerodynamic self-excited forces are spatially distributed on the deck of the bridge, while the general formulation of the self-excited flutter forces is based on applying the self-excited forces on the center of elasticity of the deck [19]. According to the present research, however, the spatial distribution of self-excited aerodynamic forces is not only inefficient, but also it might be against general assumptions of the flutter formulation, which might lead to less accurate results. Although the idea behind the uniform distribution of wind loads is to consider the effects of the spatial distribution of aerodynamic forces over the cross-section of the bridge deck, our recent studies indicate that the aerodynamic performance is not sensitive with regard to this factor. Therefore, spatial discretization of aerodynamic loads on the entire deck system decreases the efficiency of the solution procedure. The latter will be discussed more in detail later in this paper when providing the numerical examples. Last but not least, Cheung et al.’s [18] proposed method might not perfectly and efficiently consider the structural interactions between decks, piers, towers, and cables because the compatibility between the displacements of the interaction joints is obtained through an iterative process. However, this technique could be effective for simple types of structures under static/quasi-static forces. In the case of complex structures like long-span cable-stayed bridges or structures under dynamic and aerodynamic excitations such as non-uniform seismic excitations or aerodynamic self-excited forces, such simplification can no longer be adopted. All the above-mentioned concerns have been eliminated in the present study by introducing the Integrated Finite Strip Flutter formulation, as follows: the shell spline finite strip models the bridge deck, the one-dimensional column strip models the towers and piers, and the cable strip models the bridge cables. Transition section elements in each of these strips provide the connectivity and boundary conditions between different structural components of a cable-stayed bridge. In this research, first the IFSM is briefly reviewed, and then the general formulation of the flutter equation of motion is presented. The integrated finite strip flutter methodology along with the derivation of the aerodynamic stiffness and damping matrices, as well as the structural property matrices, are fully explained. In addition, a simplified method has been developed to handle the boundary conditions corresponding to the flutter problem. Thirdly, multi-mode and full-mode aeroelastic approaches are integrated within the finite strip treatment. Finally, a simple eigenvalue analysis based on linearization of the flutter equation of motion is proposed through which the critical flutter frequencies and consequently critical wind speeds can be obtained by solving a standard eigenvalue analysis. The proposed finite strip flutter solution is applied to two numerical examples where the flutter behavior of a thin-walled flat plate, as well as the long-span cable-stayed Kap Shui Mun Bridge, are investigated. The accuracy of the method in predicting critical flutter speeds and frequencies, as well as the convergence of the currently developed finite strip method for flutter analysis, is examined. The results show that IFSM and the proposed flutter analysis scheme is a powerful and efficient tool for the aerodynamic analysis of a plate in general and for bridge structures as a system. Rapid convergence and small storage requirements are highlights of the proposed method. 2. Review on integrated finite strip method The Integrated Finite Strip Method (IFSM) [15,16] is the latest version of the finite strip method and is able to fully simulate and analyze the bending and vibration of long-span cable-stayed bridges. Not only the bridge deck but also all other structural parts of a long-span cable-stayed bridge including towers, piers, linked beams, and cables, as well as bearings, can be modeled in the environment of the spline finite strip method by using different types 147 H. Naderian et al. / Computers and Structures 212 (2019) 145–161 of strips, including shell spline strips, columns strips, and cable strips. The key parameter of the originality of the IFSM is the socalled transition section, which is based on unequal spline strips. Transition sections provide compatibility of displacements and vibrations in the intersection joints of different components, and they also handle the boundary conditions of bearings. The compatibility can be achieved by selecting an extremely small length for a spline transition section [15]. Another important feature of the IFSM is that the structural interactions between different components, for example between towers and cables, can be modelled. This capability is very important, especially for dynamic analyses of long-span bridges like seismic analysis and flutter analysis, in which there is considerable interaction and load effect transmission between different structural elements. More details on the development of the IFSM have been addressed by Naderian et al. [15]. 3. Aerodynamic flutter equation of motion The self-excited aeroelastic forces acting on a plate or a bridge structure are caused by wind effects. This means that flutterresultant forces can be considered as external forces. Similar to general dynamic systems, the equilibrium of the effective forces, including the inertia force, the damping force, the elastic force, and the applied flutter force, will form the flutter equation of the motion. For a structure for which displacements are defined by a vector fdg including n degrees of freedom, matrix analysis can be used to formulate the flutter equation of motion as: n o n o ½Ms Š d€ þ ½C s Š d_ þ ½K s Šfdg ¼ fF aer g ð1Þ where ½Ms Š, ½C s Š, and ½K s Š are the global structural mass, damping, and stiffness matrices respectively with the order of n  n, while n o fF aer g is the global self-excited aerodynamic force vector, and d_ n o and € d are the vectors of the first and second derivatives of dis- placement respectively. The self-excited forces fF aer g are dependent on the vertical (w), lateral (u), and rotation ðhÞ displacements of the structure and their first derivatives, and it can be expressed as: n o fF aer g ¼ ½C aer Š d_ þ ½K aer Šfdg ð2Þ in which ½C aer Š and ½K aer Š are the n  n global aerodynamic damping and stiffness matrices respectively. The self-excited aerodynamic force associated with the second derivative of displacement, called aeroelastic inertial forces, is omitted, as the air density is significantly lower than that of the structural materials. Substituting Eq. (2) into Eq. (1) yields the flutter equation of motion as: n o n o n o ½Ms Š d€ þ ½C s Š d_ þ ½K s Šfdg ¼ ½C aer Š d_ þ ½K aer Šfdg ð3Þ 4.1. General Both the structural and aerodynamic properties of different components of a bridge can be derived by the integrated finite strip method. In the case of a long-span cable-stayed bridge, only the self-excited flutter forces acting on the bridge deck are considered in the analysis. In the following sub-section, the integrated finite strip discretization as well as the derivation of the structural and aerodynamic matrices in the environment of IFSM are presented. 4.2. Integrated finite strip discretization Considering the bridge deck as a shell structure, both in-plane and out-of-plane degrees of freedom need to be considered. A flat shell spline strip is shown in Fig. 1, where each knot of a nodal line has four degrees of freedom, three translational (u, v, w), and one rotation ðhÞ. The vector of the displacement parameters of a shell spline strip centered at ym is given as: fdgm ¼ ½uim ; mim ; wim ; him ; ujm ; wjm ; hjm ŠT ð4Þ in which ½MŠ ¼ ½Ms Š, ½C Š ¼ ½C s Š ½C aer Š, and ½K Š ¼ ½K s Š ½K aer Š, where ½M Š, ½C Š, ½K Š are the system mass, damping, and stiffness matrices respectively. The aerodynamic damping ½C aer Š and stiffness matrices ½K aer Š are not necessarily symmetrical because of the nonconservative nature of the aerodynamic forces. Therefore, the system stiffness ½K Š and system damping ½C Š matrices are asymmetric, which results in coupling effects between the flutter response frequencies. In other words, the asymmetry of the aerodynamic properties enables the flutter response of the structures to be coupled between multiple mode shapes. ð5Þ It must be noted that the parameters of Eq. (5) are in fact spline coefficients and not the physical degrees of freedom. The integrated finite strip method is based on the use of unequally spaced B3-spline functions, through which the transition section is developed for connecting different structural components together. Moreover, it allows the locations of the supports and the concentrated load to coincide with the knots on the nodal lines, thus leading to more accurate results. Moreover, the introduction of unequally spaced interior knots allows one to describe accurate responses in the region of high stress gradients, or at the locations of abrupt geometric changes, by spacing the knots more closely. In a flutter analysis problem of a flat shell or a bridge deck, the unequally spaced B3-spline function centered at ym may define the longitudinal shape of function as: 8 0 > > > > > > Am ðy ym > > > > < Am ðy y m Um ðyÞ ¼ > Bm ðymþ2 > > > > > B ðy > > m mþ2 > > : 0 y < ym 3 2Þ 3 2Þ 3 ym 1 Þ3 þ C m ðy yÞ þ Dm ðymþ1 yÞ3 yÞ3 2 ym 2 6 y < ym ym 1 6 y < ym 1 ð6Þ ym 6 y < ymþ1 ymþ1 6 y < ymþ2 ymþ2 6 y in which  Am ¼ ðymþ1  Bm ¼ ðymþ2 C m ¼ ðymþ2 Dm ¼ ðymþ2 or n o n o ½MŠ €d þ ½C Š d_ þ ½K Šfdg ¼ 0 4. Integrated finite strip flutter formulation ym 2 Þðym ym 2 Þðym 1  ym 2 Þ  ymþ1 Þ ym 1 Þðymþ2 ym Þðymþ2  ym 2 Þ ðymþ2 ym 1 Þðymþ1  ym 2 Þ ðymþ1 ym 2 Þðymþ1 1 1 ym 1 Þðym ym 1 Þðymþ1 ym 1 Þðym 1 ym Þðymþ2  ym 2 Þ  ymþ1 Þ 1 1 ð7Þ The section lengths of the B3 spline function can be chosen to be equal for flutter analysis of simple cases like a long-span flat shell where there is no need to employ unequally spaced splines. The length of the plate strip L is divided into m sections of equal length h. A typical local B3-spline function of equal spaced is defined as: 8 ðy ym 2 Þ3 > > > > 3 2 2 3 > > > h þ 3h ðy ym 1 Þ þ 3hðy ym 1 Þ 3ðy ym 1 Þ 1 <

Um ðyÞ ¼ 3 h3 þ 3h2 ymþ1 y þ 3hðymþ1 yÞ2 3 ymþ1 y 3 6h > > > > ðy > yÞ3 > mþ2 > : 0 ym 2 6y6y ym 1 6 y 6 ym ym 1 ym 6 y 6 ymþ1 ymþ1 6 y 6 ymþ2 otherwise ð8Þ 148 H. Naderian et al. / Computers and Structures 212 (2019) 145–161 Fig. 1. Shell spline finite strip. The membrane displacement functions u and v and the flexural displacement function w are expressed as the product of transverse polynomials and longitudinal B3-splines: u¼ rþ1 X ðN 1 U1m ðyÞuim þ N2 U5m ðyÞujm Þ ð9Þ ½Uim Š ¼ ½ U 1 U0 U1 U2 .. .. . . Um 2 Um 1 Um Umþ1 Š ð14Þ ½U1m Š, ½U2m Š, ½U5m Š and ½U6m Š are related to the displacements u and v of nodal lines i and j respectively, while ½U3m Š, ½U4m Š, ½U7m Š and ½U8m Š are related to the displacement w. In the short form, Eq. (12) is expressed as: m¼ 1 v¼ rþ1 X ff g ¼ ½Nt Š½UŠfdg ðN1 U2m ðyÞv im þ N2 U6m ðyÞv jm Þ ð10Þ m¼ 1 w¼ rþ1 X ðN3 U3m ðyÞwim þ N 4 U4m ðyÞhim þ N 5 U7m ðyÞwjm m¼ 1 þ N6 U8m ðyÞhjm Þ ð11Þ where r is the total number of longitudinal sections on a nodal line. In the matrix form, Eqs. (9)–(11) can be rewritten as: 8 9 2 3 N1 0 0 0 N2 0 0 0 > = <u> 6 7 ff g ¼ v ¼ 4 0 N1 0 0 0 N2 0 0 5 > ; : > 0 0 N3 N4 0 0 N5 N6 w 2 38 9 ½U1m Š > > uim > > > > 6 7> > > 6 7> ½U2m Š > > v im > > 6 7> > > > 6 7 > > ½U3m Š 6 > wim > 7> > > > > 6 7> > = 6 7< him > ½U4m Š 6 7  6 7 > 6 7> ½U5m Š > ujm > > 6 7> > > > 6 7> > > v 6 7> ½U6m Š > > jm > > 6 7> > > >w > 6 7 > > ½U7m Š jm > 4 5> > > > > ; : hjm ½U8m Š ð15Þ where ½N t Š is a 3  8 matrix and ½UŠ is an 8  8 (m + 3) matrix, while fdg is an 8 (m + 3)1. Thus, fdg is an 8 (m + 3)1 matrix. 4.2.1. Column and cable strips The cantilever behavior of the piers and towers of a long-span cable-stayed bridge can be simulated by the so-called threedimensional (3D) Column Strip (CS3), which is a vertical shell spline strip fixed at one end, for providing the support boundary conditions, and free at the other end, as shown in Fig. 2. In some cases, a one-dimensional (1D) column strip (CS1) is more applicable and is actually easier to be employed for modelling piers and towers. In a CS1, each knot belonging to the nodal line has three translational degrees of freedom. Consequently, the displacement function can be defined as: u¼ rþ1 X um Um ðzÞ ð16Þ m¼ 1 w¼ rþ1 X wm Um ðzÞ ð17Þ m¼ 1 ð12Þ The transverse shape functions adopted in Eq. (12) are cubic Hermite polynomial functions for vertical displacement variation and linear interpolation for in-plane displacements as: N1 ¼ 1 X; N2 ¼ X; N3 ¼ 1 N 5 ¼ ð3X 2 2X Þ 2 XÞ N 6 ¼ xðX 3X 2 þ 2X 3 ; N4 ¼ xð1 2X þ X 2 Þ; 3 ð13Þ In Eq. (12), ½U1m Š to ½U8m Š are the row matrices, and each of them includes (m + 3) local B3-splines as: Fig. 2. Column Strip (CS3). 149 H. Naderian et al. / Computers and Structures 212 (2019) 145–161 v¼ rþ1 X v m Um ðzÞ ð18Þ m¼ 1 For the displacement-strain relationships, only the bending in the vertical and transverse directions and the axial stress are considered, while the shear stress and torsional moment are assumed to be negligible, as the amounts of these forces are very low in CS1. In order to model the bridge cables in the finite strip environment, the cable strip is employed. The cable strip is a simplified version of the CS1. Generally, cables can resist only axial tension stresses and not compression. For this reason, in the cable strip model, only the axial stress defines the strain-displacement relationship. 4.2.2. Connectivity provisions In contrast with the finite element method in which auxiliary bridge structural elements like piers and cables can be easily connected to the bridge deck, the concept ‘‘element” in the longitudinal direction is not defined in the conventional finite strip methods. In other words, it is impossible to insert an extra component at an intermediate point within a spline strip. In order to eliminate this impediment, a special transition section has been developed within the IFSM [15]. The methodology has been applied to connect different components of complex structures such as long-span cable-stayed bridges [15]. Therefore, deck-pier and deck-tower connectivity as well as cable-deck and cabletower connectivity are modelled by employing these newly developed transition section elements. The transition section is created by using unequally spaced B3spline functions. A typical transition section connecting two different components is shown in Fig. 3. Without losing generality, it is assumed that the length of the normal and the transition sections are H and Ht respectively. In Fig. 3, the vertical line is a nodal line on the column strip, while the horizontal line is a nodal line on the shell spline strip of the deck. This schematic can also be applied to the cable strip and shell spline strip joints in order to ensure connectivity. The vertical and horizontal lines overlap at knots 3 and 8 of the shell and column strips, respectively, and in order to achieve identical displacements at knots 3 and 8, the ratio of Ht =H should be infinitely small. Depending on the complexity of the structure, different values of h can meet the required accuracy. Numerical studies have shown that under normal circumstances good accuracy can be obtained with the ratio of Ht /H = 0.001 [15]. Using the developed transition section in the spline finite strip procedure, compatibility for the displacements of different components of the structure is satisfied. It must be noted that the bearings can be modeled as special boundary conditions for the transition section. For example, in order to model a fixed bearing, which allows rotation but restricts translations, the displacements of knots 3 and 8 must be equal in order to achieve compatibility, which can be implemented through the transition section as detailed above. 4.3. Finite strip structural property matrices After the discretization of the entire bridge system is performed using the IFSM, the principle of minimum potential energy can be applied to derive the mass matrix [m] and stiffness matrix [k] of the shell strip, column strip, and cable strip, as follows: ½mŠ ¼ ½kŠ ¼ Z Z ½NŠTi q½NŠi dV ð19Þ ½BŠTi ½DŠ½BŠi dV ð20Þ in which q is the density of the strip, [D] and [B] are the elastic matrix and the strain matrix respectively, and [N] is the shape function matrix, which is the multiplication of the transverse functions and the longitudinal spline functions as ½Nt Š½UŠ. It must be noted that in the case of a one-dimensional column strip (CS1) and a cable strip, the integral is performed over the length of the strip. Similar to the FEM, the IFSM strip properties are converted to the knots of the nodal lines along the strip during the model formulation process. However, the number of knots required is significantly reduced due to the semi-analytical nature of the IFSM. The structural matrices of the strips can be assembled using the standard assembly techniques, and the global structural stiffness matrix ½K s Š and global structural mass matrix ½Ms Š of the entire bridge system are formed after assembling all structural components. The structural damping matrix ½C s Š is usually a function of the structural stiffness and mass matrices, and it can be obtained from the following formula, referred to as ‘‘classic Rayleigh damping”: ½C s Š ¼ a½Ms Š þ b½K s Š ð21Þ where a and b are the Rayleigh damping factors, which can be determined by having two structural damping ratios n associated with two specific frequencies. The damping ratio n of the nth free vibration mode of a structure is given by: nn ¼ 1 x aþ nb 2xn 2 ð22Þ where xn is the angular frequency of the nth mode. Assuming the same damping ratio n for two different modes, a and b can be calculated. 4.4. Finite strip aerodynamic property matrices 4.4.1. Aerodynamic self-excited forces In order to derive the aerodynamic damping and stiffness matrices of a flat shell or a bridge deck using the finite strip method, first, the self-excited flutter forces need to be recognized. According to the linear aerodynamics, three types of uniformly distributed line loads acting along the center of elasticity of a moving structure, including lift force Lf , drag force Df , and pitching moments M f , as illustrated in Fig. 4, define the self-excited flutter forces per unit length of the structure as: Lf ¼ Fig. 3. Typical transition section element.   1 w0 Bh0 w Bu0 u qair U 2 B K w H1 þ K w H2 þ K 2w H3 h þ K 2w H4 þ K w H5 þ K 2w H6 U U U 2 B B ð23Þ   1 u0 Bh0 u w0 w þ K 2w P 3 h þ K 2w P 4 þ K w P 5 þ K 2w P 6 Df ¼ qair U 2 B K w P 1 þ K w P2 U U U 2 B B ð24Þ 150 H. Naderian et al. / Computers and Structures 212 (2019) 145–161 0 0 0 1 w  Bh  w u u M f ¼ qair U 2 B2 K w A1 þ K w A2 þ K 2w A3 h þ K 2w A4 þ K w A5 þ K 2w A6 U U U 2 B B   ð25Þ where u, w, and h are the vertical bending, lateral bending, and torsional displacements, as shown in Fig. 5, while u0 , w0 ; and h0 are the first derivatives of the corresponding displacements. K w is the nondimensional reduced frequency given by K w ¼ Bxf =U, where xf is the flutter angular frequency xf , B is the deck width, and U is the mean wind velocity. The terms Hi , P i , and Ai (i = 1, 2, 3, . . . , 6) are the aerodynamic derivatives associated with the self-excited lift Lf , drag Df , and moment Mf aerodynamic forces respectively. The flutter derivatives are all functions of the reduced frequency K w and are usually measured from experimental tests performed on the cross-section of the structure. Also, qair refers to the air density. When the cross-section of the deck is a slab girder, the center of elasticity is equal to the geometric shear center, while in the case of a box girder bridge, the aerodynamic forces need to be mapped on the nodal lines of the strip. This procedure, however, is out of the scope of the current research, as the deck of the bridge is assumed to be a slab girder. Nevertheless, this can be achieved simply by adjusting the displacement functions. Also, it must be noted that the wind velocity acting on the cross-section of the plate (deck) might come from different angles of attack. It is the responsibility of the aerodynamic design engineer to investigate which wind attack angle(s) can cause the maximum aerodynamic force and subsequently instability of the structure. 4.4.2. Spline discretization of aerodynamic forces In conventional analytical methods and the finite element method, the uniformly self-excited forces are linearly lumped along the structure, whereby one-half of each element load acts upon each element-end. Herein, the B3-spline interpolation is employed to distribute the aerodynamic forces along the plate. The spline discretization of the aerodynamic forces leads to a more accurate distribution because of the higher continuity degree (C2). In the finite strip formulation, it is more convenient to locate a nodal line passing through the aerodynamic line loads, acting on the centre of elasticity of the deck. Applying the principle of minimum potential energy in the IFSM environment, the self-excited aerodynamic forces originating from the wind velocity as an external wind-induced force can be distributed on the knots along the nodal line of a spline strip. Integrating the forces over the elements, one obtains: ff aer g ¼ Z T ½NŠ qdA ð26Þ where q represents the self-excited aerodynamic forces given by: q ¼ Lf þ Df þ M f ð27Þ Fig. 5. Aerodynamic displacement components. Substituting Eq. (27) into (26) yields: ff aer g ¼ Z ½/ŠT ½Nt ŠT Lf dA þ ½/ŠT ½Nt ŠT Df dA Z þ ½/ŠT ½Nt ŠT M f dA Z ð28Þ where ½/ŠT is the longitudinal B3 spline function of Eqs. (6) or (8). ½N t Š is the transverse shape function of Eq. (13), where its integraRb tion 0 ½Nt ŠT dx over the width of strip b corresponding to the aerodynamic forces Lf , Df , and Mf is given by the following matrices respectively: 2 0 0 0 0 0 0 0 0 3 7 6 40 0 0 0 0 0 0 05 0 0 1 0 0 0 0 0 3 1 0 0 0 0 0 0 0 7 6 40 1 0 0 0 0 0 05 0 0 0 0 0 0 0 0 2 2 0 0 0 0 0 0 0 0 3 7 6 40 0 0 0 0 0 0 05 0 0 0 1 0 0 0 0 ð29Þ ð30Þ ð31Þ The integration of the spline representations of Eqs. (6) or (8) over the length of the strip L for each degree of freedom (displacement parameter) results in m + 3 values corresponding to the knots  of a nodal line. This integral can be called /f . Assuming the Fig. 4. Self-excited aerodynamic forces [1]. 151 H. Naderian et al. / Computers and Structures 212 (2019) 145–161 location of the aerodynamic loads on the first nodal line of the strip i, and after mathematical rearrangement, Eq. (28) is expanded to the following matrix formulation: P 4 ½/fdiag Š 6 6 ½0Š 6  6 H ½/fdiag Š 6 6 6  6 1 2 2 6 BA6 ½/fdiag Š ff aer g ¼ qair U K w 6 2 6 6 6 6 6 4 2 P 1 ½/fdiag Š 6 6 ½0Š 6  6 H ½/fdiag Š 6 5 6  6 BA5 ½/fdiag Š 1 þ qair UBK w 6 6 2 6 6 6 6 6 4 2 ½0Š P 6 ½/fdiag Š ½0Š ½0Š ½0Š H4 ½/fdiag Š BP3 ½/fdiag Š ½0Š BH3 ½/fdiag Š ½0Š BA4 ½/fdiag Š B2 A3 ½/fdiag Š ½0Š P 5 ½/fdiag Š ½0Š ½0Š ½0Š H1 ½/fdiag Š ½0Š BA1 ½/fdiag Š B2 A2 ½/fdiag Š ð32Þ    in which /fdiag is the diagonal form of /f + 3)  (m + 3), and is given by: 2   6 /fdiag ¼ 6 4 /f with the order of (m 3 1 .. BP 2 ½/fdiag Š ½0Š BH2 ½/fdiag Š 38 9 u > > i > > 7> > > 7> v > i > > 7> > > > > 7> > wi > 7> > > > > > 7< > 7 hi = 7 7> u > j > 7> > > > 7> > vj > > 7> > > 7> > > > > 7> > w > j 5> > > > : ; hj 38 0 9 > > ui > > 7> > v0 > > 7> > i> > 7> > 0> > > > 7> > w > i> 7> > > > > 7> 0 7< hi = 7 7> u0 > 7> > j > > > 7> > v 0j > > 7> > > > 7> > 0> > 7> > > w j> 5> > > 0 > > : ; hj . /fmþ1 7 7 5 ð33Þ 4.4.3. Finite strip aerodynamic stiffness and damping matrices The first and second matrices in Eq. (32) are the aerodynamic stiffness and damping matrices of a spline strip respectively. As there are m + 3 knots on each nodal line, the number of each displacement parameter in Eq. (32) will be m + 3. Applying the spline finite strip derivation, the aerodynamic stiffness and damping matrix of a shell spline strip can be obtained by: 2 P4 ½/fdiag Š ½0Š P 6 ½/fdiag Š BP3 ½/fdiag Š 6 ½0Š ½0Š ½0Š 6 ½0Š 6  6 H ½/ Š ½0Š H ½/ Š BH ½/ Š 6 6 fdiag fdiag fdiag 4 3 6 6 BA ½/ Š ½0Š BA ½/ Š B2 A ½/ Š 1 6 fdiag fdiag 4 3 ½kaer Š ¼ qair U 2 K 2w 6 6 fdiag 6 2 6 6 6 6 6 4 P 1 ½/fdiag Š 6 6 ½0Š 6  6 H ½/fdiag Š 6 5 6  6 BA5 ½/fdiag Š 1 ½caer Š ¼ qair UBK w 6 6 2 6 6 6 6 6 4 2 ½0Š P 5 ½/fdiag Š BP 2 ½/fdiag Š ½0Š ½0Š ½0Š ½0Š H1 ½/fdiag Š BH2 ½/fdiag Š ½0Š BA1 ½/fdiag Š B2 A2 ½/fdiag Š   3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 ð34Þ 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 ð35Þ Except for the spline strip with self-excited forces, the aerodynamic stiffness and damping matrices of the rest of the strips in the entire IFSM model is zero, as there is no self-excited force acting on other elements of the model. The global aerodynamic stiffness ½K aer Š and damping ½C aer Š matrices are formed using the conventional assembly techniques. Finally, the aerodynamic equation of motion for the entire structural system in the IFSM environment is formulated in the general form of Eq. (2). 4.5. Boundary conditions The strip in the FSM must be defined with preset boundary conditions. There are several techniques that can be applied to model the boundary conditions of the strip [13]. Despite the advantages of the spline finite strip method over the finite element method in terms of computational efficiency, handling of a complex amended scheme of local splines for considering the end and internal boundary conditions, keeps this solution untidy. In other words, current amendment schemes for boundary conditions are unable to be generalized, and dealing with boundary conditions using standard techniques like penalty functions is complicated. Herein, a straightforward method for modelling boundary conditions based on replacing the spline displacement parameters with physical degrees of freedom is proposed. This will result in a general unified formulation of otherwise very complex and tedious amended schemes of local splines in the vicinity of the boundary supports and at any internal support. This makes the method more versatile and adjustable with other numerical techniques. For an equally spaced B3-spline function, one can rewrite the physical displacements uphy ; v phy ; wphy and the rotation hphy at the end knots of a nodal line i of a spline strip as: 1 ðu 1 þ 4u0 þ u1 Þ 6 1 v phy0 ¼ ðv 1 þ 4v 0 þ v 1 Þ 6 1 wphy0 ¼ ðw 1 þ 4w0 þ w1 Þ 6 1 hphy0 ¼ ð h 1 þ h1 Þ 2h 1 uphym ¼ ðum 1 þ 4um þ umþ1 Þ 6 1 v phym ¼ ðv m 1 þ 4v m þ v mþ1 Þ 6 1 wphym ¼ ðwm 1 þ 4wm þ wmþ1 Þ 6 1 hphym ¼ ð hm 1 þ hmþ1 Þ 2h uphy0 ¼ ð36Þ where h is the spline section length. Dividing a spline strip into m equal sections will give 8(m + 3) spline parameters which define the displacement function of the shell spline strip. Similarly, the displacements and rotations of the intermediate knots are also presented by spline parameters. The transformation matrix corresponding to each degree of freedom ui v i wi hi uj v j wj hj is expressed by: 2 1 2h 6 61 66 6 6 6 6 6 6 6 6 6 6 T¼6 6 6 6 6 6 6 6 6 6 6 6 4 0 1 2h 4 6 1 6 1 6 4 6 3 1 6 .. . 1 6 4 6 1 6 1 6 4 6 1 6 1 2h 0 1 2h 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 ð37Þ 152 H. Naderian et al. / Computers and Structures 212 (2019) 145–161 where T is a 8  8 (m + 3) matrix. The global transformation matrix of a flat shell spline strip is presented as: T 2 3 T 6 6 6 6 6 6 6 Tt ¼ 6 6 6 6 6 6 4 T T T T T T 7 7 7 7 7 7 7 7 7 7 7 7 7 5 ð38Þ ½K st Š ¼ ½T t ŠT ½K s Š½T t Š ð39Þ ½Mst Š ¼ ½T t ŠT ½M s Š½T t Š ð40Þ After transferring all degrees of freedom to the physical coordinate system, for each restrained degree of freedom a corresponding zero value is imposed on the physical displacement vector, which means the corresponding rows and columns in the structural and aerodynamic properties matrices are eliminated. The proposed amended scheme is more applicable to the equal section splines, while for the unequal section splines the penalty function is used. 5. Flutter eigenvalue analysis Eq. (4) is an eigenvalue problem which can be solved using the techniques of either frequency domain or time domain analyses. In the present research, the flutter eigenvalue problem is solved in the environment of the frequency domain method. The flutter response frequencies of the structure are primarily dependent on the natural frequencies of the system, and different schemes can be used for solving the flutter equation like PK, PK-F method, Jacobi diagonalization, QL or QR transformation, subspace iteration, etc. Herein, a simple method is proposed which is based on linearizing the quadratic flutter Eq. (4). The philosophy of linearization of the flutter eigenvalue problem is to allow one to use the eigenvalue n o solution of an undamped system. Assuming z ¼ d_ , one can reformulate Eq. (4) as: )    ( z_ ½C Š ½K Š z ½M Š ½0Š n o ¼ d_ ½I Š ½0Š fdg ½0Š ½IŠ z_ o n d_ ) 8n o9 > < €d > =  kfwg  ¼ n o ¼ kekt > fwg : d_ > ; ð47Þ In the full-mode flutter analysis, the flutter eigenvalues of an original n degrees of freedom system, rather than an m order system modified by m natural modes, is calculated. The full-mode aeroelastic analysis employs the vector iteration method to solve partial eigenvalues of Eq. (43). The iterative analysis must converge to the most dominant value of k, which is equal to the highest frequency xf . For a complicated structural system like a long-span cable-stayed bridge with a massive number of degrees of freedom, full-mode eigenvalue flutter analysis requires a large amount of time, which is computationally very inefficient. Research has shown that difference between the results of full-mode flutter analysis and multi-mode analysis is slight [1]. The flutter phenomenon usually occurs under the effect of the first few torsional and/or heave modes of the structure. For that reason, full-mode analysis might not be considered as efficient for the eigenvalue flutter analysis. 5.2. Multi-mode finite strip aeroelastic analysis A globally accepted algorithm for solving the eigenvalue flutter problem, especially for linear elastic structural systems, is the flutter modal analysis in which the aerodynamic response of the structure is discretized into m degrees of freedom instead of high-order n degrees of freedom. Multi-mode flutter analysis is an approximate approach in which the system flutter oscillation mode is assumed to be the combination of a few natural modes of the target structure. Herein, m is the number of natural modes involved in the flutter oscillation. When multi-mode flutter analysis needs to be performed, the flutter equation of motion (Eq. (4)) must be presented in the modal space, in which the structural and aerodynamic matrices are generalized as: ½M  Š ¼ ½wŠT ½M Š½wŠ ð48Þ ½C  Š ¼ ½wŠT ½C Š½wŠ ð49Þ or ð41Þ ½C  Š ¼ ½wŠT ð½C s Š ½C aer ŠÞ½wŠ ¼ ½wŠT ½C s Š½wŠ ½wŠT ½C aer Š½wŠ ð50Þ and or  ( 5.1. Full-mode finite strip aeroelastic analysis Therefore, the transformed structural matrices can be obtained by: 8n o9  > ( n _ o ) € >  ½M Š ½0Š < d = ½C Š ½K Š d n o ¼ ½0Š ½IŠ > ½IŠ ½0Š : d_ > ; fdg n o X_ ¼ ½K  Š ¼ ½wŠT ½K Š½wŠ ð42Þ ð51Þ or ½K  Š ¼ ½wŠT ð½C s Š ½C aer ŠÞ½wŠ ¼ ½wŠT ½C s Š½wŠ ½wŠT ½C aer Š½wŠ ð52Þ In short form, Eq. (42) is presented by: n o ½AŠ X_ ¼ ½BŠfX g ð43Þ where ½MŠ ½0Š ½ AŠ ¼  ½0Š ½BŠ ¼  ½C Š ½K Š ½I Š ½0Š fX g ¼  z fdg ½I Š   ¼ ð44Þ  (n o) d_ fdg ð45Þ ¼   kfwg kt e fwg ð46Þ 5.3. Solution procedure In the case of full-mode flutter analysis, generalization of the property matrices is not necessary because all natural frequencies of the system must be considered, and the natural mode shapes are not required to be decoupled. For both full-mode and multi-mode flutter analysis, the eigenvalue analysis of Eq. (4) yields a set of complex conjugate pairs of eigenvalues ki ¼ li  ixi (i ¼ 1; 2; . . . ; m=2Þ. The real part li is in fact the logarithmic decrement of the structural system, while the imaginary part xi is the flutter frequency response of the structure at a certain wind speed. The real part li of at least one complex conjugate pair must be positive in order to identify the occurrence of flutter instability. In 153 H. Naderian et al. / Computers and Structures 212 (2019) 145–161 other words, when the real part of the eigenvalue becomes zero for the first time at the lowest wind speed, the critical flutter condition occurs. In performing the flutter eigenvalue analysis, the critical flutter speed U f and critical flutter frequency xf can be solved simultaneously from the real and imaginary values of the resultant eigenvalue. The range of wind speeds U min to U max in increments of DU, as well as the iteration of the oscillation frequency from the initial value of the corresponding natural mode, are important factors in the investigation of flutter susceptibility. In the case of fullmode analysis, the number of eigenvalues k is equal to 2n, while for multi-mode analysis this number is equal to 2 m. The eigenvalue equation of Eq. (4) (or Eq. (43)) is solved iteratively with the incremental value of the wind speed U. The aerodynamic derivatives are obtained from the flutter angular frequency xf , and the critical wind speed and flutter frequency are traced when the aeroelastic logarithmic decrement of the system becomes zero. The intersection of the flutter frequency curves of two or multiple modes determine the critical condition. 6. Numerical examples In order to illustrate the effectiveness of the proposed finite strip solution in the flutter analysis of long-span plates and bridge structures, in this section two examples of frequency-domain flutter analysis are presented. First, classical thin airfoil aerodynamics in the category of both multi-mode and full-mode analysis approaches is employed for the numerical flutter investigation of a thin flat shell. In the second case study, the aerodynamic performance of the Kap Shui Mun Bridge is evaluated. 6.1. Flutter of long-span flat shell including plate length L = 300 m, plate width B = 40 m, plate thickness t = 1 m, modulus of elasticity E ¼ 3:5  103 MPa, plate mass density q ¼ 2:5  103 kg=m3 , Poisson’s ratio t ¼ 1=6, and the air density qair ¼ 1:225 kg=m3 : The plate is divided into two strips of equal width in the transverse direction, and 10 equal-length sections in the longitudinal direction. The finite-strip discretization of the plate is illustrated in Fig. 6. The angle of the wind speed is selected to be zero degrees. First, a free-vibration spline finite-strip analysis is performed in order to obtain the natural frequencies and mode shapes of the plate. The developed finite strip program is also capable of 3Dplotting the deformed shape of the plate in different mode shapes. In Table 1 and Fig. 7, the first ten natural frequencies as well as the mode shapes of the simply supported flat shell are presented. From the deformed shape of the plate, it can be concluded that the second and fifth modes correspond to the first symmetric heave and torsional modes with natural frequencies of 0.0189 Hz and 0.2050 Hz respectively. As can be seen in Table 1, the targeted natural frequencies are in excellent agreement with the frequencies of the same plate, as presented in Ref. [18], with a difference of less than 5%. Theodorsen’s function [20] is used to calculate the aerodynamic flutter derivatives of a plate with a uniform thin cross-section as: C ðK w Þ ¼ F ðK w Þ þ iGðK w Þ where F ðK w Þ ¼ 1 GðK Þ ¼ A long-span plate with simply supported ends is selected as a case study to investigate the accuracy of the proposed finite strip solution. For verification purposes, the material and geometric properties are the same as those addressed in reference [18], 0:165  2 1 þ 0:0455 kw   0:165 0:0455 kw  2 1 þ 0:0455 kw 0:335  2 1 þ 0:3 kw ð54Þ   0:335 0:3 kw  2 0:3 1 þ kw ð55Þ in which kw ¼ 0:5K w . The values of Theodorsen’s functions F ðK w Þ and GðK w Þ against the kw are illustrated in Fig. 8. Using Theodorsen’s functions F ðK w Þ and GðK w Þ, the flutter derivatives Hi and Ai (i = 1, 2, 3, 4) for a plate section can be analytically presented as follows: F Kw   p 1 G F þ H2 ¼ þ Kw 2 Kw 2   p KwG H3 ¼ 2F 2 K2 w  p 4G H4 ¼ 1þ Kw 2 H1 ¼ Fig. 6. Simply supported long-span flat shell [18]. ð53Þ 2p ð56Þ Table 1 Natural frequencies of the long-span thin flat shell. Mode number 1 2 3 4 5 6 7 8 9 10 Frequency f (Hz) FSM Ref. [18] 0 0.0189 0.0755 0.1702 0.2050 0.3035 0.4150 0.4768 0.6351 0.6939 – 0.019 – – 0.194 – – – – – Angular frequency x (rad/s) Mode shape 0 0.1185 0.4744 1.0694 1.2879 1.9069 2.6076 2.9961 3.9902 4.3600 – Symmetric heave 1 Antisymmetric heave 1 Symmetric heave 2 Torsional 1 Antisymmetric heave 2 Torsional 2 Symmetric heave 3 Torsional 3 Antisymmetric heave 3 154 H. Naderian et al. / Computers and Structures 212 (2019) 145–161 Mode 6 Mode 1 Mode 7 Mode 2 Mode 8 Mode 3 Mode 9 Mode 4 Mode 5 Mode 10 Fig. 7. Finite strip mode shapes of the simply supported thin flat shell. Fig. 8. Real and imaginary parts of Theodorsen’s circularity function C ðK w Þ ¼ F ðK w Þ þ iGðK w Þ. H. Naderian et al. / Computers and Structures 212 (2019) 145–161 155 Fig. 9. Aerodynamic flutter derivatives of the thin flat shell. Fig. 10. Torsional flutter response frequencies of the simply supported long-span plate. A1 ¼ A2 ¼ A3 ¼ A4 ¼ pF 2K w  p Kw p K 2w  1 G F þ þ 8 2K w 8 ! 2 Kw F KwG þ 8 64 2 ð57Þ pG 2K w The variations of the theoretically calculated flutter derivatives Hi and Ai of the plate section against the reduced velocity are illustrated in Fig. 9. In the environment of the spline finite strip method and using an iterative frequency domain analysis, the response flutter frequencies of the plate for different wind speeds are obtained and expressed in Fig. 10. In addition, the variation of the logarithmic decrement against the wind speed is shown in Fig. 11. Only the torsional mode has been considered in this study using both multimode and full mode flutter analysis. Also, the structural damping is neglected in this example. According to Fig. 11, the logarithmic decrement becomes zero at a wind speed of 113 (m/sec), and shows the occurrence of flutter instability. The corresponding response frequency for the critical wind speed of 113 (m/sec) is 0.1055 (Hz), which represents the critical flutter frequency. In Table 2, the critical flutter wind speed 156 H. Naderian et al. / Computers and Structures 212 (2019) 145–161 aerodynamic forces allocated to each knot in a spline nodal line is slightly less than that of the present finite strip formulation. Therefore, less wind energy is input into the structure, and this initiates the occurrence of flutter at a higher speed. The other reason behind the difference between the critical wind speed of the present formulation and that obtained by Chang et al. [18] originates in the spline distribution of the aerodynamic forces in the current flutter finite strip method versus the trapezoidal distribution of the same forces employed in the work of Cheng et al. [18]. In the spline distribution, more forces are allocated to the middle knots than to the knots close to the two ends of the plate, while in the trapezoidal distribution, except for the two ends of the nodal lines, the same amount of forces are allocated to the rest of the knots. and frequency of the present plate are compared with the results reported by other researchers. The comparison between the results shows a maximum 2% and 11% difference between the critical flutter frequencies and critical flutter wind speeds respectively, which confirms the accuracy of the present finite strip formulation and the corresponding flutter analysis scheme. Also, there is almost no difference between the numerical finite strip results of the single-mode and full-mode analyses, while the computational time for these analyses was 23 s and 1620 s respectively, which shows a dramatic difference. It can be concluded that there is no need to include all the vibration modes in the flutter analysis. The reason for the higher critical wind speed in the case of the PK-F finite strip flutter analysis [18] could be the spatial distribution of the self-excited loads, in which the magnitude of the Fig. 11. Logarithmic decrement of the simply supported long-span plate. Table 2 Flutter analysis results of the simply supported long-span thin flat shell. Method Single-mode Full-mode FEM SFSM PK-F van der Put formula Reference Critical flutter frequency f (Hz) Frequency error% Critical wind speed Uf (m/sec) Wind speed error% This paper 0.1055 – 113.1 – This paper 0.1055 0 113.2 0 [21] 0.1034 +2 125.4 11 [18] 0.1045 +0.9 123.5 9 [18] – – 121.9 7 Table 3 Convergence study for critical flutter wind speed and frequency. Number of strips Number of sections 4 1 2 3 4 8 12 16 20 f (Hz) Uf (m/sec) f (Hz) Uf (m/sec) f (Hz) Uf (m/sec) f (Hz) Uf (m/sec) f (Hz) Uf (m/sec) 0.1071 0.1069 0.1069 0.1069 102.3 102.5 102.6 102.9 0.1060 0.1059 0.1059 0.1058 110.1 110.4 111.4 111.6 0.1058 0.1055 0.1055 0.1056 112.8 113.1 113.1 113.1 0.1055 0.1055 0.1055 0.1055 113.1 113.1 113.1 113.1 0.1055 0.1055 0.1055 0.1055 113.1 113.1 113.1 113.1 Fig. 12. Kap Shui Mun Bridge. H. Naderian et al. / Computers and Structures 212 (2019) 145–161 157 Fig. 13. Geometrical properties of the deck: (a) top view; (b) front view. It is worth noting that a number of parametric studies were performed in the present research, as summarized in Table 3, in order to investigate the convergence of the flutter frequency and flutter wind speed results against the spline finite strip mesh size. The numerical results of Table 3 show that the number of sections in the longitudinal direction of the strip directly influences the convergence of the results, while the number of strips has a lower impact on the accuracy and convergence rate of the solution. The reason for this is that the greater the number of knots in the longitudinal direction of the spline strip, the more uniform distribution of the aerodynamic self-excited forces along the center of elasticity of the deck. This results in faster and better convergence of the results. the natural frequencies, are obtained by the IFSM and are presented in Table 5. The finite strip results are compared with those obtained from the finite element (SAP 2000) analysis [15] as well as with the field measurements [22]. The very good agreement between the results validates the accuracy of the IFSM in frequency analysis. In Figs. 16–19, the first symmetric and antisymmetric heave and torsional mode shapes of the Kap Shui Mun Bridge are illustrated. The corresponding natural frequency of the first symmetric torsional and heave modes is 0.7526 (Hz) and 0.4250 (Hz) 6.2. Aerodynamic performance of the Kap Shui Mun bridge For the second case study, the integrated finite strip method is employed to investigate the flutter performance of the Kap Shui Mun Bridge, as illustrated in Fig. 12. The geometric and material properties of the deck are expressed in Fig. 13 and Table 4 respectively. The finite strip modeling of the west and east towers, as well as the 3D integrated finite strip modeling of the entire cable-stayed bridge system, are displayed in Figs. 14 and 15 respectively. Detailed information regarding the integrated finite-strip modeling of the Kap Shui Mun Bridge system can be found in Ref. [15]. 6.2.1. Free vibration analysis The flutter frequency analysis is performed by selecting the natural frequencies of one or a number of specific mode shapes as the start point of the response frequency. Therefore, free vibration analysis is a necessary step before flutter analysis. In this regard, a free vibration analysis is performed to obtain the natural frequencies of the Kap Shui Mun Bridge and their corresponding mode shapes. The first ten mode shapes of the bridge, including Fig. 14. Towers models of the Kap Shui Mun Bridge. Table 4 Material properties of the deck. Properties Main span Side spans Modulus of elasticity (kPa) Mass density (kg/m3) Poisson’s ratio Moment of inertia (vertical) (m4) Moment of inertia (vertical) (m4) 2.00  108 3880 0.3 191 2530 3.00  107 3630 0.2 363 5560 Fig. 15. Three-dimensional model of the Kap Shui Mun Bridge. 158 H. Naderian et al. / Computers and Structures 212 (2019) 145–161 Table 5 Modal characteristics of first ten modes of the Kap Shui Mun Bridge. Mode number Natural frequency f (Hz) Nature of mode shape IFSM (1) FEM [15] (2) Field tests [22] (3) 1 2 3 4 5 6 7 8 9 10 0.2113 0.2409 0.4250 0.5217 0.7526 0.8523 0.9305 1.0032 1.1381 1.1391 0.2061 0.2338 0.4226 0.5160 0.7179 0.8500 0.9257 1.0023 1.1048 1.1058 – – 0.39 0.49 0.83 0.66 0.90 1.07 – – Tower Tower Symmetric heave mode of the deck Lateral bending of the deck Symmetric torsional mode of the deck Antisymmetric heave mode of the deck Lateral bending of the deck Heave mode of the deck Tower Tower Fig. 16. First symmetric heave mode of the deck (0.4250 Hz). Fig. 18. First symmetric torsional mode of the deck (0.7526 Hz). Fig. 17. First antisymmetric heave mode of the deck (0.8523 Hz). respectively, while this value is 1.3419 (Hz) and 0.8523 (Hz) respectively for the first antisymmetric torsional and heave modes. 6.2.2. Kap Shui Mun bridge flutter derivatives In the design of Kap Shui Mun Bridge, because of the critical location of the bridge, typhoon conditions have been considered in which the wind loading is much stronger (might need to consider of wind speed to up 95 m/s) [23]. The Kap Shui Mun Bridge was tested for aerodynamic derivatives in order to determine the critical wind velocities for the bridge deck [24]. The aerodynamic flutter performance of Kap Shui Mun Bridge has been investigated through the experimental flutter derivatives of the bridge as reported by Jon Raggett at the West Wind Laboratory at California in collaboration with Robert Scanlan at Johns Hopkins University [25]. The static aerodynamic coefficients and the aerodynamic flutter derivatives have been derived through wind tunnel tests. The flutter derivatives, obtained experimentally, are reproduced in Figs. 20 and 21. Herein, the deck width B is 35.7 m. All the flutter derivatives are related to the final main span bridge deck configu- Fig. 19. First antisymmetric torsional mode of the deck (1.3419 Hz). ration both with and without traffic and in which the road curb is also considered. The vertical incidence angle of the wind, sometimes called the angle of the wind attack, in the present study is -2.5 degrees. 6.2.3. Flutter analysis of the Kap Shui Mun bridge After obtaining all structural and aerodynamic properties of the the Kap Shui Mun Bridge by using the integrated finite strip method, the flutter analysis can be performed. The damping ratio of the bridge system is chosen to be 0.02. In the previous example, it was shown that full-mode flutter analysis takes much more storage and time comparing with the multi-mode flutter analysis, while the results of both techniques are almost the same. The full-mode flutter analysis for a massive structure like the Kap Shui Mun Bridge with thousands of degrees of freedom will take a mas- H. Naderian et al. / Computers and Structures 212 (2019) 145–161 159 Fig. 20. Flutter derivatives A1 to A3 of the Kap Shui Mun Bridge against reduced velocity. sive amount of storage and time. Therefore, multi-mode flutter analysis is applied in this example for aeroelastic analysis of the Kap Shui Mun Bridge. The first symmetric torsional mode of the deck is taken as the critical mode shape for the flutter. The integrated finite-strip flutter eigenvalue analysis results are plotted in Figs. 22 and 23, where the flutter frequencies as well as the logarithmic decrements of the Kap Shui Mun Bridge for different wind speeds are calculated and displayed. The results show that both the flutter frequency and logarithmic decrement decrease when the wind speed increases, which could be a sign of flutter instability. However, the logarithmic decrement never becomes zero, which means that the flutter instability does not happen within the standard values for wind speeds and even in typhoon conditions, for which the critical wind speed could be 95 m/sec. Therefore, this proves that the aerodynamic flutter design of the Kap Shui Mun Bridge meets the requirements for aeroelastic stability. 6.3. Efficiency analysis Using a personal computer with an Intel Core2 Duo CPU (1.66 GHz) and 3 GB physical memory, the computer times required by IFSM and FEM for two test cases with similar grid conditions were compared. For the first case study, the computational time in IFSM and FEM was 4000 s and 4400 s, respectively, and in the second case study it was 12,300 s and 13,600 s, respectively, which shows superiority of IFSM in terms of computational cost. In other words, around 10% percent decrease in computation time can be witnessed in finite strip method compared to finite element analysis. Due to the semi-analytical properties in the longitudinal direction, the number of sections for each strip in the IFSM model could be significantly reduced without losing accuracy. Consequently, the time for computation could also be reduced, and greater efficiency is achieved. It must be noted that this case study 160 H. Naderian et al. / Computers and Structures 212 (2019) 145–161 Fig. 21. Flutter derivatives H1 to H3 of the Kap Shui Mun Bridge against reduced velocity. Fig. 22. Torsional flutter response frequencies of the Kap Shui Mun Bridge. was performed on a simplified model. In the case a complicated long-span cable-stayed bridge model, the efficiency of the proposed integrated finite strip solution will be more appreciated due to the need for fewer iterations and less degrees of freedom in comparison with FE analysis [15]. For instance, when a continuous real time analysis is required in which a lot of iterative analyses are performed, the advantages of the proposed methodology in time saving will be more illustrated. H. Naderian et al. / Computers and Structures 212 (2019) 145–161 161 Fig. 23. Logarithmic decrement of torsional flutter of the Kap Shui Mun Bridge. 7. Conclusions The integrated finite strip method was extended to the area of wind engineering, where aerodynamic flutter analysis can be performed in a minimal computational time, providing high accuracy and efficiency. Due to the semi-analytical nature and narrow bandwidth of the dynamic properties of the system, the results converged rapidly, with a few degrees of freedom. The aeroelastic damping and stiffness matrices were derived by using the integrated finite strip method, and the total potential energy of a shell strip was obtained from the algebraic summation of the in-plane and out-of-plane deformations. In the current finite strip solution, the boundary conditions were modelled by using a simple technique based on a physical coordinate system. Also, the flutter problem was solved through a very straightforward eigenvalue frequency analysis. Both multi-mode and full-mode flutter analyses approaches have been outlined in the IFSM environment in order to evaluate the critical flutter speed as well as the critical flutter frequency of the structural system. The numerical results show that performing the full-mode flutter analysis is much more time-consuming than the multi-mode analysis, while the resultant critical flutter frequencies and corresponding flutter wind speeds do not change significantly. Therefore, it is suggested to use the multi-mode technique for flutter analysis. Also, the results show that the uniform distribution of self-excited flutter forces on the centre of elasticity of a cross-section is reasonably good when compared with the spatial distribution of aerodynamic forces on the entire surface of the deck. Last but not least, the method is capable of performing a three-dimensional flutter analysis. The successful extension of the finite strip method to aerodynamic flutter analysis can open a new door for further development of this methodology in different areas of wind engineering, such as buffeting analysis of a bridge, wind-vehicle-structure (WVS) interactions, and the aerodynamics of composite hybrid FRP bridge structures. Acknowledgements The authors wish to gratefully acknowledge the financial support of the Canada NSERC Industrial Postgraduate Scholarship. References [1] Ge YJ, Tanaka H. Aerodynamic flutter analysis of cable-supported bridges by multi-mode and full-mode approaches. J Wind Eng Ind Aerodyn 2000;86:123–53. [2] Ding QS, Chen AR, Xiang HF. Coupled flutter analysis of long-span bridges by multimode and full order approaches. J Wind Eng Ind Aerodyn 2002;90:1981–93. [3] Huang Ming-Hui, Lin Yuh-Yi, Weng Ming-Xi. Flutter and buffeting analysis of bridges subjected to skew wind. J Appl Sci Eng 2012;15(4):401E413. 401. [4] Chad Gibbs S, Wang Ivan, Dowell Earl. Theory and experiment for flutter of a rectangular plate with a fixed leading edge in three-dimensional axial flow. J Fluids Struct 2012;34(2012):68–83. [5] Wu Teng, Kareem Ahsan, Ge Yaojun. Linear and nonlinear aeroelastic analysis frameworks for cable-supported bridges. Nonlinear Dyn 2013;74(3):487–516. [6] Kreis ES, Andr´e JC. A numerical inquiry into the flutter phenomenon in longspan bridges. Latin Am J Solids Struct 2005;2(2005):321–37. [7] Matsumoto Masaru, Matsumiya Hisato, Fujiwara Shinya, Ito Yasuaki. New Consideration on flutter properties basing on SBS -Fundamental flutter mode, similar Selberg’s formula, torsional divergence instability, and new coupled flutter phenomena affected by structural coupling. BBAA VI international colloquium on: bluff bodies aerodynamics & applications Milano, Italy, July, 20-24 2008, 2008. [8] Phan Duc-Huynh, Nguyen Ngoc-Trung. Flutter and buffeting control of longspan suspension bridge by passive flaps: experiment and numerical simulation. Int’l J Aeronaut Space Sci 2013;14(1):46–57. https://doi.org/ 10.5139/IJASS.2013.14.1.46. [9] Canor Thomas, Caracoglia Luca, Denoël Vincent. Application of random eigenvalue analysis to assess bridge flutter probability. J Wind Eng Ind Aerodyn 2015;140:79–86. [10] Xie Gui-hua, Yin Jie, Liu Rong-gui, Chen Bei, Cai Dong-sheng. Experimental and numerical investigation on the static and dynamic behaviors of cable-stayed bridges with CFRP cables. Compos Part B: Eng 2017;111:235–42. [11] Yang Yongxin, Zhou Rui, Ge Yaojun, Zou Xiaojie, Zhang Lihai. Flutter characteristics of twin-box girder bridges with vertical central stabilizers. Eng Struct 2017;133:33–48. [12] Lee N, Lee H, Baek C, Lee S. Aeroelastic analysis of bridge deck flutter with modified implicit coupling method. J Wind Eng Ind Aerodyn 2016;155:11–22. [13] Cheung MS, Li W, Chidiac SE. Finite strip analysis of bridges. 1st ed. London: Spon; 1996. [14] Cheung YK, Tham LG. The finite strip method. USA: CRC Press; 1998. [15] Naderian Hamidreza, Cheung Moe MS, Shen Zhenyuan, Dragomirescu Elena. Integrated Finite Strip Analysis of long-span cable-stayed bridges. Comput Struct 2015;158:82–97. [16] Naderian Hamidreza, Cheung Moe MS, Shen Zhenyuan, Dragomirescu Elena. Seismic analysis of long-span cable-stayed bridges by an integrated finite strip method. ASCE J Bridge Eng 2015. ISSN 1084-0702. [17] Cheung YK, Fan SC, Wu CQ. Spline finite strip in Structure Analysis. In: Proceedings, the international conference on finite element method, Shanghai. p. 704–9. [18] Lau DT, Cheung MS, Cheng SH. 3D flutter analysis of bridges by spline finite strip method. J Struct Eng 2000:1246–54. https://doi.org/10.1061/(ASCE) 0733-9445(2000)126:10(1246). [19] Scanlan RH, Lin WH. Effects of turbulence on bridge flutter derivatives. J Engrg Mech, ASCE 1978;104(4):719–33. [20] Simiu E, Scanlan RH. Wind effects on structures. 3rd ed. New York: Wiley; 1996. [21] Cheng SH. ‘‘The 3D finite element flutter analysis of long-span bridge”, Master thesis. Shanghai, China: Tongji University; 1994. [22] Lau CK, Mak WP, Wong KY, Man KL, Chan WY, Wong KF. Structural performance measurement and design parameter validation for Kap Shui Mun Cable-Stayed Bridge. Adv Steel Struct 1999;2(1). [23] Yung Benedict Kai Kwan. A critical analysis of the Kap Shui Mun Bridge. Proceedings of bridge engineering conference, 23 April 2008, England, 2008. [24] Scanlan RH, Stroh SL, Raggett JD. Methods of wind response investigation employed for the Kap Shui Mun Bridge. J Wind Eng Ind Aerodyn 1995:1–11. [25] Raggett Jon, Scanlan Robert. Flutter derivatives of Kap Shui Mun Bridges using wind tunnel tests. West Wind Laboratory, California & Department of Civil Engineering, Johns Hopkins University; 1997.