Innovative approach with harmonic modes and finite element modeling for the nonlinear dynamic analysis of a suspension bridge

The nonlinear characteristics of a general multiple degree-of-freedom suspension bridge under harmonic excitation are studied with an innovative approach consisting of the harmonic modes, the finite element modeling and the incremental harmonic balance method. The geometric nonlinearity induced by large-amplitude vibration is considered. The harmonic mode and mean kinetic energy of the system are proposed to describe the nonlinear responses. The harmonic mode can be expressed as a linear combination of different linearized system modes. The resonance responses of the structure are characterized with the participation ratio of these system modes. Its composition under harmonic excitations can easily be described in terms of the different harmonic mode shapes and then the system modes. The proposed method is further applied to analyze, in detail, the properties of the harmonic modes and the evolution of the steady-state nonlinear responses with variations in the excitation frequency.

Abbreviations a 1 * a 12 Physical coefficients a ik Modal phase of the kth system mode a jk, , b jk Coefficients of Fourier series a k Scaling factor of the kth system mode b ik Modal amplitude of the kth system mode A, B Node Lumped mass matrix of the link element P ik Parameter related to contribution of different system modes in each harmonic mode q Dynamic cement response q 0 .
Dynamic displacement cable related to x 0 K Stiffness displacement response R Error vector t Time u 1, u 2 Horizontal displacements of node v 1, v 2 Vertical displacements of node V(t) Vibration reconstructed by system modes W Real vibration shapes Y Vibration shape fitted by system mode shapes a, b Coefficients of Rayleigh damping DðÁÞ Interval of physical variable x Circular frequency of excitation x 0 Frequency close to x s Non-dimensional time u i Mode shape vector of the ith mode

Introduction
The dynamic analysis for of a large-scale structure is usually difficult due to the complexity of system. Discretization is always necessary before the model can be applied for the analysis. There are basically two kinds of discretization approaches, namely modal discretization and space discretization. The first approach treats the structure as a distributed-parameter system (continuum model). The system vibration is decoupled into time-invariant space functions (mode shapes) and time-dependent scaling coefficients. The system dynamic response can then be represented usually by the responses of the few lower-order modes. The modal information can be easily investigated by this model [1][2][3][4]. More recently, Gwon and Choi [5][6][7] conducted a series of studies with the continuum model of a three-span suspension bridge with special focuses on the effects of hanger extensibility, boundary conditions of the girder and the configuration of the main cables.
The nonlinear vibration of suspension bridge has also been widely studied based on the continuum model. Abdel-Ghaffar and Rubin [8,9] presented the nonlinear free vertical-torsional coupled vibrations of suspension bridges. Abdel-Rohman and Nayfeh [10] adopted a nonlinear continuum model to investigate the effectiveness of an active damper. Lazer and McKenna [11] developed a continuum suspension bridge model governed by the coupled nonlinear wave and beam equations. Most of these studies adopted the multiple scale method. However, it is limited up to considering the coupling effects of two or three modes [8][9][10][11]. They are, therefore, more popularly used for preliminary studies.
The second approach discretizes the structure in space with a number of finite elements. Finite element method can be more accurately and flexibly model the target structure. Thus, it has been extensively used in analyzing the dynamics of various suspension bridge [12][13][14]. West et al. [12] studied the effects of sag of cable, moment of inertia of deck, and cable support stiffness on the vertical vibration characteristics of a single-span suspension bridge. Finite element modeling has also been widely adopted in recent years for the analysis of wind-induced structural response [13,14]. However, finite element modeling approach of suspension bridge can seldom be used for systematic nonlinear dynamic analysis. That is because such models usually have elements, which means large number of degrees of freedom (DoFs), and it is almost impossible to apply multiple scale method in such circumstances.
The incremental harmonic balance (IHB) method [15,16] has been used to obtain solution for a multi-DoF system. It has since been improved to analyze the quasi-periodic motion of nonlinear systems [17][18][19]. Hui, et al. have investigated the nonlinear behavior of a 6-DoF sectional model of a cable-stayed bridge [20][21][22]. Although the IHB method has the capacity of handling the multi-DoF structure, investigation of nonlinear dynamics of structure with a large number of DoFs has not been reported. Besides the algorithm of solution finding, analyzing the obtained results of the complex structure also requires more effective approach, as the classical examination of response by checking some specific DoFs may not be sufficient to interpret the system response.
This study addresses the limitations of existing methods in the nonlinear dynamic analysis of a structure with a large number of DoFs. A harmonically excited finite element model of suspension bridge, which is composed of 52-DoFs, is studied using the IHB method. For better result interpretation and understanding, a new approach is proposed to describe the evolution of the nonlinear dynamic responses with changing excitation frequency and amplitude.
The structure of this paper is organized as follows: Sect. 2 derives the suspension model, and necessary validations are carried out. Section 3 introduces the IHB method which is used to calculate the structure response. Basic response property is examined in Sect. 4. In Sect. 5, the harmonic mode, the mean kinetic energy and system modes are proposed for interpreting the bridge response, and the approach is illustrated in the analysis. Conclusions are given in Sect. 6.

