Chapter 6 Seismic analysis Contents 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 6.2 Mathematical approach . . . . . . . . . . . . . . . . . . . . . . 130 6.3 Elastic analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 130 6.3.1 Direct Response History Analysis: DRHA . . . . . . . . . . 131 6.3.2 Modal Response History Analysis: MRHA . . . . . . . . . . 133 6.3.3 Modal Response Spectrum Analysis: MRSA . . . . . . . . . 136 6.3.4 Previous studies . . . . . . . . . . . . . . . . . . . . . . . . . 140 6.3.5 Comparison of the results . . . . . . . . . . . . . . . . . . . . 141 6.4 Inelastic analysis . . . . . . . . . . . . . . . . . . . . . . . . . 150 6.4.1 Modal Pushover Analysis: MPA . . . . . . . . . . . . . . . . 151 6.4.2 Extended Modal Pushover Analysis: EMPA . . . . . . . . . 158 6.4.3 Coupled Nonlinear Static Pushover: CNSP. A new method . 161 6.4.4 Discussion of the results . . . . . . . . . . . . . . . . . . . . 166 6.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177 6.1 Introduction Several seismic analysis procedures are currently available for the designer, varying their complexity, accuracy and computational cost. Addressing the magnitude of the errors which may be introduced by the selected integration scheme, and where they come from, is paramount in order to ascertain the degree of reliability of the results. Thisvaluableinformationwillassistintheselectionoftheoptimumanalysis strategy, based on the seismic risk and the importance of the structure. This chapter includes a thorough review of the current practice in analysis methodologies, applied both to the elastic and inelastic seismic responses of the studied cable-stayed bridges, establishing the level of errors introduced in each case and proposing a new pushover procedure, appropriate for the nonlinear seismic be- haviour of three-directionally excited structures with strong modal couplings. The response of a structure under ground motion excitations is presented from a mathematical point of view in section 6.2. The rest of the chapter deals with dif- ferent procedures to face this problem and their inherent errors or implementation 130 Chapter 6. Seismic analysis particularities; section6.3simpli(cid:28)esthepicturebyconsideringtheelasticproperties of the materials throughout the structural response, detailing both spectral, modal dynamics and direct response history analysis; section 6.4 addresses the full non- linear problem, contrasting the results given by normative and advanced pushover strategies with the reference solution obtained by means of nonlinear response his- tory analysis, and proposing new methodologies adapted for cable-stayed bridges. The chapter is closed with the main conclusions in section 6.5. Alltheproceduresemployedinthissectionstartfromthedeformedcon(cid:28)guration due to the self-weight of the structure (see chapter 3 for more details). 6.2 Mathematical approach: the system of dynamics The coupled system of di(cid:27)erential equations, which governs the dynamic response of a structure subjected to an earthquake excitation, was already introduced in section 2.2.3. Here, it is extended for convenience to a 6N-degree of freedom model; composedofN nodeswith6degreesoffreedomeach(threedisplacementsandthree rotations). The system of dynamics is obtained by establishing the d’Alembert balance of forces (including inertia, damping and deformation-related forces): mu¨+cu˙ +f (u,u˙) = mιu¨ (6.1) S g − Where u(t) [6N 1] is the relative displacement vector, m and c [6N 6N] are × × respectively the mass and damping matrices of the structure, whereas f relates the S forceanddisplacementvectors. Finally, ι[6N 3]isthein(cid:29)uencematrixconnecting × the degrees of freedom of the structure and imposed accelerogram directions u¨ (t), g which generally, and neglecting the foundation rotations (section 5.2), is a vector with three components u¨T(t) = (u¨X,u¨Y ,u¨Z), where u¨j is the ground acceleration g g g g g in j-direction (j = X,Y,Z, see (cid:28)gure 5.2 for the description of the axis system and accelerogram components). The multiple support seismic excitation is not consid- ered(section3.6)andhencetheaccelerogramsareassumedequalalongthesupports of the bridge. The e(cid:27)ect of seismic devices is not taken into account in expression 6.1 nor in the rest of the present chapter, the interested reader is referred to chapter 9. 6.3 Elastic analysis If it is assumed that the structure behaves in the linear elastic range during the seismic response; f = ku in expression (6.1), being k the elastic sti(cid:27)ness matrix of S the structure, leading to the following traditional form of the system of dynamics: mu¨+cu˙ +ku = mιu¨ (6.2) g − Since the simpli(cid:28)ed deck-tower connection prevents contacts and no seismic de- vices are considered in this chapter, providing that the materials remain in the 6.3. Elastic analysis 131 elastic range, the response is completely linear1. Di(cid:27)erent strategies to address the solution of this system of equations are presented below. 6.3.1 Direct Response History Analysis: DRHA DirectResponseHistoryAnalysis(DRHA)isequivalenttononlineardynamics(NL- RHA) (section 2.2.3) in elastic analysis; it is based on the direct integration of the coupled system of dynamics with linearized sti(cid:27)ness matrix (expression (6.2)) em- ploying step-by-step numerical integration schemes. The most important advantage of NL-RHA, the accurate consideration of the full nonlinear problem, is largely lost in this case since no material nonlinearities are present. Nonetheless, many other sources of nonlinear behaviour may be considered even in the elastic range and are captured in a precise way with this procedure, i.e. non-classical damping (for exam- ple in structures with general viscous dampers), geometric nonlinearity, contacts, radiation of energy at the foundations, among others. There are several implicit step-by-step time-integration algorithms available to resolve the system of dynamics, but all of them introduce phase errors which are increased if the vibration period considered is reduced, or if the step-time is in- creased2 (see (cid:28)gure 2.4 in section 2.2.3). The practical highest frequency of interest in the proposed cable-stayed bridges is 20 Hz (section 5.4.4) and the step-time has been considered ∆t = 0.01 s after an optimization study, which implies that the the phase error is approximately 100 % in these higher modes after 10 s of earthquake excitation (precisely when the peak response is typically recorded3), which is no- ticed in the errors of direct dynamic analysis (section 6.3.5). The phase error has been proved to be more important in terms of accelerations and velocities than con- sidering displacements, the bene(cid:28)cial e(cid:27)ect of the damping has been also observed elsewhere [Abaqus 2010]. The broadly accepted integration algorithm proposed by Hilber, Hughes and Taylor (HHT) [Hilber 1977] is employed in this thesis because of its valuable ability to dampen the higher order vibrations, which may cause numerical instabilities and have small contributions in the global response. The HHT operator replaces the actual equilibrium equation (expression (6.2)) with a balance of inertial forces at the end of the time step and a weighted average of the static forces at the beginning and at the end of the time step: 1Strictlyspeaking,inthissectiontheresponseiselastic,butthelinearityislostduetogeometric secondordere(cid:27)ects,thatarealwaysincludedintheanalysis. However,thissourceofnonlinearity wasobservedtobesmallforthestudiedbridges,andtheproceduresbasedonmodaldecomposition linearize the sti(cid:27)ness matrix prior to the extraction of modal properties and the seismic analysis. Consequently,theresponsecouldbeconsideredlinearandelasticifmaterialnonlinearitiesarenot recorded. 2AccordingtoAbaqususer’smanual[Abaqus 2010],astep-timelowerthan10−5sisrequiredto achievelessthan5%phaseerrorinasti(cid:27)responseof62Hzafter10secondsofdynamicresponse. 3Usually,thepeakseismicresponseisrecordedaroundt=10s,duetothemodulatingfunction employed in the accelerograms ((cid:28)gure 5.3). 132 Chapter 6. Seismic analysis (cid:0) (cid:1) mu¨ +(1+α ) cu˙ + ku α (cu˙ + ku ) = F (6.3) |t+∆t a |t+∆t |t+∆t − a |t |t |t˜ F is the force vector and α the parameter controlling the arti(cid:28)cial numerical a damping introduced by HHT. Symbols and mark the beginning and the ·|t ·|t+∆t endofthestep, whereas isevaluatedsomewherebetweenthosepoints, depending ·|t˜ on the numerical damping; t˜= t+(1+α )∆t, where ∆t is the time-step. a The operator de(cid:28)nition is completed by the Newmark’s formulae for displace- ment and velocity integrations: ∆t2 (cid:2) (cid:3) u = u +∆t u˙ + (1 2β ) u¨ +2β u¨ (6.4a) |t+∆t |t |t 2 − a |t a |t+∆t (cid:2) (cid:3) u˙ = u˙ +∆t (1 γ ) u¨ +γ u¨ (6.4b) |t+∆t |t − a |t a |t+∆t HHT is an implicit single-step, second-order accurate algorithm which is uncon- ditionally stable in linear analysis if the parameters α , β and γ satisfy: a a a 1 1 1 β = (1 α )2; γ = α ; α 0 (6.5) a a a a a 4 − 2 − − 3 ≤ ≤ Considering α = 0, the procedure reduces to the trapezoidal rule and the a maximum accuracy is theoretically achieved, since excessive damping of important higher modes is avoided and the phase error is minimized (see (cid:28)gure 2.4). Speci(cid:28)c sensitivity analysis conducted in one of the proposed cable-stayed bridges4 have concluded that the numerical damping factor has a small in(cid:29)uence on the extreme seismic response, and it was observed that eliminating this arti(cid:28)cial dissipation the convergenceisachievedwiththesameagilityasincludingnon-zerovalues. However, Abaqus [Abaqus 2010] suggests α = 0.05 in linear and nonlinear analysis, and a − thisvaluewasselectedforalltheapplicableanalyseshereinafter(indirectdynamics, DRHA, and its extension in the nonlinear range, NL-RHA). The main drawback of this procedure is undoubtedly the large computational cost. Takingintoaccountthatnocontactissimulatedwiththesimpli(cid:28)edconnection between the deck and the towers (section 3.2.3) and that seismic devices are not considered in this chapter, direct dynamic analysis seems hardly justi(cid:28)able if the material response is maintained completely elastic. However, this procedure has been considered in the present section in order to verify its accuracy and to validate its use in subsequent inelastic analyses. As it was stated in chapter 5, constant damping is not possible in direct dynam- ics. Rayleigh damping with the parameters a = 7.48 10−2 and a = 6.32 10−4 0 1 · · (section 5.4.4) has been considered in this methodology, whereas the synthetic ac- celerograms obtained with the same damping dependence on the frequency (section 5.4.4) are employed in direct analysis (DRHA and NL-RHA). 4‘H’-shaped towers with LP =200 m on rocky soil. 6.3. Elastic analysis 133 6.3.2 Modal Response History Analysis: MRHA An e(cid:30)cient way to solve the system of dynamics consists in decoupling expression (6.2) by means of the modal transformation, expanding the displacement vector in terms of modal contributions: 6N (cid:88) u(t) = φ q (t) (6.6) i i i=1 Where φ [6N 1] is a vector with the normalized displacements of each degree i × of freedom associated with the i-mode, q (t) being its generalized coordinate. φ is i i obtained from the real eigenvalue problem5: kφ = ω2mφ (6.7) i i i Substitutingexpression(6.6)in(6.2),andpre-multiplyingeachtermbyφT gives: n 6N 6N 6N (cid:88)(cid:0)φTmφ q¨ (t)(cid:1)+(cid:88)(cid:0)φTcφ q˙ (t)(cid:1)+(cid:88)(cid:0)φTkφ q (t)(cid:1) = φTmιu¨ (t) (6.8) n i i n i i n i i − n g i=1 i=1 i=1 It could be demonstrated that both the mass matrix (m) and the linear sti(cid:27)ness matrix (k) are orthogonal with respect to the eigenvectors φ (i = 1, ,6N) of the i ··· structure. If,inaddition,thedampingisconstantforallthevibrationmodes,orifit is represented by Rayleigh or Caughey distributions, a classical6 damping matrix is obtained, which is also orthogonal to the eigenvectors [Chopra 2007] [Clough 1993]. Because of these orthogonality properties, all the terms in each summation vanish except the i = n components. Therefore, if linear elastic behaviour and classical damping is assumed7, the system of dynamics may be simpli(cid:28)ed up to its transfor- mation in 6N independent SDOF di(cid:27)erential equations, each one representing the dynamic response of a mass-spring-damper system with frequency ω and damping n factor ξ , coincident with the values of the corresponding nth vibration mode. n M q¨ (t)+C q˙ (t)+K q (t) = φTmιu¨ (t) (6.9) n n n n n n − n g Where M , C and K are respectively scalars de(cid:28)ning the generalized modal n n n mass, damping and sti(cid:27)ness of the SDOF associated with the n-mode. Therightpartofexpression(6.9)representsthee(cid:27)ectiveseismicforceassociated with the nth mode, which is applied to the corresponding SDOF. Such force may be expressed in terms of its spatial distribution along the structure, de(cid:28)ned by s = mι 5Damping is ignored in the extraction of vibration modes. 6If dampers are incorporated (chapter 9), the damping matrix is no longer classical, leading to a non-diagonal matrix which is not generally orthogonal [Chopra 2007]. Nonetheless, modal analysismaystillbeemployediftheseismicdevicesarelinearviscoelasticorviscous(cid:29)uiddampers, employing in this case complex mode shapes [Villaverde 2009]. 7Geometricnonlinearitiesareincludedinthestaticsteppriortotheextractionofmodalprop- erties and the modal dynamic analysis, where the sti(cid:27)ness matrix is linearized. 134 Chapter 6. Seismic analysis (independentoftime), andintermsofitsdependencewiththetimebymeansofthe accelerogram. The seismic excitation is generally a vector with three components, however, since the behaviour is completely linear, the response may be studied separatelyineachcomponent,andcombinedafterwardstoobtainthetotalresponse. Therefore, hereinafter in this paragraph, only one component of the accelerogram is considered u¨ (t), and consequently u¨ (t) = u¨ and ι = 1 in expression (6.1) (the g g g generalization of the mathematical background for three-directional excitations is included in section 6.4). Expanding the spatial distribution of the seismic excitation in the base of eigen- vectors which represent the vibration modes: 6N 6N (cid:88) (cid:88) s = mι = s = Γ mφ (6.10) i i i i=1 i=1 Pre-multiplying both sides by φT and considering again the orthogonality of the n mass matrix, the modal participation factors Γ are obtained: n φTmι φTmι Γ = n = n (6.11) n φTmφ M n n n Substituting expression (6.10) in equation (6.9) and dividing by M : n q¨ (t)+2ξ ω q˙ (t)+ω2q (t) = Γ u¨ (t) (6.12) n n n n n n − n g If we consider the change of variable: q (t) = Γ D (t) (6.13) n n n Introducing this new variable in equation (6.12) yields the classical expression of a SDOF with an accelerogram imposed at the base: D¨ (t)+2ξ ω D˙ (t)+ω2D (t) = u¨ (t) (6.14) n n n n n n − g The response of each one of the 6N di(cid:27)erential equations is integrated by means of numerical step-by-step procedures8 to obtain q (t). Once these values are known, n the contribution of the nth mode to the relative displacements of the structure is obtained by reversing the change of variable: u (t) = φ q (t) = φ Γ D (t) (6.15) n n n n n n The nth mode contribution to the response of the structure in time-domain may be obtained by assuming that an equivalent force vector (f (t)) is applied, and n varied, in order to maintain the deformation at every instant. f (t) = ku (t) = kφ Γ D (t) (6.16) n n n n n 8For example the interpolation of the excitation or Newmark’s methods [Chopra 2007]. 6.3. Elastic analysis 135 Now, taking into account the eigenvalue problem (6.7), the expression (6.16) is simpli(cid:28)ed: f (t) = ω2Γ mφ D (t) = s (cid:2)ω2D (t)(cid:3)(t) = s a (t) (6.17) n n n n n n n n n n (cid:124) (cid:123)(cid:122) (cid:125) (cid:124) (cid:123)(cid:122) (cid:125) sn an Generalizing the expression above, any element of the modal response (forces, displacements, etc.), r (t), may be obtained as the product of two factors; (i) the n structuralresponsewhenthespatialdistributionoftheseismicforce(s )isstatically n applied (modal static response, rst); and (ii) the pseudo-acceleration recorded in n time domain for the n-mode SDOF system subjected to the earthquake, a (t): n r (t) = rsta (t) (6.18) n n n Finally, the total response of the structure for the studied accelerogram compo- j nent (u¨ , where j = X,Y or Z) is obtained by aggregating directly the individual g modal contributions during the earthquake: 6N 6N (cid:88) (cid:88) r (t) = r (t) = rsta (t) (6.19) n n n n n=1 n=1 Thegenericstructurewith6N degreesoffreedomhasconsequently6N vibration modes. Nevertheless, only the (cid:28)rst modes have a signi(cid:28)cant contribution in the global response and, therefore, the analysis may be conducted considering just the (cid:28)rst J vibration modes (J < 6N), truncating the summation and obtaining: J J (cid:88) (cid:88) r (t) r (t) = rsta (t) (6.20) n ≈ n n n n=1 n=1 It is neither reasonable nor interesting to include all the vibration modes; the modal contribution of higher-order modes was studied in chapter 4, whereas the higher interesting frequency was discussed in some detail along chapter 5. From the analysis point of view, unacceptably small time steps should be imposed to avoid the (cid:28)ltration of signals with very high frequencies, furthermore, very small elements need to be employed in the FEM discretization to accurately represent the mode shape of these complex modes (they have many points with zero displacements). Without losing sight of the negligible contribution for extremely high-order modes (Γ is very reduced), removing them from the calculation is advisable. However, n a more restrictive condition than the one imposed by Eurocode 8 [EC8 2005a] is considered in modal dynamic analysis conducted here; the frequency of the last mode included is f = 35 Hz. lim If the structure is three-axially excited, the responses obtained by considering each accelerogram component separately are directly aggregated in time domain to obtain the total response. Modal dynamic analysis (MRHA) is far less computationally expensive than direct integration (DRHA) introduced above since J independent static analyses 136 Chapter 6. Seismic analysis are performed (to obtain rst), besides J dynamic analyses of SDOF equations (to n obtain a (t)), being J < 6N. Unfortunately, due to simpli(cid:28)cations assumed to n decompose the system of dynamics in expression (6.9), modal dynamic analysis is strictly valid if the materials remain completely elastic during the earthquake and if there are not any other source of nonlinearity included in the model. The errors introduced with modal dynamics in elastic analysis are basically; (i) the error causedbythenumericalintegrationprocedureoftheSDOFmotion; (ii)thepossible (cid:28)ltering of higher modes by the step-time adopted; and (iii) the vibration modes neglected (between J and 6N). Compared with other procedures, the sources of errors are less detrimental here, since the employed integration scheme is exact if the step-time of the analysis and the accelerogram are coincident9 (∆t = 0.01 s). Furthermore, the results remain almost unchanged if an increased number of modes is included in the analysis, as it has been observed in this section. Therefore, modal dynamic analysis (MRHA) is considered as the reference methodology for comparison purposes in the elastic analysis performed in this chapter. In modal dynamics, a constant modal damping ratio ξ = 4 % (n = 1, ,J) is n ··· considered, whereas the synthetic accelerograms obtained with this constant damp- ing level (section 5.4.3) are employed in the analysis. 6.3.3 Modal Response Spectrum Analysis: MRSA In structural design, the most interesting result is the peak value of the response r , o notitstime-domainhistoryr(t). ConsideringaSDOFsystemwithfrequencyω and n damping ratio ξ , its maximum response under a speci(cid:28)c earthquake may be known n exactly by means of the acceleration response spectrum Sa = Sa(ω ,ξ ), being n n n Sa the pseudo-acceleration spectral ordinate10 or, in other words, the extreme n relative displacement of the SDOF during the earthquake (Sd ) times the square of n the corresponding frequency. Returning to the study of the system of dynamics, and assuming the response perfectly linear and elastic (without any seismic device), we arrived previously to expression (6.18), repeated here for convenience: r (t) = rsta (t) = rst(cid:2)ω2D (t)(cid:3) (6.21) n n n n n n The maximum response of nth mode is: (cid:104) (cid:105) r = rst ω2max D (t) = rst(cid:2)ω2Sd (cid:3) = rstSa (6.22) n,o n n t | n | n n n n n 9Abaqus [Abaqus 2010] states that the integration scheme in MRHA is exact if the seismic ex- citation(whichisproportionaltothegroundacceleration,expression(6.1))varieslinearlybetween twoconsecutivestepsand,therefore,theprocedureisexactifthestep-timeintheanalysisisequal to the step of the accelerogram. 10Strictly speaking, it is employed the pseudo-acceleration spectrum San = ω2Sdn, not the acceleration spectrum max u¨T , however it is referred simply as the acceleration spectrum since | | max u¨T,n Sanconsideringthetypicaldampingfactorsinearthquakeengineering[Chopra 2007]. | |≈