Bridge model and validation
The two-dimensional (planar) finite element model of a steel suspension bridge with 10 hangers is shown in Fig. 1. It is 11 m long with 1/10 sag-to-span ratio. The 10 hangers are uniformly distributed along the span with the shortest one 0.2 m long. The elastic modulus of material and mass density are 210 Gpa and 7.85 9 10 3 kg/m 3 , respectively. The cross-sectional areas of the main cable, hanger, and the bridge girder are, respectively, 0.005 m 2 , 0.001 m 2 and 0.01 m 2 . The moment of inertia of the bridge girder about its horizontal axis is 8.33 9 10 -8 m 4 , and the main cable and hangers are assumed flexible. As the model is a planar model, only the deformations and forces within the plane will be considered in this study.
The cable between two adjacent hangers or end support is modeled with a straight link element, ignoring its bending stiffness [17,23,24]. The bridge girder is modeled with Euler-Bernoulli beam elements, and hanger is also modeled as a link element. The deck has hinge supports at both ends. Considering that the contribution of the nonlinear stiffness of bridge deck to the response is not significant comparing with the cables with nonlinear stiffness, the stiffness of the bridge deck is assumed as linear in this study, as its contribution to the overall stiffness of bridge is much smaller than the cable. Diagonal mass matrix is adopted for convenience for both the beam and link elements. The nonlinearity of the system comes from the geometric nonlinearity of the main cable and hangers affecting the stiffness of the suspension bridge.
As cable is the most critical component of the bridge, it is modeled at first in this study. By setting the sag-to-span ratio of the cable and considering the static balance function of cable under self-weight of cable and deck, the initial internal force of each element and the shape of cable at static equilibrium position (within the x-z plane) can be determined.
Since the cable between two adjacent hangers or end support is modeled with a straight link element, the in-plane (x-z plane) stiffness matrices of a link element with an initial internal axial force F c as shown in Fig. 2  where k = EA c /l 0 , E is the elastic modulus of the material, A c is the sectional area of cable, and l 0 is the original length of the element with no axial force. The displacements can be expanded into Taylor series up to the third order in this paper, as third-order nonlinearity has been found significantly contributing to an accurate dynamic analysis of such cablesupported structures [23]. Equations (1) and (2) can be rewritten as: . and the linear and nonlinear stiffness matrices of the cable element can be written as: where matrices K c1 , K c2 and K c3 represent, respectively, the linear, second-and third-order nonlinear stiffness matrices of cable, respectively. The stiffness matrices of the hanger element can also be obtained similarly by setting d = 0 (Fig. 2). Besides the stiffness matrix, the lumped mass matrix of the link element is shown as follows: where m A and m B indicate the masses at points A and B, respectively. Usually the entire mass of the element is divided equally between m A and m B . After obtaining all the linear and nonlinear stiffness matrices and mass matrix of every element, the corresponding matrices of the entire cable can be easily obtained by assembling those of individual cable elements in the finite element modeling (FEM) method. The nonlinear stiffness of the cable is then validated by comparing the deformation of cable calculated based on the nonlinear stiffness and that obtained based on shape-finding algorithm (adopting the static equilibrium function as mentioned above), which is the exact solution. A uniform static force is applied on every node of the cable as shown in Fig. 3a. The variation of mid-span deformation along with the force is shown in Fig. 3b. The positive value of force indicates the downward direction. It can be clearly checked from the exact solution that cable becomes more and more stiff with the increase in force in the downward direction, but the much more severe softening phenomenon can be observed in the opposite direction, when the force is applied upward. The asymmetric nonlinear stiffness of cable arises from the curved shape of cable. When the cable deforms downward under uniform force, the axial force increases, and the stiffness increases accordingly; when it deforms upwards, the axial force decreases, and the cable will lose its stiffness gradually. At the same time, it can be found that when up to the thirdorder nonlinear terms, the calculated formation is almost identical to the exact solution. However, if only up to the second-order nonlinear terms are adopted, the big differences can be easily observed at relatively larger deformation. Such results indicate that the thirdorder nonlinear terms of the stiffness are essential for large deformation analysis.
The finite element model of the bridge model is then established, and the first four modal frequencies of the linearized system are, respectively, 7.73, 12.40, 19.91 and 28.02 rad/s. The first four mode shapes of the structure are shown in Fig. 4. Such results are identical to that obtained from the commercial finite element analyzing software, indicating the reliability of the model. The first and third modes are noted antisymmetric, while the second and fourth modes compose of symmetric vibration. Rayleigh damping C = aM ? bK 1 is adopted with a = 0.003, b = 0.0012. M and K 1 are the system mass and linear stiffness matrices, respectively. The damping ratios of the first four modes are therefore obtained as 0.48, 0.76, 1.20 and 1.70%, respectively, which fall in the usual range with steel bridge.

Incremental harmonic balance method
Considering that the incremental harmonic balance method has been proved to be reliable and accurate in seeking the response solution of nonlinear systems so it is adopted in this study [20][21][22]. Details of the incremental harmonic balance (IHB) method can be found in [16]. Since the third-order nonlinear stiffness of cable is considered, the dynamic equation of motion for the damped system can be written as: where F is a vector of the amplitude and location of external harmonic excitations. Variable x is the circular frequency of excitation.  x 0 , which is close to x. Thus, let q 0 be the dynamic displacement responses corresponding to x 0 . A frequency and displacement close to this system state can be expressed as: Equation (7) is substituted into Eq. (6) to get where R is a vector of error and s = xt, and the terms with the prime denote the differentiation with respect to s.
Various studies have proved that the dynamic response of nonlinear system under harmonic excitation can be expressed in the form of a Fourier series, whose frequency components are integral multiples of the excitation frequency [16,20,25,26]. In this study, five frequency components are adopted in this study: The component with the excitation frequency is denoted as the first harmonic term. The components whose frequency is two to four times that of the applied excitation are named as the second to fourth harmonic terms, respectively, while the constant term is referred to as the stationary term in this study.
Equation (12) is then substituted into Eqs. (8) and (9). Galerkin process is then applied to one cycle with s between null to 2p such that the incremental linear equation is as follows: where Solution of Eq. (13) can be numerically obtained with the Newton-Raphson iterative method or the arclength iterative method. It should be noted that the initial set of solutions of the tracking process needs to be determined carefully with trial-and-error method.

Effect of nonlinear stiffness of hangers
In most nonlinear dynamic analysis of suspension bridges, which are based on modal-based discretization, the hangers are always modeled with membrane uniformly distributed along the cable connecting the bridge deck [27,28]. The deformation of membrane during vibration is ignored, which means its effect on the dynamic behavior of the bridge structure is ignored.
In this study, the large deformation of hangers is considered by taking into account their nonlinear stiffness. Such effect is firstly examined based on the deformation under static loads. Anti-symmetric loads are applied to the bridge deck as shown in Fig. 5a. The deformation of the bridge is shown in Fig. 5b when the force amplitude at each node is about 557 N. The deformation of Node 4 with or without considering the nonlinear stiffness of hangers is compared in Fig. 6a with varying load amplitude. It can be clearly noted that when the nonlinear stiffness of hanger is considered, the system shows clear softening feature with increasing force. It is also clear that the inclination of hangers is responsible for this phenomenon, because there will not be any nonlinear stiffness if the hangers are kept vertical. Figure 6b shows the inclination of hanger connecting Nodes 4 and 15 with varying force amplitude. These results indicate that when the bridge girder's largest deformation large (larger than 0.04 m in this case), the inclination of hanger will have significant effects on the responses of the bridge.
To verify this point, four harmonic forces with the same amplitude (180 N) and frequency are applied simultaneously and vertically on the bridge deck at Nodes 3, 4, 9, and 10 to excite the first system mode,    Fig. 5. Two models are prepared with one considering the nonlinear stiffness of hangers but not the other. Figure 7 shows the frequency-response curves of the first two harmonic terms obtained at Node 4 from the two models under the same set of excitations. Differences in the responses of the two models can be observed when the response amplitude of the first harmonic term is over 0.03 m. It can be found that when the nonlinear stiffness of hangers is considered, the system shows stronger softening feature. At the same time, the peak response amplitude of the first harmonic term is much higher when the nonlinear stiffness is not considered. The pattern of the two curves in the first two harmonic terms also shows significant differences. These may indicate that the inclusion of nonlinear stiffness of the hangers in the nonlinear analysis of the system is necessary when large vibration of the deck is induced. Therefore, both the nonlinear stiffness of the main cable and hangers will be taken into account in the following studies for an accurate nonlinear response analysis of the bridge deck.

Vertical responses
Two identical 225 N harmonic excitations are applied vertically at Nodes 6 and 7 of the bridge structure to excite the second system mode. The excitation frequency varies between 8.0 and 18.0 rad/s with increments determined by the arc-length method. This range covers the second modal frequency of the system, and the system responses are calculated with the IHB method. The characteristics of the nonlinear responses are reported in the following sections. Figure 8 shows the frequency-response curves at the vertical DoFs of Nodes 3 and 6 located at 1/4 span and mid-span of the model, respectively. The two curves exhibit different patterns with increasing excitation frequency. The plots of stationary, first harmonic and third harmonic components feature a dominant peak at 12.18 rad/s indicating primary resonance of the second mode with larger response at mid-span. Another dominant peak can be observed at 14.00 rad/s in the second harmonic component at Node 6 but not at Node 3. The responses around 12.00 rad/s at Node 3 are much larger than those from Node 6. Moreover, there are two peaks at 11.66 rad/s (a) Stationary term (b) 1st harmonic term (c) 2nd harmonic term (d) 3rd harmonic term (e) 4th harmonic term The peaks in different harmonic components indicate some kinds of resonance phenomena, but many of them cannot be clearly recognized. The first harmonic component dominates the response, while other components are much smaller. Additionally, different positions may have different response characteristics as shown making it difficult to estimate the contribution of different frequency components to the overall responses. The above results suggest that, though any existing analysis technique can be used to find the nonlinear responses, it is difficult to explain the obtained results with reasons. An analysis strategy using existing techniques needs to be developed to analyze and interpret in simple terms the nonlinear dynamic response. Figure 8 shows that the response of a nonlinear multi-DoF system composes several harmonic components. The vibration of the whole system cannot be depicted fully by the responses from a few DoFs. Therefore, the harmonic mode for describing the nonlinear characteristics of the responses is proposed to address this issue. It is defined as the shape of vibration of the system at every DoF of the structure for a specific frequency component. Since each harmonic term includes cosine and sine components as discussed in Eq. (10), each harmonic mode would include two components vibrating at the same frequency but with p/2 phase difference. There is also one stationary mode (SM) corresponding to the shape formed by the stationary terms in the responses at each DoF. The harmonic mode shapes along the bridge deck are plotted in Table 1.

Harmonic mode
According to the results shown in Table 1, the properties of the harmonic mode are summarized as follows: 1. A harmonic mode is not necessarily identical with the system mode shape. They may be different as shown in Fig. 4 and Table 1; 2. Other than the phase difference, the mode shapes of the cosine and sine components of the harmonic modes are usually different; 3. The mode shapes of different harmonic modes may be quite different when the system vibrates in steady state. 4. The harmonic mode shape varies with the excitation frequency. The lower modes are less affected, but the higher modes are relatively more sensitive to the excitation frequency. The dashed line represents the mode shape from the cosine terms, and the solid line denotes the mode shape from the sine terms. All of the harmonic mode shapes in this paper will be similarly represented

Quantitative evaluation of response
The harmonic mode shapes in Table 1 show that the response at a DoF as shown in Fig. 8 can hardly provide a proper description on the participation of different harmonic terms in the overall response, because the vibration at different DoFs is very different even with the same excitation frequency. Thus, the vibration-induced system kinetic energy is calculated as an index on the participation of different harmonic modes to the overall response. It is calculated as summation of the kinetic energy from each DoF of the system as: where n is the total number of DoFs of the system. According to Eq. (10), the velocity at each DoF consists of four frequency components, which are the integer times of the excitation frequency. The kinetic energy at the ith DoF of the system can then be expressed as: There are square terms of each frequency component as well as the products of different frequency components (cross-frequency terms) in Eq. (15).
According to the properties of trigonometric function, the mean values of the cross-frequency terms over time will be equal to null, and the mean of the square terms is 0.5. Therefore, the mean kinetic energy of response at a single DoF can be written as: The mean kinetic energy (MKE) of the whole system can then be expressed as: Equation (17) gives the MKE as the summation of MKE of all frequency components in the response. This provides a means to quantify the contribution of each frequency component to the overall response. It should be noted that the kinetic energy induced by the stationary term is null. Figure 9 shows the frequency-MKE curves of different harmonic modes. An inspection of Fig. 4 and Table 1 shows the significance of each harmonic mode response without the need to check the response amplitude at different DoFs.

The orthogonal basis of vibration
Since the system modes of a structure are orthogonal to each other, any vibration may be expressed as a linear combination of the system modes. Therefore, the system modes will be adopted in this study for the decomposition of the nonlinear responses. The mode shape vectors are orthogonal to each other, satisfying the following requirement: (a) The 1st harmonic mode (b) The other harmonic modes where u i (i = 1,..,m) is the mode shape vector of the ith system mode, and m is the total number of the participating mode shapes. Table 2 shows the mode shapes of the lowest 9 modes and the 12th mode. It can be seen that with the increase in order, the local behavior becomes increasingly prominent. All the mode shapes are normalized to have a unit largest displacement for evaluating the participation of different modes. The sine and cosine components of a harmonic mode can then be expressed as: where a k (k = 1,..,m) is the scaling factor of the kth system mode. It is found that 12 th mode is a critical mode for describing the system vibration especially the symmetric vibration. This might due to the fact that the 12th mode reflects the most basic symmetric pattern of the model, although it has a relative high modal frequency. The harmonic mode shapes at x = 12.18 rad/s shown in Table 1 are adopted for examining the proposed approach, and the error between the real vibration shape (W) and that from curve fitting from several selected system mode shapes (Y). The error is defined in Eq. (20). Figure 10 shows the error of the fitted curves for different harmonic modes. It can be noted that when a smaller number of system modes are involved, i.e., mode 1:2, and 12, the errors are quite large especially in the higher harmonic modes. The error decreases significantly with an increase of number of system modes in the curve fitting. When nine modes are adopted, the errors of all fitted curves become numerically small close to null. This observation indicates that many system modes might be excited simultaneously in different harmonic modes (vibration frequencies). However, higher system modes are more likely to be excited with increasing order of the harmonic mode.
The cosine and sine components of the harmonic mode shape are noted to compose of the same system modes but with different scaling factors obtained from Eq. (20). It is then easy to combine the two components through a simple manipulation of trigonometric functions as: ¼ À2:2 Â 10 À3 u 12 cosð2xt þ 3:03Þ þ 9:17 Â 10 À4 u 2 cosð2xt À 0:44Þ þ 1:5 Â 10 À3 u 4 cosð2xt À 2:94Þ: The above discussions show that a single harmonic mode can be viewed as similar to the modal vibration of a linear system, except that the latter consists of vibration at a single frequency x, whereas the former has vibration components vibrating at several frequencies including the excitation frequency and its multiples, as x, 2x, 3x, … and with a static component represented by the SM.
Results shown in Fig. 11 indicate that the first eight system modes plus the 12 th mode are required for an accurate representation of the harmonic modes in the current study.
Hence, the displacement response of the ith harmonic mode can be generally expressed in the following form as: where b ik (k = 0, 1,…, 8,12) and a ik (k = 0, 1,…, 8,12) are the modal amplitude and phase of the kth system mode. Coefficient b ik is a scaling factor without unit, and b ik C 0 for any i and k.

Contribution of system modes
A change of excitation frequency x will also lead to a change in the modal amplitudes of different system modes and the harmonic mode shape. Figure 12 shows the modal amplitudes of different harmonic modes versus x. The second system mode is noted dominating the first harmonic mode, with the maximum amplitude at around 0.05. The fourth system mode dominates the fourth harmonic mode, especially when x = 11.66 and 12.49 rad/s. The contribution of different system modes in each harmonic mode can be quantitatively evaluated as: . . .; 6; 12 ð Þ : ð23Þ Figure 13 shows the contribution of different system modes to a harmonic mode. One system mode is noted participating in different harmonic modes, suggesting that it can contribute to vibration at different frequencies simultaneously under a monofrequency harmonic excitation. This marks the big difference between the vibration of a nonlinear system and linear system. The second system mode contributes 98.9% of the first harmonic mode when x = 12.18 rad/s as shown indicating a primary resonance of the second mode. The second harmonic mode shows an indistinct peak at x = 14.0 rad/s in Fig. 9b. This mode is contributed up to 69.1% by the fourth system mode. Since the fourth modal frequency of the linearized system is 28.02 rad/s, it is reasonable to claim that super-harmonic resonance of the fourth system mode may have been induced when x = 14.0 rad/s. The modal participation ratios corresponding to the two peaks at 11.66 and 12.49 rad/s in the fourth harmonic mode in Fig. 9b are also calculated. The fourth system mode dominates the fourth harmonic mode at these two excitation frequencies, although the vibrating frequencies of the fourth Fig. 11 Comparison of harmonic modes and a combination of system modes harmonic mode in these two cases do not have any distinct relation with the fourth modal frequency of the linearized system. This may still suggest some kind of resonance related to the fourth system mode.

Conclusions
A new strategy is proposed to analysis and interpret the nonlinear responses of a multi-DoF system. The incremental harmonic balanced method and finite element modeling techniques are adopted.
The harmonic mode and mean kinetic energy of the system are proposed to represent the nonlinear response. The harmonic mode shapes are represented as linear combinations of the system mode shapes similar to the modal responses of a linear system. Different resonance responses of a nonlinear system can be categorized in terms of the contribution of each system mode to the harmonic mode. Other features of the nonlinear system can then be conveniently labeled and referenced. This general form of the harmonic mode and mean kinetic energy may be applicable for the interpretation of nonlinear vibration of a general multi-DoF structure.
The proposed approach is illustrated with the analysis of a 52-DoFs finite element model of suspension bridge. The geometric nonlinear stiffnesses of cable and hanger induced by large vibration amplitude are considered. The interaction of the first few modal resonance responses of the nonlinear system is investigated. In checking the system responses within a wide frequency range, a very simple excitation is found to induce many different types of resonance responses in a high-dimensional system. Although some of them can be categorized with the help of the proposed approach, many of them are yet to be studied further to reveal the true behavior of a realistic suspension bridge structure. For instance, although the response in this study can be seen as combination of several system modes, it remains unclear about why the only selected modes contribute to a specific case. In addition, real engineering structure might be more complicated, the present  IHB method might be difficult to be applied for analysis, due to its low efficiency.