UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
NONLINEAR DYNAMIC ANALYSIS OF BRIDGE PIERS By CESAR FERNANDES, JR. A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA Copyright 1999 by Cesar Femandes, Jr. I would like to dedicate this dissertation to my parents, Cesar and Dalva, my brothers, Magno and Marcus, and to my Fiancee, Leandra. I could not have reached such an accomplishment without their help. ACKNOWLEDGMENTS I would like to thank my parents, Cesar and Dalva, my brothers, Magno and Marcus, and all my family, for their unconditional support in all phases of my graduate studies. I also would like to thank my fiancee, Leandra, for being so patient and supportive during all my work. I also would like to thank all the faculty from the Civil Engineering Department at the University of Florida, and the members of my supervisory committee Dr. Hoit, Dr. McVay, Dr. Fagundo, Dr. Hays, and Dr. Wilson for always having their doors open to answer my questions. Finally I would like to thank all the graduate students from the workstation lab, specially, Mark and Wirat, for their support and friendship. TABLE OF CONTENTS page ACKNOWLEDGMENTS ......................................... ....................... iv LIST O F TA B L ES ......................................................................... ............................ viii L IST O F FIG U R E S ............................................................................ ...........................x A B STRA CT ............................... ............................................... xvi CHAPTERS 1 IN TR O D U C TIO N ....................................................................... ............................1 B background ................................................................... ........................................... Literature R review ........................... .... ....................................................................3 Lim stations ................................................................................ ......................... 7 O rganization................................... .................................................................... 7 2 NONLINEAR DYNAMIC ANALYSIS ........................ .......................9 Theory ................................................ ....................................... ......................... 9 Equations of Motion, Mass, and Damping Matrices ................................................... 10 Equations of Motion for Ground Motion............................. .....................14 Mass Matrices, Consistent and Lumped ................................. ......................18 Mass Matrix for the Uniform 3DBeam Element..........................................18 C onsistent............................................ .................................................18 Lum ped ................................................... ............................................ 20 Mass Matrix for the Shell Element.................................. .........................22 C onsistent............................................ ................................................. 22 Lum ped .............................................. .................................................. 25 Remarks about the mass matrix ................................. ....................26 D am ping............................................................................... ..................................27 Estimating Modal Damping Ratios............................. ... ........................29 M ass Condensation ........................................ ... ............................................. 30 TimeHistory Analysis. Direct Integration Methods................................................34 Numerical Evaluaton of Dynamic Response. Newmark's Method...........................36 Choice of tim e step At.......................... ..... ................................................ 40 N onlinear Problem s ....................................................................... ...................... 41 Analysis of the Nonlinear Response using Newmark's Method ........................41 Nonlinear Dynamic Analysis Algorithm ..............................................................47 3 DISCRETE ELEMENT MODEL AND MATERIAL HYSTERESIS .....................51 Discrete Element Derivation....................... ........... .........................51 Element Deformation Relations.......................... ....... .......................52 Integration of Stresses for Nonlinear Materials......................... ......................54 Elem ent End Forces ............................................ ............................................. 58 Elem ent Stiffness .................................................................... .........................59 Secant and Tangent Stiffness of the Discrete Element..............................................59 H ysteresis M odels................................................ ............................................ 60 M material M odels.............................................. .................................................62 Uniaxial Mild Steel Model.................................................................63 Uniaxial Monotonic Concrete Model Used in FLPIER.......................................64 Proposed Models for the Uniaxial Inelastic Cyclic Behavior of Concrete..................67 Rational M odel................................................ ...........................................67 L loading ....................................................................... ...............................68 U nloading........................................... .................................................. 72 Reloading ............................ ..... .................. ........................... 75 B ilinear M odel ............................................ ...............................................77 Strain Rate Effect.............................................. ...............................................78 Confinem ent Effect.......................... ... ..................... ........................... 79 4 M ODAL ANALYSIS .............................................................. ........................83 Natural Vibration Frequencies and Modes ................................. ...........................83 Modal and Spectral Matrices ...................... ....................................86 Normalization of Modes ....... ....................................................... .......................87 M odal Equations ..................................................................... ......................... 88 Elem ent Forces.................................................. ............................................... 92 Modal Equations for Ground Motion .................................. ........................92 Response Spectrum Analysis................................ ........................................... 95 Modal Combination Rules.............................................................96 How FLPIER handles Modal Analysis...................... ... ...............................99 5 MULTIPLE SUPPORT EXCITATION......................... .................................101 6 SOIL STRUCTURE INTERACTION ............................... .........................108 U coupled M ethod.............................................. ...........................................109 Coupled Method...................................... .................. ................ 10 Cyclic Behavior of Soil........................................................ .......................... 12 Cyclic D degradation ............................................................... .........................113 Strain R ate E ffect.................................................................. ......................... 14 Radiation D am ping............................................................... ......................... 15 vi 7 PREDICTIONS OF RESPONSE ..................................... ....................... 116 Example 1 Steel Section 1 ........................................................116 Example 2 Steel Section 2...................................................... ......... 118 Example 3 Circular Reinforced Concrete Column 1 .............................................. 121 Example 4 Circular Reinforced Concrete Column 2........................................... 125 Example 5 Rectangular Reinforced Concrete Column........................................ 128 M onotonic Tests........................ ..................... .............. ............129 C yclic Test SO ............................................... ........................................... 136 C yclic Test S ............................................. .............................................142 Cyclic Test S2.......................... .... ................... ............................ 146 Cyclic Tests S3 and S4 ..................... ....................................... 148 Cyclic Test S10.............................................................. .........................153 Example 6 Piles in Sand.................................................. ............... 159 T est SP ........................................................................... ................................159 Tests PG2 and PG3 .................. ....................................................... 162 Example 7 Mississippi Dynamic Test...................................................................166 8 CON CLU SION S .................................................................. ......................... 172 APPENDICES A M A SS UNITS............................................. ..................................................175 B GAUSS QUADRATURE ..........................................................178 C FLPIER Manual for dynamic analsyis........................... .. ....................184 REFERENCES ....................................................................................................184 BIOGRAPHICAL SKETCH ...........................................................191 vii LIST OF TABLES Table page 21. Natural frequencies of a uniform cantilever beam: ConsistentMass Finite Element and exact solution ................................................................. ......................... 21 22. Natural frequencies of a uniform cantilever beam: LumpedMass Finite Element and Exact Solution.................................... .....................................................21 23. Recommended damping rations for structures ................... ..........................30 31. Curves coefficients ............................................................... .......................... 73 71. Results for Exam ple 2............................................................ ....................... 119 72. Design details for Example 3..........................................................121 73. Model parameters for Example 3 ................................ .........................121 74. Design details for Example 4....................... .... ................................125 75. Loading for each test .......................................................................128 76. Parametric tests, units are KN/mm2....................................................... .............131 77. Parametric tests for confinement under different strain rates.................................135 78. Parameters for test SO units are KN/mm2.......................... ........................ 139 79. Parameters for test S1 units are KN/mm2....................... .............................143 710. Parameters used in test S3 and S4...................... ....................149 712. Some soil properties for CSP1.........................................................160 713. Masses (tons) for pile group tests...................... ... ...............................162 A 1. M ass density units.............................................................................. .............. 177 B1. GaussLegendre abscissas and weights ................................. ............................182 LIST OF FIGURES Figure page 11. Bridge pier com ponents........................................................ ...........................2 21. Tower subjected to ground motion after Chopra (1995).........................................15 22. Support motion of an Lshaped frame. a) Lshaped frame; b) influence vector f: static displacements due to Dg=1; c) effective load vector after Chopra (1995).................................................. 17 23. 3D B eam elem ent ............................................................... .................................19 24. True 9node rectangular element.................................................22 25. Mapping for a true rectangular 9node shell element...............................................22 26. Shell element of uniform thickness .................................................. ...................23 27. Lumped mass matrix at the nodes of true rectangular 9node shell element. Numbers shown are fractions of the total element mass at each node.....................26 2.8. Full and condensed versions of the structure............................ ........................34 29. Average acceleration ............................................................ ......................... 37 210. Linear acceleration ............................................................. .........................37 211. Secant and Tangent approaches. After Chopra(1995)...................................... 43 212. NewtonRaphson Method............................. ..... ........................45 31. Representation of discrete element. After Hoit et al., 1996 ....................................51 32. Discrete element displacements. After Hoit et al. 1996 ..........................................54 33. Various components of total strain in the section. After Hoit et al., 1996 ................55 34. Rectangular section with integration points. After Hoit et al., 1996.......................56 35. Circular section with integration points. After Hoit et al, 1996 ..............................57 36. Secant and tangent material stiffness............................ ........................ 60 38. Elasticperfectly plastic model for mild steel......................................................64 39. FLPIER concrete points......................... ............................ ...........................65 310. Concrete strains ................................................................... ........................66 311. Concrete stresses.................................... ...................................................... 66 312. Envelop curve for concrete.......................... ...... ........................68 313. Typical compression loading........................ ...... ..........................70 314. Typical loading in tension ...................... ......................................71 315. Typical unloading in tension............................................................73 316. Typical unloading in compression............................ .........................74 317. Compression unloading with gap ....................... .. ................................74 318. Compression unloading with no gap................................. ......................75 319. Typical loading, unloading and reloading in compression........................... ...76 320. Concrete behavior with gap............................... .... ..........................76 321. Bilinear model for concrete..........................................................77 322. Stressstrain curve for concrete in FLPIER........................................................78 323. Confined and unconfined concrete models response ............................................. 79 324. Core width for different cross sections........................... ....................... 81 325. Confined concrete model...................... ....... ..........................82 41. Generalized SDF system for the nth natural mode............................................90 42. Typical response spectrum ...................... .....................................96 43. M odal analysis of pier ........................................................... ........................ 99 5.1. M multiple support m option ................................................... ........................... 102 52. 2D frame submitted to multiple support motion ................................................. 105 53. Pile subjected to multiple support excitation.....................................................106 61. a)Coupled model; b)Uncoupled model.............................. ........................11 62. Typical py curve................................................................. ...... ........................ 12 63. Cyclic soil m odel...................... ......... ..............................................................113 71. Exam ple 1 com puter m odel.....................................................................................117 72. Comparison FLPIER x reference for cyclic loading ............................................. 117 73. Exam ple 2 com puter m odel.................................................................................... 118 74. FLPIER x Hays, shear force comparison .................................. ........................119 75. FLPIER x Hays, bending moment comparison......................................... ..........120 76. Shear comparison Example 3 x FLPIER........................ .................................123 77. Moment comparison Example 3 and FLPIER.....................................................124 78. Dynamic model comparison................................ .......................124 79. Static shear comparison Example 4....................................................... 126 710. Dynamic moment comparison Example 4. Hoop spacing = 5 in........................127 711. Dynamic moment comparison Example 4. Hoop spacing = 2 in........................127 712. Computer model for Test 3..........................................................130 713. Column capacity using FLPIER.......................... ...... ....................130 714. E and E, changed........................................ ..............................................132 715.fc and f changed ............................................................... ....................... ...... 132 716. E andf changed.................................................................... .......................133 xii 717. All parameters changed .................... .........................133 718. Confinement and Strain rate effect............................. .......................135 719. Imposed displacement history for the first 90 seconds........................................137 720. Comparison FLPIER x Test, original properties..............................................138 721. Comparison FLPIER x Test, modified properties................................................138 722. Imposed tip displacement history for test.......................... ............ ..................140 723. Com prison tsOl ........................................ ..............................................140 724. Com prison ts02............................................................................................141 725. Com prison ts03 ............................................................. ................................... 141 726. Comparison ts04........................ ................. ............................................142 727. Imposed displacement in X direction..................................................................142 728. Imposed displacement in Y direction......................... .............................143 729. Com prison tsl 1........................................ ..............................................144 730. Comparison tsl c ....................................................................................144 731.Com prison ts12 ........................................... .............................................. 145 732. Com prison tsl2c ............................................................ ..........................145 733. Imposed displacements in X direction....................... ................................ 146 734. Imposed displacements in Y direction......................... ...............................147 735. Comparison test S2 x FLPIER ................................... ........................147 736. Imposed X displacement for test S3.................................................................... 149 737. Imposed Y forces for test S3 ...................................................149 738. Comparison Test S3 x FLPIER ts31 ............................. .........................150 739. Comparison Test S3 x FLPIER ts32 .............................. .........................150 xiii 740. Imposed displacement in X direction for test S4..................................................151 741. Imposed load in Y direction for test S4.......................... ............................151 742. C om prison ts4 ............................................................... ................................. 152 Fig 743. Com prison ts42....................................... ............................................ 152 744. Com prison ts43 ............................................ ............................................153 745. Imposed force in Xdirection for S10...................................................................154 746. Imposed force in the Zdirection for S10 ........................ ............ ............ 155 747. Com prison tslO0 x test ................ .................................. ......................... 155 748. Com prison ts 02 x test .................................................. ... ....................... 156 749. Com prison tsl03 x test ..................................................... ....................... 156 750. Comparison ts104 x test .........................................................157 751. Single pile in sand............................................................ .......................... 160 752. Top ring acceleration........................................................ ......................... 161 753. Top displacement comparison........................ ...... ......................161 754. Moments comparison ................. ............... ........................................162 755. 2 x 2 pile group in sand ..........................................................163 756. 3 x 3 pile group in sand ............................................................................163 757. Top of pile lateral displacement comparison, 2 x 2 group ..................................64 758. Top of pile bending moment comparison, 2 x 2 group .........................................164 759. Top of pile bending moment comparison, 3 x 3 group .......................................165 760. M ississippi test structure .................................................... .......................167 761. Load history for Mississippi test ................................ .........................168 762. Displacement history for Mississippi test.......................... ......................... 168 xiv 763. Comparison test x FLPIER with CPT soil properties..........................................169 764. Comparison test and FLPIER with SPT soil properties ......................................169 765. Comparison test and FLPIER with original soil properties.................................170 B1. (a) Trapezoidal approximation using the abscissas 1 and 1. (b) Trapezoidal approximation using abscissas x, and x2 after Mathews, 1987............................179 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy NONLINEAR DYNAMIC ANALYSIS OF BRIDGE PIERS By Cesar Femandes, Jr. August 1999 Chairman: Marc Hoit Major Department: Civil Engineering Bridge piers are often subjected to lateral loading that is not neglectible when compared to vertical loads. Such loading conditions may include wind, water and earthquakes. In order to develop a time domain analysis for the nonlinear dynamic response of piers and their foundations the computer program FloridaPier was modified. FloridaPier is a nonlinear finiteelement program, developed at the University of Florida, designed for analyzing bridge pier structures composed of pier columns and cap supported on a pile cap and piles with nonlinear soil. The program was developed in conjunction with the Florida Department of Transportation (FDOT) structures division. The piers and the piles are modeled as nonlinear 3D beam discrete elements. These elements use the true material stressstrain curves (steel and concrete) to develop its behavior and stiffness modeling. Nonlinear dynamic capability was added to these elements by adding to the stressstrain curves the ability to represent loading, unloading and reloading behavior, typical of dynamic loading. In addition, a mass matrix for this element was also implemented. The mass matrix for the piles cap was also developed. The cap is modeled as linear shell elements. The stiffness matrix for this type of element is considered to remain unchanged during the dynamic analysis. A nonlinear model was also added for the soil. The py springs generated by the program now have the capability of loading, unloading and reloading, just like the steel and concrete. Such models can then be applied to dynamic analysis of reinforced concrete structures subjected to seismic, impulsive, or wind loads. In addition, being the current stateofpractice, modal analysis was also implemented. Seismic load can now be applied to the linear structure considering the nonlinear behavior of the foundation. A description of the analysis used to model the nonlinear dynamic behavior of bridge piers, as well as the implemented nonlinear material behavior, is presented in this dissertation. CHAPTER 1 INTRODUCTION Background Extensive research efforts have been directed to the nonlinear response of structures subjected to extreme load events. These extreme load events could be an earthquake or hurricane for a building, a ship impact for a bridge, or the effects of waves and wind action for offshore oil platforms. Traditionally, large factors of safety have been used in such cases, resulting in overconservative design and cost ineffectiveness. On the other hand, an unsafe design could result in catastrophic human and economic losses. Because a more sophisticated nonlinear dynamic analysis is computationally expensive, these structures are designed using factored static loads to account for the dynamic effects. This procedure is acceptable for very low frequency vibrations, however the introduction of nonlinearity, damping, and pilesoil interaction during transient loading may significantly alter the response. Because in recent years computers have become much faster and cheaper, it has become possible to consider, and consequently to study, the dynamic nonlinear behavior of structures considering many factors neglected in the past. In this dissertation, the computer program Florida Pier (Hoit et al., 1996), which will be referred simply as FLPIER from now on, has been modified to allow the nonlinear dynamic analysis of bridge piers. FLPIER is a computer program based on the Finite Element Method developed by Drs. Hoit, Mcvay, and Hays at the University of Florida for the nonlinear static analysis of bridge piers. Nonlinear aspects of structural analysis, such as material and geometric nonlinearity, as well as structuresoil interaction can be incorporated into the analysis leading to more accurate results. It can model all the components of a bridge pier and its foundation, such as pier, pier cap, piles cap, piles, and soil, as shown in Fig.11. The pier, pier cap, and piles can be represented using discrete elements that can incorporate the effects of material and geometric nonlinear behavior. More details about the discrete element are found in Chapter 3. The piles cap is modeled as linear 9 node shell elements. The lateral soil resistance is modeled as nonlinear py springs, while the axial resistance is modeled as nonlinear tz springs. Pier Cap Bridge Springs e Piers Substructure Pile Cap Piles Soil layers Fig. 11. Bridge pier components FLPIER is now used by many DOTs throughout the United States because of its reliability and ease of use. Unlike other general Finite Element programs, like ADINA and SAP, where modeling and analyzing can be time consuming, in FLPIER it is easy and fast for the user to perform these tasks thanks to a user friendly interface for model generation. The modification of soil or structure parameters in the model is not difficult either. The results can also be seen through a graphic interface that is currently being updated. In the modified dynamics version, resulting from this research, speed and ease were maintained, allowing the user to easily perform the nonlinear dynamic analysis and modify parameters in the soil or structure if necessary. Although the program is more suitable for the analysis of bridge piers, other types of structures can also be modeled. The new contribution for the field is the proposed concrete model. This model was implemented in the FLPIER code to allow the nonlinear dynamic analysis of reinforced concrete sections. Literature Review Over the last years different analytical models have been proposed for the analysis of reinforced concrete structures. Models for these types of structures, which are under primarily flexural and axial loads, can be classified as: (i) Simple or lumped models. (ii) Discrete models. (iii) Fiber models. (iv) Finite element models. Singledegreeoffreedom (SDOF) models belong to the first class of analytical models. In this class of models it is assumed that the structure's response to an earthquake is dominated by its first natural frequency, allowing the system to be represented as a SDOF system with lumped mass and stiffness properties (Crandall (1956), Craig (1981), Paz (1985), and Chopra (1995)). A more general representation for multidegreeof freedom systems (MDOF) is derived using the concept of shear building. In this model the stiffness of each story is represented by nonlinear springs, and the beams are considered to be infinitely rigid. Despite its simplicity and satisfactory performance in predicting the maximum response, this class of models does not provide enough data for more detailed seismic analysis. Furthermore for more complicated frames the assumption that the beams are infinitely rigid may not be correct. In the second class of analytical models, the discrete models, there is a correspondence between the analytical model and the actual structure. In such models, a linear elastic element and a nonlinear spring represent the structural elements. The most common case is that of a nonlinear spring attached to both ends of a linear beam element. Atalay (1975), Clough (1966), Nakata et al. (1978), Park (1984), and Takeda (1970), among others, have extensively used this class of models to analyze the behavior of reinforced concrete structures. In these models a set of predefined rules defines the hysteretic behavior of the nonlinear springs. These rules are usually obtained from laboratory experiments with real scale specimens. It is mainly the difference between these rules that distinguishes the models. Although these models give satisfactory results, its main disadvantage is the fact that the nonlinear spring's rules are based on experiments that may not correspond to the actual structural member, or type of loading, that they are representing. The fiber models have been used in the study of reinforced concrete (Ala Saadeghvaziri (1997), Hajjar et al. (1998), Park et al. (1972), and Zeris and Mahin (1991a), Zeris and Mahin (1991b)) and steel members (Baron and Venkatesan (1969), Chen and Atsuta (1973)). These models are based on the finite element approach, and are better suited for members and structures under complex loading histories. In these models the crosssection is divided into segments. Each segment can then be divided in one or more fibers. Each fiber is assumed to obey a uniaxial stressstrain relationship. From the integration of the stresses of each fiber over the cross section, the element forces can be calculated, and from the evaluation of the stiffness of each fiber the overall element stiffness can be obtained. Once the element forces and stiffness are obtained the analysis is carried out using standard Finite Element Method procedures. Therefore only the stressstrain relationships for concrete and reinforcing steel in the case of reinforced concrete sections, or steel, in the case of steel sections, are necessary to describe the properties of each section of the element. This makes these models very effective under complex loads. The main difference among all the fiber models are the rules adopted for the uniaxial behavior of the different materials that make the cross section. In the case of most civil engineering structures, steel and concrete, but other materials can also be used if the stressstrain relationships are known. In the specific case of concrete the model's backbone is the envelope curve obtained from a monotonic test. This curve limits the concrete stresses in any loading phase. In some models the compression envelope curve for concrete is represented by the wellknown Hognestad parabola (Ala Saadeghvaziri (1997), Park et al. (1972)). Another approach is to use multilinear curves to define concrete behavior in compression (Zeris and Mahin (1991a), Zeris and Mahin (1991b)). In the case of dynamic analysis, the unloading and reloading rules are particular for each author. In the case of Ala Saadeghvaziri (1997), and Zeris and Mahin (1991a), Zeris and Mahin (1991b), unloading under compressive stress has a slope equal to the initial Young's modulus of the material. However Hajjar et al. (1998) and Park et al. (1972), have their own more complicated expressions for the unloading curve. The reloading curves are very specific for each model, and the reader is referred to the mentioned references for more information. The tension strength of concrete can be neglected (Zeris and Mahin (1991a), Zeris and Mahin (1991b)), or assumed to be equal to the concrete tensile strength, (Ala Saadeghvaziri (1997), Park et al. (1972)) in which case its slope is assumed to be equal to the initial slope of the compression side. The unloading and reloading criteria are again specific for each model, and the reader is referred to the above references for more information. Confinement, strainrate, and stiffness degradation effects are also particular for each model. In the case of steel, no distinction is made between steel sections and reinforcement steel. The rules are valid for both cases. For the stressstrain curves, the Baushinger's effect can be considered (Park et al. (1972), Baron and Venkatesan (1969)), or ignored, in which case a bilinear or trilinear relationship (elastoplastic with kinematic or isotropic hardening) is used (Ala Saadeghvaziri (1997), Chen and Atsuta (1973), Zeris and Mahin (1991a), Zeris and Mahin (1991b)). The major disadvantage of these models is that they are computationally very expensive, but with the recent advances in computer technology, this class of models has become more popular because of its versatility. This is the analytical model used in this work. FLPIER already incorporates a nonlinear discrete element, which uses fiber modeling at two points along the element's length to characterize its nonlinear behavior. New stressstrain curves were introduced to allow the nonlinear discrete element to perform dynamic analysis. The details can be found in Chapter 3. The last class of analytical methods is the finite element method. In this class of methods different elements are used to represent the structural members, such as truss members to represent the reinforcing steel, and plane stress elements to represent the concrete. The cracking typical of concrete represents a computational difficulty for these models, requiring the development of more sophisticated elements. The development of such elements is a challenge based on complex elasticity and plasticity theories. Like fiber modeling, this class of methods is computationally expensive and time consuming. Limitations It is very difficult for a model to incorporate all the aspects inherent to nonlinear dynamic analysis, and the model presented here is no exception. The first limitation is the fact that all the theory developed in Chapter 2 for nonlinear dynamic analysis is based on small displacement theory. The second limitation comes from the fact that the effect of shear deformation is not included in the constitutive models. It was also not included in the original derivation (Hoit et al., 1996). Problems of local buckling are also outside the scope of this work. Organization In Chapter 2 all the theory necessary for the formulation of the problem is presented. In Chapter 3 the derivation of the discrete element is described and the 8 adopted constitutive models for concrete and steel are presented. Then Chapter 4 discusses the actual stateofdesign procedure for modal analysis. In Chapter 5 the concept of multiple support excitation is introduced. In Chapter 6 the soil structure interaction and the dynamic soil behavior are explained. In Chapter 7 the response predicted by FLPIER is compared to various literature results. Finally in Chapter 8 the conclusions and suggestions for future work are discussed. CHAPTER 2 NONLINEAR DYNAMIC ANALYSIS Theory In a static problem the frequency of the excitation applied to the structure is less than one third of the structure's lowest natural frequency. In this case the effects of inertia can be neglected and the problem is called quasistatic. For such problems the static equations [K]{D} = {R} are sufficiently accurate to model the response, even though the loads R} and displacements {D vary (slowly) with time. The static loads {R) may result from surface loads and /or body forces. On the other hand if the excitation frequencies are higher than noted above or if the structure vibrates freely, the inertia effects must be considered in the analysis. The inertia effects are accounted for by the mass matrix, written as [m] for an element and [M] for a structure, which is a discrete representation of the continuous distribution of mass in a structure. The effects of damping, if important, are accounted for by the damping matrices [c] and [C]. The dynamics problems can be categorized as either wave propagation problems or structural dynamic problems. In wave propagation problems the loading is often an impact or an explosive blast. The excitation and the structural response are rich in high frequencies. In such problems we are usually interested in the effects of stress waves. Thus the time duration of analysis is usually short and is typically of the order of a wave transversal time across a structure. A problem that is not a wave propagation problem, but for which inertia is important, is called a structural dynamics problem. In this category, the frequency of excitation is usually of the same order as the structure's lowest natural frequencies of vibration. A typical example of a wave propagation problem would be that of analyzing the stresses in a pile when it is grounded. It is important not to exceed the allowable stresses in order not to damage the pile. Earthquake analysis of structures is a typical structural dynamics problem, where the inertia forces govern the response of the structure. Problems of structural dynamics can still be subdivided into two broad classifications. In the first one, we are interested in the natural frequencies of vibration and the corresponding mode shapes. Usually, we want to compare natural frequencies of the structure with frequencies of excitation. In design, it is usually desirable to assure that these frequencies are well separated. In the second classification, we want to know how a structure moves with time under prescribed loads, like under impacts, blasts or wind loads, and/or motions of its supports, like in the case of an earthquake. We are interested in the timehistory analysis. The two most popular methods of dynamic analysis are modal methods and direct integration methods. Equations of Motion, Mass, and Damping Matrices The equations that govern the dynamic response of a structure will be derived by requiring the work of external forces to be absorbed by the work of internal, inertial and viscous forces, for any small kinematically admissible motion (i.e., any small motion that satisfies both compatibility and essential boundary conditions). For a single element, this work balance becomes I{.} I{F}dV + {Su} {(Q}dS + (8 {}, (p ', S Eq. 2. 1 = ({8 {}+{8u} p{i}+{6.} c,{u}) V Ve where {8u} and {(e} are respectively small arbitrary displacements and their corresponding strains, {F} are body forces, {(} are prescribed surface tractions (which are typically nonzero over only a portion of surface Se), {P}, are concentrated loads that act at total of n points on the element, {8u}, 'is the displacement of the point at which load {p}i is applied, p is the mass density of the material, Cd is a materialdamping parameter analogous to viscosity, and the volume integration is carried out over the element volume V,. Using usual Finite Element notation, we may write the continuous displacement field {u}, which is a function of both space and time, and its first two time derivatives, as {u} =[N]{d} {u} = [N]{d} {ii} = [N]{d Eq.2.2 In Eqs. 2.2 the so called shape functions [N] are functions of space only, and the nodal DOF {d} are functions of time only. Thus Eqs. 2.2 represent a local separation of variables. Combination of Eqs. 2.1 and 2.2 yields {8d})' f[B {ao}dV+ Jp[N]'[N]dVd}j + kd [N]T[N]dVd}j ve v' V Eq. 2.3 I[N]j{F}dV [Nj]Tf{dS {p}, 0 Ve Se I in which it has been assumed that the concentrated loads {p}j are applied only at the nodal points locations. Since {Sd} is arbitrary, Eq. 2.3 can be written as [m{} + [c]{d} + {r"} = {r"'} Eq. 2.4 where the element mass and damping matrices are defined as [m]= p[N] [N]dV Eq. 2.5 [c] = cd[N]T[N]dV Eq. 2.6 and the element internal forces and external loads vector are defined as {rin" = [B]T{a}dV Eq. 2.7 {r"'} = [N]' {F}dV + f[N]' {D}dS + ( {p} Eq. 2.8 Ve Se =1 Equation 2.4 is a system of coupled, secondorder, ordinary differential equations in time and is called a finite element semidiscretization because although the displacements, {d}, are discrete functions of space, they are still continuous functions of time. Methods of dynamic analysis focus on how to solve this equation. Modal methods focus on how to uncouple the equations, transforming the NDOF coupled system into N uncoupled SDOF systems, each of which can be solved independently of others. More details about this method of analysis can be found in Chapter 4. Direct integration methods discretize Eq. 2.4 in time to obtain a sequence of algebraic equations. Structure matrices [M], [C], and {R"'t} are constructed by standard Finite Element Method procedures, i.e. conceptual expansion of element matrices [m], [c], and {rn"'} to "structure size" followed by addition of overlapping coefficients, in the same way it is done for assembling stiffness matrices in static problems. However, the exact manner in which {R'"} is computed depends on the dynamic analysis procedure adopted. When Eqs. 2.5 and 2.6 are evaluated using the same shape functions [N] as used in the displacement field interpolation (Eqs. 2.2), the results are called consistent mass and consistent damping matrices. These matrices are symmetric. On the element level, they are generally full, but on the structure level, they have the same sparse form as the structure stiffness matrix. When p and Cd are nonzero, consistent matrices [m] and [c] are positive definite. Using the mass matrix for example, the kinetic energy S{rd}[m]{d} is positive definite for any nonzero {d}. In typical structural analysis we are more interested in dry fiction and hysteresis loss, than in viscous damping. It is still not well understood how the damping mechanisms develop in structures, so from a practical standpoint Eq. 2.6 does not correctly represent structural damping. The internal force vector, Eq. 2.7, represents loads at nodes caused by straining of material. Equations 2.4 and 2.7 are valid for both linear and nonlinear material behavior; that is, in Eq. 2.7, {a} could be a nonlinear function of strain or strain rate. For linearly elastic material behavior, {a } = [E][B]{ d} and Eq. 2.7 becomes ri"'} =[k]{d) Eq. 2. 9 where the usual definition of the stiffness matrix holds, that is, [k]= f[B]T[E][B]dV Eq. 2.10 ee e e when Eq. 2.10 is used, Eq. 2.4 becomes [m]{d} +[c]{d} +[k]{d) = {r')} Eq. 2. 11 which can be interpreted as saying that external loads are equilibrated by a combination of inertial, damping, and elastic forces. For the assembled structure, from Eq. 2.11 we get the equation of motion for linear systems, [M]{D} +[C]{D) + [K]{D} = {R'} Eq. 2.12 where {R"et} corresponds to loads {R} of a static problem, but is in general a function of time. Or, returning to Eq. 2.4, equations of the assembled structure can be written in the alternative form [M]{}D +[C]{b} + {Rt} = (R'} Eq. 2.13 which does not require that the material be linearly elastic and represents the equation of motion for nonlinear systems. Equations of Motion for Ground Motion It is now opportune to derive the equations of motion for structural systems subjected to ground motion. Consider the tower shown in Fig. 21, modeled as a cantilever beam with concentrated masses at the nodes The displacement of the ground is denoted by Dg, the total (or absolute) displacement of the mass mj by DLj and the relative displacement between this mass and the ground by Dj. At each instant of time these displacements are related by D (t) = D, (t) + D, (t) Eq. 2.14 Such equations for all the n masses can be combined in vector form: {D}'(t)= {D)(t)+ {D),(t)[1] where [1] is a vector of order n with each element equal to unity. n D Rigidbody Motion D, Fig. 21. Tower subjected to ground motion after Chopra (1995) Only the relative motion [D] between the masses and the base due to structural deformations produce elastic and damping forces (i.e. the rigidy body component of the displacement of the structure produces no internal forces). Thus Eq. 2.13 is still valid, however the inertia forces are related to the total acceleration {f' and from Eq. 2.15 we can write {b}'(t) = {b}(t) + {D (t)[1] Eq. 2.16 and substituting this value back into Eq. 2.13 we obtain [M]{b +[C]{) + R'"') = } Eq. 2.15 Eq. 2.17 The external force vector now becomes the effective earthquake forces vector, and is given by P, = [M]{i,} Eq.2.18 A generalization of the preceding formulation is useful if all the DOF's of the system are not in the direction of the ground motion, or if the earthquake excitation is not identical at all the structural supports (see Chapter 5 for more details). In this general approach the total displacement of each mass is expressed as its displacement DYj due to static application of the ground motion plus the dynamic displacement D, relative to the quasistatic displacement: {D}'(t)= {D) + {D}(t) Eq.2. 19 The quasistatic displacements can be expressed as {D}'(t)= {D}) (t), where the influence vector i represents the displacements of the masses resulting from static application of a unit ground displacement; thus Eq. 2.19 becomes {D}'(t) = {D}(t) + {D}, (t) Eq. 2.20 The equations of motion are obtained as before, except that Eq. 2.20 is used instead of Eq. 2.15, resulting in [MN]{b+ [Ctiv + {R" }= [Mf )D}, Eq. 2.21 Now the effective earthquake forces are PI j= [M]ebj Eq. 2.22 To help illustrate the concept, consider the inverted Lshaped frame with lumped masses subjected to horizontal ground motion shown in Fig. 22. Assuming the elements to be axially rigid, the three DOF's are as shown. Static application ofDg= 1 results in the displacements shown in Fig. 22. Thus f={l 1 0}T in Eq. 2.21, and Eq. 2.22 becomes Pg (t)= [MiDg (t) m,+m, [0 Eq. 2.23 12=1 U 13=0 U3 SI= 1 (m2 +m3)D(t) Stationary base c) Fig. 22. Support motion of an Lshaped frame. a) Lshaped frame; b) influence vector f: static displacements due to Dg=1; c) effective load vector after Chopra (1995) Note that the mass corresponding to b, = 1 is m2+m3, because both masses will undergo the same acceleration since the connecting beam is axially rigid. The effective forces in Eq. 2.23 are shown in Fig 22. Observe that the effective forces are zero in the vertical DOF's because the ground motion is horizontal. Mass Matrices. Consistent and Lumped A mass matrix is a discrete representation of a continuous distribution of mass. A consistent element mass matrix is defined by Eqs. 2.5 that is, by [m]= Jp[N] [N]dV. It is termed "consistent" because [N] represents the same shape v functions as are used in the displacement field interpolation, and to generate the element stiffness matrix. A simpler formulation is the lumped mass matrix, which is obtained by placing particle masses mi at nodes I of an element, such that Y m, is the total element mass. Particle "lumps" have no rotary inertia unless rotary inertia is arbitrarily assigned, as is sometimes done for the rotational DOF of beams and plates. A lumped mass matrix is diagonal but a consistent mass matrix is not. The two formulations have different merits, and various considerations enter into deciding which one, or what combination of them, is best suited to a particular analysis procedure. The mass matrix for a 3D uniform beam element, which is used to model the pier and piles, and for a shell element, which is used to model the pile's cap, are developed next. Mass Matrix for the Uniform 3DBeam Element Consistent The formulation found here is given by Przemieniecki (1968). As a local coordinate system, consider the system shown in Fig.2.3. The origin is at node 1 with the ox axis taken along the length of the beam and with the oy and oz axis as the principal axes of the beam cross section. The matrix N for this element consists of twelve displacements, six deflections and six rotations, that is, Eq. 2.24 ty dd5 d3 d, .. dd7 d dd 7 d10 x Fig. 23. 3D Beam element Using the engineering theory of bending and torsion and neglecting shear deformations, we can easily show that the matrix N in the relationship {u}=[N]{d} is given by Eqs. 2.25. 1 6(4 _2) 6(4 2) (1 4 34 2)L4 (1+443 2)Lr 6( _2) 6( ~ 0 ( 24 + 32)Le (24 32)L 0 132+2(3 0 (1 )L 0 (4 245 +3)L 0 3 2 24 0 0 (2 4)L 0 0 132 +243 (1 4)Lr ( +2 2 )L 0 0 0 322 3 34 2 s LM (2 3)L 0 NT Eq. 2.25 The nondimensional parameters used in these equations are x y z SL Eq. 2.26 where L is the length of the beam. The matrix N can then be substituted into Eq. 2.5, and performed integration over the whole volume of the element. The resulting 12 x 12 consistent mass matrix is given by 13 61  +.. 35 0 _; i 13 61 0 0 I' 4 0 0 0 0 40nc o o o  o a  9 6 210 13 1. 13 6 0 * 0 0 ?0 2 35210 1L? I 2f. o o o o o 9 a ', 3 33 1 , T 1O S 3Q4TL 0 0 0 3514,9 0 0 0 0 0 0 0 S ~0  0 0 o  9 131 13 L z T lw 13 6 2 2 0 a o o 4 0 0 0 ,&TGZ 14 o 3a 0 o 0 0 w T+i 7.7 0 M I i 12o I 1, 21+ 4o t 0 0o 0 Z 210 o G 0 o o T054+ where the matrix terms with the moments of inertia ly or I, represent rotatory inertia and the terms with the polar moment of inertia J, represent torsional inertia of the element. Lumped The lumped mass matrix for the uniform beam segment of a threedimensional frame element is simply a diagonal matrix in which the coefficients corresponding to translator and torsional displacements are equal to onehalf of the total inertia of the beam segment while coefficients corresponding to flexural rotations are assumed to be zero. The diagonal lumped mass matrix for the uniform beam of distributed mass m = pA and polar mass moment l. = pJ, of inertia per unit of length may be written conveniently as m= iil 1 1 / 0 0 1 1 1 Ij/m 0 o] Eq. 2.28 Tables 21 and 22 compare the accuracy of the consistent and lumped mass formulations using finite elements. Table 21. Natural frequencies of a uniform cantilever beam: Element and exact solution ConsistentMass Finite Number of Finite Elements, Ne Mode 1 2 3 4 5 Exact 1 3.53273 3.51772 3.51637 3.51613 3.51606 3.51602 2 34.8069 22.2215 22.1069 22.0602 22.0455 22.0345 3 75.1571 62.4659 62.1749 61.9188 61.6972 4 218.138 140.671 122.657 122.320 120.902 5 264.743 228.137 203.020 199.860 Source: Chopra (1995). Table 22. Natural frequencies of a uniform cantilever beam: LumpedMass Finite Element and Exact Solution Number of Finite Elements, Ne Mode 1 2 3 4 5 Exact 1 2.44949 3.15623 3.34568 3.41804 3.45266 3.51602 2 16.2580 18.8859 20.0904 20.7335 22.0345 3 47.0294 53.2017 55.9529 61.6972 4 92.7302 104.436 120.902 5 153.017 199.860 Source: Chopra (1995). Mass Matrix for the Shell Element Because in FLPIER the pile's cap is modeled as true rectangular 9node shell elements, we will limit the formulation to this particular type of element. Consistent Consider the true rectangular 9node shell element shown in Fig. 24 below. I b b Fig. 24. True 9node rectangular element Now consider the following mapping: Fig. 25. Mapping for a true rectangular 9node shell element It is easy to verify that the shape functions N for each node are 1 N, = ( 1)(1)p' I 1 N2 4(l + 1)( 1)P11 1 N, 4(p + 1)(T + 1)LT N4 (P 1)(W + 1)pq 1 N, =  + 1)( 1)( 1)) 1 N6 = (p + 1)(pi)(i 1)(1 + 1) 1 N,  (p + 1)(p 1)(n1 + l)Tq + 1 N9 = (+1)( 1)(I + 1)(1) Eqs. 2. 29 u 08 Centerline Fig. 26. Shell element of uniform thickness If now we consider for each node the six DOF's illustrated in Fig. 26, which define the shell element, the matrix of shape functions [N], recalling again the relationship {u}= [N]{d}, takes the form it ( Ni...N9 zeros N,. N, N = NEq. 2. 30 zeros N,...N, Ni...N,_ Note that [N]'[N] has dimensions 54 x 54, which are the correct dimensions for the mass matrix considering the shell element illustrated in Fig. 26. It is easily verified that for the true rectangular element illustrated in Fig. 24 the change of coordinates from Ip and Tr to x and y can be expressed as x = aL y= br Eq. 2.31 and recalling for convenience Eq, 2.5 [m]= Ip[N]'[N]dV Eq. 2.32 If we now consider uniform the thickness t and the mass density p over the entire element, the differential volume dV can be written dV=tdxdy and Eq. 2.5 now becomes [m]= pt [N]r[N]dxdy Eq.2.33 the differencials dx and dy can obtained directly from Eq. 2.31, dx = adit dy = bdrI Eq. 2.34 and Eq. 2.32 can be rewritten [m] = ptab J [N]'[NJ]dld Eq. 2.35 1 The integral in Eq. 2.35 can be easily evaluated by means of Gaussian Quadrature (the reader is referred to Appendix B for more details on this procedure), therefore any element of the mass matrix can be obtained using Eqs. 2.30 and Eqs. 2.35 and this completes the formulation of the mass matrix for the true rectangular 9node shell element. Lumped Lumping the mass for a beam element is a process that seems to be possible by intuition and physical insight, however, for higherorder elements, like the shell element, or elements of irregular shape, intuition can be risky. Accordingly, systematic schemes for lumping are necessary. In FLPIER the HRZ scheme is used. The HRZ scheme (Cook et al., 1989) is an effective method for producing a diagonal mass matrix. It can be recommended for arbitrary elements. The idea is to use only the diagonal elements of the consistent mass matrix, but to scale then in such a way that the total mass of the element is preserved. Specifically, the procedural steps are as follows (Cook et al., 1989): 1) Compute only the diagonal coefficients of the consistent mass matrix. 2) Compute the total mass of the element, m. 3) Compute a number s by adding the diagonal coefficients mi, associated with the translational DOF (but not rotational DOF, if any) that are mutually parallel and in the same direction. 4) Scale all the diagonal coefficients by multiplying them by the ratio m/s, thus preserving the total mass of the element. Following this procedure for the true rectangular 9node shell element with uniform thickness t, results in the following lumped mass distribution 4 1 36 1 36 36 4 16 4 0o  36 36 36 1 1 36 4 36 36 Fig. 27. Lumped mass matrix at the nodes of true rectangular 9node shell element. Numbers shown are fractions of the total element mass at each node Remarks about the mass matrix Cook et al. (1989) makes important remarks about the mass matrix. The first one is the fact that the mass matrix chosen must correctly represent the mass distribution on the element, because the product [m]{d} or [M] {D) must yield the correct total force on an element according to Newton's law F = ma when {d} represents a rigidbody translational acceleration. The second remark is about the consistent and lumped mass matrix. While the consistent mass matrix is always positive definite, the same can not be said about the lumped mass matrix. If it contains zeros or negative entries in the diagonal, then it is called positive semidefinite. This may or may not cause some matrix operations to give strange results. He also suggests that a consistent mass matrix may be more suitable for flexural problems, while lumped mass matrices usually yield natural frequencies that are less than the exact values. The main reason for the development of lumped mass schemes was to save memory space in computational calculations. In the past this was a problem, but today computers are much faster and memory has become cheap and abundant. Therefore the use of consistent mass formulation is justified, since it avoids certain instabilities in the matrices operations. Damning Damping in structures is not viscous, i.e. is not proportional to velocity; rather it is due to mechanisms such as hysteresis in the material and slip in connections. These mechanisms are not yet well understood. Moreover, these mechanisms are either too difficult to incorporate into the analysis, or they make the equations computationally too expensive. Therefore with the actual limited knowledge about damping mechanisms, viscous damping is usually adopted in most analysis. Comparisons of theory and experiment show that this approach is sufficiently accurate in most cases. Damping in structures can be considered in two ways: 1)phenomenological damping methods, in which the actual physical dissipative mechanisms such as elasticplastic hysteresis loss, structural joint friction, or material microcracking are modeled. 2) spectral damping methods, in which viscous damping is introduced by means of specified fractions of critical damping (Critical damping, for which the damping ratio is = 1, marks the transition between oscillatory and nonoscillatory response). The first class of methods requires detailed models for the dissipative mechanisms and almost always result in nonlinear analyses. In the second class of methods, experimental observations of the vibratory response of structures are used to assign a fraction of critical damping as a function of frequency, or more commonly, a single damping fraction for the entire frequency range of a structure. The damping ratio 4 depends on the material properties at the stress level. For example, in steel piping 4 ranges from about 0.5% at low stress levels to about 5% at high stress levels. In bolted or riveted steel structures, and in reinforced or prestressed concrete, S has the approximate range 2% to 15% (Cook et al., 1989). A popular spectral damping scheme, called Rayleigh or proportional damping is used to form the damping matrix [C] as a linear combination of the stiffness and mass matrices of the system, that is [C]=a[K]+ 3[M] Eq. 2.36 where a and 0 are called, respectively, the stiffness and mass proportional damping constants. Matrix [C] given by Eq. 2.36 is an orthogonal damping matrix because it permits modes to be uncoupled by eigenvectors associated with the undamped eingenproblem. The relationship between a, P, and the fraction of critical damping E is given by the following equation, S (aw + Eq. 2.37 damping constants a and P are determined by choosing the fractions of critical damping (tj and 2) at two different frequencies (ol and 02) and solving simultaneous equations for a and p. Thus: a=2(22 )/(22 2) Eq. 2.38 P = 2Bo(o(,2 co,) / (o o2) Eq. 2. 39 The damping factor a applied to the stiffness matrix [K] increases with increasing frequency, whereas the damping factor P applied to the mass matrix [M] increases with decreasing frequency. For structures that may have rigidbody motion, it is important that the massproportional damping not be excessive. Usually, the natural frequencies ot and Co2 are chosen to bound the design spectrum. Therefore ol is taken as the lowest natural frequency of the structure, and )2 is the maximum frequency of interest in the loading or response. Cook et al. (1989) suggests a value of 30 Hz as the upper frequency for seismic analysis, because the spectral content of seismic design spectra are insignificant above that frequency. Estimating Modal Damping Ratios Because damping is still an unknown subjected the estimate of damping rations still presents some challenge. Recommended damping values are given in Table 23 for two levels of motion: working stress levels or stress levels no more than onehalf the yield point, and stresses at or just below the yield point. For each stress level, a range of damping values is given; the higher values of damping are to be used for ordinary structures, and the lower values are for special structures to be designed more conservatively. In addition to Table 23, recommended damping values are 3% for unreinforced masonry structures and 7% for reinforced masonry construction. It is important to note that the recommended damping ratios given in Table 23 can only be used for the linearly elastic analysis of structures with classical damping. This implies that all the material components of the structure are still behaving in their linearelastic phase. For structures subjected to strong motions, that will lead to crushing of concrete or yielding of steel, characterizing nonlinear material behavior, hysteretic damping must be added to the analysis through nonlinear forcedeformation relationships. Table 23. Recommended damping rations for structures Stress Level Type and Condition Damping Ratio(%) of structure Working stress, no more Welded steel, prestressed 23 than about /2 yield point concrete, wellreinforced concrete (only slight cracking) Reinforced concrete with 35 considerable cracking Bolted and/or riveted steel, 57 wood structures with nailed or bolted joints At or just below yield point Welded steel, prestressed 57 concrete (without complete loss in prestress) Prestressed concrete with 710 no prestress left Reinforced concrete 710 Bolted and/or riveted steel, 1015 wood structures with bolted joints Wood structures with nalied 1520 joints Source: Chopra (1995). Mass Condensation A useful tool to decrease the number of DOF's in the static analysis of a system without losing accuracy is static condensation. In this approach some degrees of freedom of the structure are chosen as master degrees of freedom and the remaining ones are called slave degrees of freedom. The choice for the master DOF is concerned with those DOF that give a better representation of the system, i.e. in a building where the beams are much stiffer than the columns the shear DOF's would be a good choice, instead of the rotations. The condensed stiffness matrix obtained includes the effects of the slave degrees of freedom, which can be recovered at any time during analysis. The same approach can be applied to the mass matrix, however dynamic condensation (Miller, (1981), Paz (1985)) is not exact, as will be shown later in this section, but it can give good results if some rules are observed in the modeling. The approach described here is given by Meirovitch (1980, pp.371372): Let us write the equations for the potential and kinetic energy for a system in the matrix form: V = [D]T[K][D] Eq. 2.40 T= [D]T[M][b] Eq. 2.41 and divide the displacement vector [D] into the master displacement vector q2 and slave displacement vector qi, or [D] =[] Eq. 2.42 Then the stiffness matrix [K] and mass matrix [M] can be partitioned accordingly, with the result V = K,, Eq. 2.43 = Lq2 L K2, Mj2 q2 Eq. 2.44 2L2 L 2L M22JLq2i Eq.2.4I where K21=K2T and M21=M1T. The condition that there be no applied forces in the direction of the slave displacements can be written symbolically in the form aV q=q,K,, +q2K21 Eq.2.45 aq, where the equation implies equilibrium in the direction of the slave displacements. Solving Eq. 2.45 for qI, we obtain q, = K1'Kl2q2 Eq. 2.46 which can be used to eliminate q, from the problem formulation. Equation 2.46 can be regarded as a constraint equation, so that the complete displacement vector q can be expressed in terms of the master vector q2 in the form q = Pq, Eq. 2.47 where P is a rectangular constraint matrix having the form P [ ,] Eq. 2.48 in which l is a unit matrix of the same order as the dimension of q2. Introducing Eq. 2.47 into Eqs. 2.40 and 2.41, we obtain V = q2 qK2q2 T2 where the reduced stiffness and mass matrices are simply K, = P'KP = K22 K,2K 'K, Eq. 2. 49 Eq. 2. 50 Eq. 2.51 M, = PMP= M, KK' M2 M,2K1'K21+KK,' M,,K' K, Eq.2.52 The matrix MI is generally known as the condensed mass matrix. What is being sacrificed as a result of the condensation process? To answer this question, let us consider the complete eigenvalue problem, which can be separated into K,,q1 + K12q2 = (M1q1 + M12q2) K21q1 + K22q2 = (M21q, + M22q2) Solving Eq. 2.54 for q2, we have q, = (M2 K21)(K22 M22)q2 Eq. 2. 53 Eq. 2.54 Eq. 2.55 so that, introducing Eq. 2.55 into Eq. 2.53, we obtain (K,, K,, K2' K2,)q = x(M,, K,, K,' M,, M2 K,'K2, +K KK,' M,, K,'K,,)q, + L2(M,2K, 'M,,2 K,2K,'M,,K,'K,, K,, K, M,2K,2 M2, + K,2KM'M2K,2K21)q, + Eq. 2. 56 Examining Eqs. 2.56 and 2.52 we conclude that the condensation used earlier is ignoring second and higher order terms in X in Eq. 2.56, which can be justified if the coefficients of X, X3, ... are significantly smaller than the coefficient of X. For this to be true we must have the entries of M12 and M22 much smaller than the entries of K12 and K22. Physically, this implies that the slave displacements should be chosen from areas of high stiffness and low mass. Moreover the nodes that carry a timevarying load should be retained as master. FLIPIER condenses the stiffness and mass to the top of the piles. While this procedure is exact for the stiffness, note that it is not for the mass. This is because in pier structures the slave DOF's are located in areas of high mass concentration, like the pile's cap. Also note that the mass of the superstructure on the pier, is usually modeled as lumped masses at the top of the pier. In the case of earthquake loading for example note that we have the loading function acting on slave DOF's, what is not acceptable in this approach. Because of these limitations for mass condensation, this approach is not recommended for the dynamic analysis of bridge piers. FLPIER was then modified to allow a "full" analysis of the structure, where neither the stiffness nor the mass matrices are condensed to the top of the piles. A typical full and the respective condensed version of a structure are shown in Fig. 28. In a typical condensed static analysis, the stiffness and loads of the superstructure are condensed to the top of the piles. The condensed analysis is then carried out. At the end of the analysis the superstructure's forces and displacements are recovered and the analysis is terminated. Structure S Cap Piles Condensed masses Condensed Piles Fig. 2.8. Full and condensed versions of the structure TimeHistory Analysis. Direct Integration Methods In direct integration methods or stepbystep methods, a finite difference approximation is used to replace the time derivatives appearing in Eq. 2.12 or 2.13 at various instants of time. Direct integration is an alternative to modal analysis methods. For many structural dynamics and wave propagation problems, including those with complicated nonlinearities, direct integration is easier to implement. In direct integration, the approach is to write the equation of motion (2.12), at a specific instant in time, [M]({)} +[C]({b} +[K]{ D}, = (R" Eq. 2. 57 where the subscript n denotes time nAt and At is the size of the time increment or time step. The absence of time step subscripts on matrices [M], [C], and [K] in Eq. 2.57 implies linearity. For problems with material or geometric nonlinearity, [K] is a function of displacement and therefore of time as well. Accordingly, from Eq. 2.57, [M]j{D} +[C]{D}) + R"}) = {R"'} Eq. 2.58 {R'"',} is the internal force vector at time n At due to straining of material. It is obtained by assembling element internal force vectors, {ri't},, given by Eq. 2.7 using {o}n. For nonlinear problems, {R""}, is a nonlinear function of {D}n and possibly time derivatives of {D)}, if the strain rate is an issue. For linear problems, the internal force vector is given by the relationship {R"in,=[K]{D},. In Eq. 2.58 [M] and [C] are taken as time independent, although for some problems these may also be nonlinear. In this work [M] and [C] remain constant during the analysis, and the internal force vector, {R'"'}, is only a function of the displacements {D},. The nonlinear relationship for {o}, will be developed later. Different methods for direct integration of Eqs. 2.57 and 2.58 can be categorized as explicit or implicit. The first category, the explicit methods, has the form { D} = f({ D),,O {/(}, { D}_,,...) Eq. 2. 59 and hence permit {D}+,, to be determined in terms of completely historical information consisting of displacements and time derivatives of displacements at time n At and before. The main characteristic of explicit methods is the fact that the next approximation for the displacements is based only on the known previous approximations for the displacements and their respective derivatives. Note the use of the equilibrium condition at time n. The Central Difference is an example of an explicit method. The second category, the implicit methods, has the form {D) = f( { { D,.... )Eq.2.60 and hence computation of {D},+1 requires knowledge of the time derivatives of {D},+1, which are unknown. The main characteristic of the implicit methods is the fact that the next approximation for the displacements depends on unknown values of their derivatives. Note that the equilibrium condition is used at time step n+l. Newmark's and the WilsonTheta are examples of implicit methods. There is vast literature about the advantages and disadvantages of each approach, the reader is referred to Cook et al. (1989), Chopra (1995), Paz (1985), Craig (1981), or Crandall (1956), for a more extensive discussion on these approaches. Generally speaking under certain conditions the implicit methods are more stable than the explicit methods. Because an implicit method was used in this research we will limit our discussion to this class of methods. Numerical Evaluaton of Dynamic Response. Newmark's Method As mentioned earlier implicit methods are those in which the approximation for the next displacements {D},,+ depends on unknown values for its time derivatives. The main advantage of the implicit methods is the fact that most of the useful methods are unconditionally stable and have no restriction on the time step size other than required for accuracy. In 1959 N. M. Newmark developed a family of timestepping methods based assumptions for the variation of the acceleration over the time step. The first method is 37 called average acceleration and is shown in Fig. 29. Successive application of the Trapezoidal Rule leads to Eqs. 2.61 to 2.65. ui ti ti+l t At Fig. 29. Average acceleration ii() =0(",+i+) Eq. 2.61 At = "u + 2 (, + ki,) Eq. 2. 63 2 u(T) = u, + iiu + (ii, + ,) Eq. 2. 64 At2 u,,I =u, + iAt + (ii, + ii) Eq. 2.65 The second method is called linear acceleration and is illustrated in Fig.210. The Trapezoidal Rule is also used to obtain Eqs. 2.66 to 2.70. U A u,1 * t t,+1 t At Fig. 210. Linear acceleration ii() = ii, + (ii,, + ii ) Eq. 2. 66 At 12 u(T)= +iT + (ii, +u) Eq. 2.67 2At At ui1 = u + 2(i,, + +i,) Eq. 2. 68 2 3 ,C2 C u(t)= u, +th+iu,+ (u, u,) Eq. 2.69 2 6At u,, =u, +z,At+(At)2 ~i, + ij Eq. 2.70 The Newmark's family of methods can be summarized into the following two equations: ",+, = u, + [( )iii +yii,, ]At Eq. 2.71 u,W = u, + uAt + p ii, + pii/, At2 Eqs. 2. 72 where y and p are parameters that can be determined to obtain integration accuracy and stability. Note that when y = V2 and P = 1/6, this method reduces to the linear acceleration method. Newmark originally proposed as an unconditionally stable scheme, the constant average acceleration method, in which case y = /2 and P = '. This can be shown considering that Newmark's method is stable if (Chopra, 1995): At 1 1 < Eq. 2. 73 T. n 5 2 P For y = /2 and 1 = /4 this condition becomes At  < o Eq. 2. 74 This implies that the average acceleration method is stable for any At, no matter how large, as mentioned before; however, it is accurate only if At is small enough. For y = / and p =1/6, Eq. 2.73 indicates that the linear acceleration method is stable if At <: 0.551 Eq. 2. 75 T For the solution of displacements, velocities and acceleration at time i+1 we now consider the equilibrium equation also at time i+1: Miii, + Ci,,, + Ku,1 = F+,, Eq. 2. 76 Solving Eq. 2.72 for uii, in terms of ui+ and substituting in Eq. 2.71, equations for u,,+,and u,,,in terms of the unknown displacements u,+; only are obtained. Substituting these equations into Eq. 2.76, a system of equations is obtained, which can be solved to obtain u+IA: (boM + b,C + K)u,, = F,, +M(bou, +b2i, + bii,)+C(b,u, +b4ii, +bii)Eq. 2. 77 where b = ;b, = ;Ab = b =b 1;b4 = 1;b5 =At f2/2Eq.2.78 S Atr2 PAt t 20e 04 and finally all the quantities at time i+1 can be written as iii+l = b (u,,, u,)+ b7u, + b8ii, Eq. 2. 79 ui, = u, + bii, + boii, Eq. 2.80 ,i = u, + (u,, u) Eq. 2.81 where be = bo; b7 = b2; bs = b3 ; b9 = At(y) ; blo = y.At Eq. 2.82 Equation 2.77 can be now written in the condensed form: Ku,,, = p(t) Eq. 2. 83 where the effective stiffness K is given by K=boM+b,C+K Eq. 2.84 and the effective load vector P(t) is p(t) = F,, +M(bou, +b26i, +b3ii,)+C(bu, +b4u, + bsii,) Eq. 2.85 which completes the formulation for linear systems. Choice of time step At The unconditional stability of the average acceleration method may lead some analysts to adopt larger time steps because of economic needs, implying a solution that may not be accurate, because unconditional stability does not mean unconditional accuracy. Cook et al. (1989) suggests the following expression for the time step At: At <(27t/,.)/20 0.3/o, Eq. 2.86 where (o, is the highest frequency of interest in the loading or response of the structure. However, he suggests that in the case of convergence difficulties, the analysis should be repeated with a smaller time step for additional assurance of a correct solution. Nonlinear Problems In structural analysis, a problem is nonlinear if the stiffness matrix or the load vector depends on the displacements. The nonlinearity in structures can be classified as material nonlinearity (associated with changes in material properties, as in plasticity) or as geometric nonlinearity (associated with changes in configuration, as in large deflections of a slender elastic beam). In general, for a timeindependent problem symbolized as [K]{D}={R}, in linear analysis both [K] and {R} are regarded as independent of {D}, whereas in nonlinear analysis {K} and/or {R} are regarded as functions of {D}. In this dissertation we will address the nonlinearities associated with changes in material properties. Analysis of the Nonlinear Response using Newmark's Method In this section the notation for SDOF systems is used to simplify the approach, the extension to MDOF systems is done later, although the concept is exactly the same. In the numerical evaluation of the response of a dynamic system we go from time step i, where the equation of motion can be written mii + cu, + (s), = p, Eq. 2.87 to time step i +1, where the dynamic equilibrium can be written mii,, + cli,, + (fs), = p,1 Eq. 2. 88 where (fs); is the system's resisting force at time i. For a linear system (fs) = kui. However for a nonlinear system the resisiting forcefs would depend on the prior history of displacements and velocities. It is then necessary to consider the incremental equilibrium equation, the difference between Eqs. 2.84 and 2.85, can be written: mAii, + cAu, + (Afs), = Ap, Eq. 2. 89 The incremental resisting force can be written (Afs), = (k,) Au Eq. 2. 90 where the secant stiffness (kj)sec, shown in Fig. 211, cannot be determined because ui,+ is not known. If however we make the assumption that over a small time step At the secant stiffness (k,)se, can be replaced by the tangent stiffness (k,)an, then Eq. 2.90 can be rewritten: (Afs), = (ki) Au, Eq. 2. 91 The incremental dynamic equilibrium equation is now: mAii, + cAi, + (k,), Au, = Ap, Eq. 2. 92 Equation 2.92 suggests that the analysis of nonlinear systems can be done by simply replacing the stiffness matrix k by the tangent stiffness (k,)rto be evaluated at the beginning of each time step. However this procedure for constant time steps At can lead to unacceptable results for two reasons: a) The tangent stiffness was used instead of the secant stiffness. b) The use of a constant time step delays detection of the transitions in the force deformation relationship. I U, U,I, u Fig. 211. Secant and Tangent approaches. After Chopra(1995) These errors can be minimized by using an iterative procedure within each time step. The idea is to guarantee dynamic equilibrium before going to the next time step. We must then solve the equation ku,,, = h(t) Eq. 2. 93 where now the effective stiffness k, becomes k, = k, + bm + bc Eq. 2.94 For convenience we drop the subscript i in ki and replace it by T to emphasize that this is the tangent stiffness, and we can rewrite Eq. 2.94 as k, = kT + bm+b c Eq. 2. 95 The effective force P(t) is now given by k(t) = f(t),, + m(b2z, + bii) + c(b4z, + b,u,,,) r, Eq. 2.96 where r, is the internal force at time step i. Since we are using the tangent estimate for the stiffness, we must iterate within each time step to guarantee dynamic equilibrium. The first step is to apply the effective force P(t) and get the first approximation for the next displacement uki+1. Associated with this displacement is the new tangent stiffness kkr and the true force f, which includes the internal force rk, the inertia force/;, and the damping force/, at iteration k. An outof balance force defined as A p(t) k+ = (t) k / is then generated. The new effective stiffness for the system is computed and the system of equations kAuk = APk is solved to obtain the incremental displacements Auk, which is added to the previous displacement u,, giving a new estimate for the displacements, velocities and accelerations at time i+1. The dynamic equilibrium is checked again with the new values and another outof balance force is generated. The process continues until the outofbalance force is within a specified tolerance. Figure 212 helps to understand the process. This process is known as the NewtonRaphson method. If the reader is familiar with numerical methods, will notice that what is defined as the The NewtonRaphson method in some references (Cook et al. (1989), Chopra (1995) and Bathe (1996)) and above is in fact a NewtonRaphson approach. Although irrelevant to this discussion the author found important to define the NewtonRaphson method from a mathematical point of view. The NewtonRaphson method is a mathematical method for finding roots. The method relies on the fact that the functionf(x), and its first derivative f(x), and second derivative f '(x) are continuous near a root p. Note how rigorous the mathematical definition for the method is, since this was not assumed in the previous application. With this information it is possible to develop algorithms that produce sequences {xk) that converge faster top. Note however thatf(x) and its first and second derivatives must be continuous for the method to work. The convergence rate of the NewtonRaphson method is quadratic. This method is also called the tangent method because the slope of the curve is used in the formulation. The recursive formula for the method is: T(xk,)=k f(xk ,for k = 1,2,... f'(x ) Eq. 2. 97 and the sequence will converge top. I U3 Fig. 212. NewtonRaphson Method However sometimes it may be difficult if not impossible to implement Newton's method if the first derivative is complicated. The secant method may then be used. The secant method is the same as the Newton's method except that the first derivative is replaced by its approximation: the slope of the line through the two previous points. Note that when successive points get close, the method becomes unstable because division by zero may occur. It can be shown that the convergence rate is approximately 1.62, being slower than Newton's method. The recursive formula for the secant method is shown in Eq. 2.97. Note the clear connection between the mathematical and engineering terms used to define both approaches. More detailed information about both methods can be found in Mathews (1987). T(x.,,)= x. x x ( x ) Eq. 2. 98 "f' (x) f(x,) f(x,, ) f(x,) f(x,, ) Xn X,_1 The extension to MDOF systems is immediate by replacing all the scalar quantities by its respective vector equivalents. Note however that for a SDOF system the tangent stiffness kr is easily obtained. Such a simple evaluation is not available if there are MDOF. However, in practice, as it will be discussed later, the physics of the problem allows us to calculate the tangentstiffness matrix [K,]. In a MDOF context NR iteration involves repeated solution of the equations [Kt]i{AD}i+i={AR},+i, where the tangent stiffness matrix [K,] and load imbalance {AR} are updated after each cycle. The solution process seeks to reduce the load imbalance, and consequently {AD}, to zero. The internal force in the equation of motion, Eq. 2.35, is written as R" i ,= R"}+[K,]{AD} Eq. 2.99 where {AD} = {D},,, {D}, Eq. 2. 100 Combining Eqs. 2.99 and 2.100 with the equations of motion, Eq. 2.35, and the trapezoidal rule, we obtain [K" {AD) }={IR Eq. 2. 101 where K = 4 [M+ 2 [C]+[K,] Eq. 2.102 At2 At and R. = {R "}.+ {RM + [M (4 {D} + (D,) + [C{, Eq. 2.103 Note that [K,] must be predicted using {D}n (and possibly {)} if strain rates effects are important) and must be factored at least once each time step during nonlinear response. If [K,] is not an accurate prediction of the true tangentstiffness matrix from time nAt to time (n + 1)At, then the solution of Eq. 2.101 for {AD} will be in error. The error in nodal forces, that is the residual, is given by the imbalance in the equation of motion as {R~"= {R",,,, [M]{D ,, [C]jd,, {R'n }., Eq. 2. 104 where the internal force vector {R'"'}n,+ is computed elementby element. The process stops when the outofbalance force vector {(R"} is smaller than a specified tolerance. Finally, considering all the topics discussed in this section, the procedure for nonlinear dynamic analysis can be summarized in the algorithm described next. Nonlinear Dynamic Analysis Algorithm Step (1). Form the initial stiffness matrix K, mass matrix M, and damping C for the structure. Step (2). Compute the constants for the numerical method chosen: Wilson0 Newmark 6 3 1 6 bo= 6 ;b, 3 bo= ;b, (OAt)2 'OA aAt2' aAt b2 = 2b,;b = 2 b2 =l/aAt;b3 =1/2a 1 b =2;b, =6At/2 b4 =l;b5 =At(2)/2 b6 = b e /0;b, = b/0 b6 = bo;b, = b2 b, =13/0;b9 =At/2 b = b3;b9 = At(1) 6 Ab2/6 = 5At For the Newmark's Method Average Acceleration: a = and 5 = '. For the Newmark's Method Linear Acceleration: a = 1/6 and 6 = '. For the Wilson0 Method: 0 = 1.4 (0 = 1.0 for Newmark's). Step (3). Initialize iio0, and uo. Step (4). Form the effective stiffness matrix using the initial stiffness matrix K: K = b,M +b,C +K Step (5). Beginning of time step loop. Step (6). Form the effective load vector for the current time step: F,, = F, +0(F,+, F,)+ M(b2i, +b3ii,)+C(b4u, + bii,) R, Step (7). Solve for the displacement increments 6u = I^l@ / Step (8). Beginning of dynamic equilibrium loop. Step (9). Evaluate approximations for ii,u and u: ii'l = bo0u'' b22, b3,i ul'' = bu'0' b4 u, b,ii, u i = u,,8u,' Ilt+OAI Step (10). Evaluate actual tangent stiffness K' and internal forces R'K for all the elements in the structure. Step (11). Evaluate new effective stiffness matrix K, based on actual tangent stiffness K'. Step (12). Evaluate the outofbalance dynamic forces: s = F, +e(F, F,) [Mii +Ciu +R'' Step (13). Evaluate the ith corrected displacement increment: Step (14). Evaluate the corrected displacement increments: 6ui =6u''i +Au' Step (15). Check for the convergence of the iteration process: NO 'NO T, < (tolerance) Return to Step (8) YES Step (16). Return to Step (5) to process the next time step. In this dissertation the full NewtonRaphson Method is applied. In this approach the tangent stiffness matrix k, is updated for every iteration, in contrast with the Modified NewtonRaphson Method in which the tangent stiffness k, is update once every time step (it remains constant during the iterations performed within each time step). This improves convergence, but additional computational effort is required in forming a new tangent stiffness matrix k, and factorizing it at each iteration cycle. The 50 iterations within each time step proceed until the value of the outofbalance force ,', is smaller than a specified tolerance T. The mass matrix Mand damping matrix C are considered to be constant throughout the analysis. As mentioned earlier in this Chapter the correct choice of time step, which depends on the lowest natural frequency of interest for the system, enforces stability and accuracy to the analysis. Note however that for a nonlinear analysis, the case is that due to yielding of steel or crushing and cracking of concrete, a lower stiffness will be obtained during the analysis, which will result in even lower values for the lowest natural frequency of interest. The reader should have this in mind when choosing an adequate time step for a nonlinear dynamic analysis. CHAPTER 3 DISCRETE ELEMENT MODEL AND MATERIAL HYSTERESIS Discrete Element Derivation A representation of the discreteelement used in FLPIER is shown in Fig. 31. It was developed by Mitchell (1973) and modified by Andrade (1994). The center bar can both twist and extent but is otherwise rigid. The center bar is connected by two universal joints to two rigid and blocks. The universal joints permit bending at the quarter points about the y and z axes. Discrete deformational angle changes ?1, P2, 3j and Tf4 occur corresponding to the bending moments Mi, M2, M3 and M4, respectively. A discrete axial shortening (6) corresponds to the axial thrust (T), and the torsional angle T's corresponds to the torsional moments in the center bar Ms. ^, *!*t_ Universal Join (Top View) Spring (Side View) Rigid end Block Fig. 31. Representation of discrete element. After Hoit et al., 1996 Element Deformation Relations In Fig. 32, wlw3 and w7w9 represent the displacements in the x, y and z directions at the left and right ends, respectively. The displacements w4 and w1o represent axial twists (twists about the xaxis) at the left and right ends, respectively; and finally w5w6 and w1wu1 represent the angles at the left and right end blocks about the x and z axes, respectively. Based on a small displacement theory we can write: h n = w, w, (w5 + w,1) Eq.3.1 h s = w, w (2 + w,,) Eq. 3.2 The elongation of the center section of the element is calculated as follows: 8 = w, w, Eq. 3.3 The angle changes for the center section about the z and y axes are the defined as s ws w2 w6 +w12 0,2 612 Eq.3.4 h h 2 s wzw w +w, 02 51 Eq. 3. 5 h h 2 The discretized vertical and horizontal angle changes at the two universal joints are then ', = 0 w6 Eq.3.6 Y2 = w, 02 Eq. 3.7 3 = wl2 01 Eq. 3. 8 Y4 =021 WII Eq. 3.9 and the twist in the center part of the element is defined as 5 = wI0 w4 Eq. 3. 10 Therefore, the internal deformations of the discrete element model are uniquely defined for any combination of element end displacements. The curvature for small displacements at the left and right universal joints about the y and z axes are defined as follows: At the left joint =4, = /h =, 2 = / h At the right joint T3 = /h S= h h The axial strain at the center of the section is given by c 2h Eq. 3. 11 Eq. 3.12 Eq. 3.13 Eq. 3.14 Eq. 3.15 Side View {3 ^ '4' f. *4 Top Vi.ew  1 I 2 2 IS Fig. 32. Discrete element displacements. After Hoit et al. 1996 Integration of Stresses for Nonlinear Materials For a beam subjected to both bending and axial loads, it is assumed that the strains vary linearly over the area of the cross section. This assumption enables the strain components due to bending about the z and yaxes, and the axial strain to be combined using super position. Examples of these three components are represented separately in Fig. 33(ac) and combined in Fig. 33(d) also shown in Fig. 33(d) is a differential force, dF,, acting on a differential area, dAi, relationship for the material dF = a,dA, Eq. 3.16 Finally, Fig. 33(e) represents the stressstrain curve. Then, for the left joint, the relationship for the strain at any point in the cross section is End View f, 1 '12 S=E, ly, l,2Z Then to satisfy equilibrium M = fdFY=j a,YdA, A A M, = \\dF,.Z, = Jfc,Z,d4, T = d =J IadA, AS A a) Strain due to b) Strain due to zaxis bending yaxis banding ) Strssstrain relationship 9) Stressstrain relationship Eq.3.17 Eq. 3.18 Eq.3.19 Eq. 3.20 c) Strain due to axial thrust d) Combined strains Fig. 33. Various components of total strain in the section. After Hoit et al., 1996 Numerical integration of Eqs. 3.18, 3.19 and 3.20 is done using Gaussian Quadrature. To use the method of Gaussian Quadrature, the function being integrated must be evaluated at those points specified by the position factor. These factors are then multiplied by the appropriate weighting factors and the products accumulated. Fig. 34 shows a square section with 25 integration points (a 5 x 5 mesh). The actual number of defaults integration points for a square section is set at 81 (a 9 x 9 mesh). For a steel H section the default number of points is 60. For circular sections, the section is divided into circular sectors (12 radial divisions and five circumferencial divisions as shown in Fig. 35), totaling 60 points. The sections are integrated at the centroids of each sector using weighting factors of 1.0. The stress in all steel bars is evaluated at the centroid and a weighting factor of 1 is used for each bar. When a circular void is encountered in a square section, the force is first computed on the unvoided section and then the force that would be acting on the voided circular area is computed and subtracted from the force computed for the unvoided section. Circular sections with voids are divided into sectors omitting the voided portion. This method of dividing the sections into points, getting the strains and stresses at each point and then integrating to get forces is usually called fiber modeling. S S S x X S. i S Concrete tegrtion Points S Steel Rebar (xl integration) Fig. 34. Rectangular section with integration points. After Hoit et al., 1996  i Steelbar(Ixl Integration) \'\ SConcrete subdivided into sectors  Note: Integration points (lxl) for concrete are at the geometric centrods of each sector Fig. 35. Circular section with integration points. After Hoit et al, 1996 Even for a nonlinear material analysis, the torsional moment M is assumed to be a linear function of the angle of twist T, and the torsional stiffness GJ, where J is the cross section torsional constant; and G is the material shear modulus, resulting in the following expression for the torsional moment M = GJT / 2h Eq. 3.21 In this discrete approach the curvature is evaluated at two deformational joints inside the element. These points are located at the quarter points from each node along the element length. Therefore the effective position of any plastic hinge that might form in the structure is restricted to these two locations. This should not cause any practical limitation for most problems. However, it should be considered if trying to match a theoretical solution with pure plastic hinges at theoretical ends of members. These discrete joints, at which the deformations are concentrated, correspond to the integration points along the length of the element, if a conventional finiteelement solution was being made. Element End Forces The element internal forces are necessary to assemble the global internal force vector necessary for equilibrium of the equation of motion. From equilibrium of the center bar V, =(M4 M)l/h Eq. 3.22 V, =(M M3)/h Eq. 3.23 and from equilibrium of the end bars fi = T Eq. 3.24 f2 = V Eq. 3.25 3 = V2 Eq. 3.26 f4 =M Eq. 3.27 f,= M +V,2 h/2+Th/2.w5 Eq. 3.28 f =M2 +V, h/2+Th/2.w6 Eq.3.29 f7 = T Eq. 3.30 fs = V Eq. 3. 31 fA = V2 Eq. 3.32 fo = M5 Eq. 3. 33 f, =M3 +V2 h/2+T.h/2 w,, Eq. 3.34 f =M4 +V, h/2+T.h/2w12 Eq.3.35 wherefi f3 andf7 f9 are the acting end forces; andf4 f6 andfiofi2 are the acting end moments. Element Stiffness As mentioned in Chapter 2 the element stiffness may change for either nonlinear static or dynamic analysis. Therefore its is necessary to evaluate the element stiffness for each iteration in each timestep. The procedure adopted uses the standard definition of the stiffness matrix; for an element having n DOF the stiffness matrix is a square matrix [K] of dimensions n x n, in which Ki is the force necessary in the ith DOF to produce a unit deflection of thejth DOF. The stiffness computed is that obtained by one of the two methods described below. The transformation of the discrete element stiffness matrix to global coordinates and the assembly of the different components of the global stiffness matrix follow standard direct stiffness procedures. Secant and Tangent Stiffness of the Discrete Element During the iteration process, the element stiffness matrix is reevaluated in each new deformed position. For each iteration, the stiffness for each integration point along the cross section within an element is stored. Then, on 12 subsequent passes, a unit displacement is applied to each element DOF keeping all other DOF fixed and the forces corresponding to that unit displacement are calculated over the cross section of the element as described earlier. If the stiffness for each of the integration points is defined by dividing the present stress by the present strain as shown in Fig. 36, then it is called secant stiffness. On the other hand, if the stiffness for each segment is defined by the slope of the stressstrain curve at the specific point being integrated, then it is called the tangent stiffness as also illustrated in Fig.36. Note that if the tangent stiffness is used, the solution of the nonlinear problems becomes NewtonRaphson, while if the secant stiffness is used we have a secant method solution. Stress / secanA Strain Fig. 36. Secant and tangent material stiffness Hysteresis Models We have shown so far how it is possible to get the element internal forces and form the updated stiffness for each iteration in the previous paragraphs. But we have considered only the case in which the element is subjected to loading. For dynamic analysis usually the case is that the element will be subjected to two additional phases: unloading and reloading. It is important to notice that even for nonlinear static analysis these two new phases could be present due to a nonlinear redistribution of forces on the structure that could cause some elements to be overloaded and others to be under loaded. The curve that is used to describe this behavior is called a hysteresis. For the inelastic analysis, a proper selection of hysteretic models for the materials is one of the critical factors in successfully predicting the dynamic response under strong motion. Several models have been proposed in the past for reproducing various aspects of reinforced concrete behavior under inelastic loading reversals. In order to closely reproduce the hysteretic behavior of various components, a highly versatile model is required in which several significant aspects of hysteretic loops can be included, i.e., stiffness degradation, strength deterioration, pinching behavior and the variability of hysteresis loop areas at different deformation levels under repeated loading reversals. However, the model should also be as simple as possible since a large number of inelastic spring are necessary in modeling the entire structure, and additional parameters to describe a complicated hysteresis loop shape may sometimes require excessive amount of information. (d) Aq () r s0ra, (9) Ct. (h) W (I) * W (n) k Fig. 37. Models for hysteresis loops proposed by some authors. After Mo, 1994 Some of the existing popular models: Clough (1966), Fukada(1969), Ayoama (1971), Kustu (1975), Tani (1973), Takeda (1970), Park (1984), Iwan (1973), Takayanagi (1977), Muto (1973), Atalay (1975), Nakata (1978), Blakeley (1973), and Mo (1988) are shown in Fig.37. It appears that most of the available models are aimed at a particular type of component, such as for use of beams, columns or shear walls only, and therefore, fall short of the versatility required for modeling practical buildings having a large number of different components. Most of these models were obtained by performing curve fits to experimental data. This approach results in really good approximations for the behavior of the specific member being studied, but lacks versatility when applied to a real structure. The advantage of the discrete element with fiber modeling is that the elements general behavior will be governed by its material properties, instead of experimental observations, making this procedure useful for studying multiple cross section configurations under more general load histories. Material Models Rather than trying to develop a new element that would model this behavior using elasticity and plasticity theories, it was thought to change the discrete element used in FLPIER to develop such behavior. Based on various studies and experiments (Chen (1982), Park and Paulay (1975), Mo (1994), Roufaiel and Meyer (1987), Park et al. (1972), Agrawal et al. (1965), Sinha et al. (1964a), Sinha et al. (1964b), Kwak (1997), Magdy and Mayer (1985), Soroushian et al. (1986), Nilsson (1979), Ozcebe (1989), Peinzien (1960), and Tseng (1976)) about the uniaxial cyclic behavior of concrete and steel, and the behavior of R/C members, the stressstrain curves of these materials were modified to work for dynamic analysis. Shear deformations were not included because it was not included in the works cited. Moreover it was also not included in the formulation of the discrete element. By using this approach it was possible to overcome some of the difficulties found in the various hysteresis models. Although computationally more expensive, this fiber modeling approach gave the modified discrete element the versatility not found in many models. A description of these modified curves follows next. Uniaxial Mild Steel Model Because the cyclic behavior of steel is very dependent on the heat from which it was produced, it was decided that rather than trying to predict the steel behavior based on exponential curves (Agrawal et al. (1965), Shen and Dong (1997)), to use the bilinear representation for the steel behavior. This represents a safe lower bound solution that is adequate for most of the construction steel. So the mild steel reinforcement is assumed to be perfectly elasticplastic (no hardening) and similar in both tension and compression as shown in Fig. 38. The parameters needed for the mild steel bilinear model are the modulus of elasticity E, and yield stress f. The rules for this model are as follows (refer to Fig. 38): Loading is represented by segment ab, the tangent stiffness is the initial modulus of elasticity for steel E, and the stress is given by = E,.E . Yielding: If E > e, yielding occurs and is represented by segment bc. The tangent stiffness E, is equal to 0 and the stress a = fy, the yielding stress for steel. The residual strain e, is given by e, = e c . Unloading is represented by segment cd, the tangent stiffness is E, and the stress is given by a = (E s,).E,. 64 Reloading is represented by segment ef, the tangent stiffness is E, and the stress is given by a = (E e,).E,. For any of these phases the secant stiffness Es, is given by E, = Eq. 3.36 C b c f r e Fig. 38. Elasticperfectly plastic model for mild steel Uniaxial Monotonic Concrete Model Used in FLPIER The concrete model used in FLPIER is generated based on the values of the concrete strength fc and modulus of elasticity of concrete Ec, input by the user. The compression portion of the curve, which is is based on the work of Wang and Reese (1993), is highly nonlinear and has a maximum compressive stress f which is related to but not always equal to the compressive strength of a standard test cylinder, f'. Based on experimental research, f, is taken to be 85% of fc, the maximum cylinder compression stress. The tension side of the curve is based on the tension stiffening model proposed by Mitchell (1973). This procedure assumes an average tensile stressstrain curve for concrete. The stress strain relation of concrete in tension is very close to linear with cracking occurring at a small rupture stress f,, The high stresses actually experienced at tensile cracks in the concrete will not be reproduced by the model. However the average response over a finite length of beam will be adequately represented. Based on the user input the program will generate the concrete curve as a series of points connected by straight lines as shown in Fig. 39. Values in between the points are obtained by interpolation of the extreme points in the interval. For values offc = 6 ksi and Ec = 4615 ksi the strain and stress values for the concrete stressstrain curve generated by the program can be seen in Figs. 310 and 311 respectively. / /  13r2,,S Fig. 39. FLPIER concrete points StressShaSbn for Concr le Sress 0.16E03C V~arl~o lil _I /0 /0.6 09901 03.00E 02 0.132iEC S0C.1650 02 'f^ fe02 33 OE03 1E03 03 IL3.n 25/ 4/1999 Fig. 310. Concrete strains StressS!to fr^ Conorete Siress 0 5009 S 163 It3.; 25/ 4/1999 Fig. 311. Concrete stresses  000,I 0010 S 0004  n.ZCc DOne iolve ~ I r~r/~ 0 Proposed Models for the Uniaxial Inelastic Cyclic Behavior of Concrete Two concrete models were implemented in FLPIER. One, called the rational model, is based on the work of Sinha et al. (1964a), the other is a bilinear model, similar to the mild steel model presented earlier. These two models are discussed in more details next. Rational Model Figure 39 shows the default envelop of the stress strain curve supplied by the program, which is a function of f' and E, input by the user. This curve is the backbone of the model, because it limits the value of the stress in concrete during all phases in the analysis. The compression portion of the concrete curve is highly nonlinear and is defined by the Hognstead parabola up to a stress equal to 0.85fc. Beyond this point, a straight line is adopted, connecting this maximum to another point with residual stress of 0.20fc and strain equal to 4 c,, as shown in the Fig. 312. For strains greater than 4 ,,, the stress is kept at a constant 0.20fc (residual stress) as suggested by Chen (1982). For the tension portion the curve is assumed linear up to a stress of f, and then has a tensionsoftening portion as shown in Fig. 312. The tensionsoftening portion attempts to account for the untracked sections between cracks where the concrete still carries some stress. The value of f, is based on the fixed value of e, shown in the figure and the modulus of elasticity E,, input by the user. For English units this will give a value of f, = 75J5 f. The rules for the rational model are described next. 0.30 f7" fA" Fig. 312. Envelop curve for concrete Loading Compression Compression follows the curve described above, a typical loading in compression phase is illustrated in Fig. 313. The equations that define this phase are: When =f =24 4{I Eq.3.37 The tangent modulus of elasticity is defined as: If E, = E, Eq. 3.38 If >e>e then Ec= f  Eq.3.39 This correction was necessary because the derivative of the Hongnestead parabola gives very high values of the tangent modulus for low values of the compressive strain e, which caused instabilities to the model. Based on research (Chen (1982), Park and Paulay (1975)) it was found that at about 30% off'c it is reasonable to assume that concrete still have the initial tangent modulus, Ec,. Equation 3.39 is just the derivative of Eq. 3.37 with respect to E. When the strain e > co, the concrete enters a phase called softening. In this research this phase is defined as: If c < c < 4 e, then 0.8 f," 0.8f:" E, = 0 Eq. 3.41 (4e, s,) and if a >4 E, f= 0.20f" Eq. 3. 42 and E, =0 Eq. 3.43 stress (ksi) 2 2.5 0.002 0.0018 0.0012 0.001 0.0008 strain in/in 0.0002 0 Fig. 313. Typical compression loading Tension A typical loading in tension is shown in Fig. 314. Tension is also represented by the envelope curve described earlier. The following equations define loading in tension If 6 < , f =E Eq. 3. 44 E, = E, Eq. 3.45 If c, < e < saf then f =fr +( E 03 ) E Y ii Eq. 3. 46 S f, Ec 6,  W If Ef < s < en, then ( 05f,  S0.5f, and finally ife > 6. then f=o E, =0 0.6 0.3 stress (ksi) 0.2 0.1 0.1 0.0001 strain in/in strain in/in Fig. 314. Typical loading in tension Eq. 3.47 Eq. 3.48 Eq. 3.49 Eq. 3. 50 Eq. 3.51 0 0005 0 000 0.00 Unloading Typical unloading phases in tension and compression are illustrated if Fig. 315 and Fig. 316, respectively. The expressions for compression and tension unloading are based on the work of Sinha et al. (1964a) defined for different concrete mixes. Due to the lack of more experimental data, it is assumed that for stronger concrete (fc > 4 ksi), the response will be that of 4 ksi concrete. If new experiments are done for stronger concrete these changes can be easily incorporated into the program. The family of unloading curves is represented by second order curves. From experiments for various concrete strengths a good fit was obtained with the formula: J a+H= ( X)2 Eq. 3.52 X where H and J are experimental constants whose values for the three mixes used are given in Table 31. For stress in units of ksi and strains in units of in/in, and X is a particular parameter, different values of which represent different members of the family. To determine the value of the parameter X, for a curve passing though a certain point in the stressstrain plane oaE,, the coordinate values are substituted into Eq. 3.52 which is then solved for the required value of X, leading to the expression: o,+H + H, 2 2) X=E, + H T + H _V Eq. 3.53 2J 2J The tangent modulus Ec for unloading is taken as the slope for two consecutive points in the unloading curve: E, Eq. 3. 54 6,, E Table 31. Curves coefficients fc (psi) H 3000 0.07 3750 0.09 4000 0.10 J 0.95 0.52 0.61 K 3.42 2.52 4.61 L 1.26 1.03 1.01 There are two choices for when the unloading in compression crosses the strain axis. The first option assumes that a gap is formed in compression, and concrete will totally unload until reloading in tension, as illustrated in Fig. 318. The second option assumes that no gap is formed, and concrete will go straight into the tension reloading phase from compression, as shown in Fig. 317. When unloading from tension a gap is always formed, and concrete will not go into compression until the gap is closed, as shown in Fig. 315. 0.6 0.5 oas : i, .. , 04_. stress 0, o 2 i.. : .. o : .... I i , 0.1 0 Q.01 0.0001 0 0.0001 Fig. 315. Typical unloading in tension 0.0003 0.0005 0.0007 0.0009 strain (inrin) 0.5 1 stress (ksi) 1.5 2 2.5 i 0.0006 0.0005 0.0004 0.0003 0.0002 0.0001 strain (in/in) Fig. 316. Typical unloading in compression 1 0.5 ... .. _ (ksi) 0.5....... 1.5 2 ..... 0.0006 0.0005 0004 0.0003 0.0002 0.0001 strain (in/in) Fig. 317. Compression unloading with gap 0.2 0.4 : i , 0.6 ..... 0.8 stress 1 (ksi)  1.2 1.6 2 0.0006 0.0005 0.0004 0.0003 0.0002 0.0001 0 strain (in/in) Fig. 318. Compression unloading with no gap Reloading A typical reloading in compression is illustrated in Fig. 319. The reloading curves in either compression or tension are represented by a family of converging straight lines; accordingly, an expression: a + K = Y(e + L) Eq. 3. 55 was chosen, in which K and L are experimental constants, and Y, is a parameter. The value of Y can be found, similarly to the value of X in the previous section. By substituting the coordinates o,,E, of a known point on the reloading curve into Eq. 3.55; and solving for Y, which is also take as the tangent modulus, Ec, one obtains: o, +K Y = +Eq. 3. 56 + L 0 1.5 2 0.006 0.0005 0 004 0.0003 0.0002 0.0001 0 strain (in/in) Fig. 319. Typical loading, unloading and reloading in compression When concrete goes into tension, if it goes back to compression it will reload only when the gap resulting form the previous cycle is closed. The process can be visualized in Fig. 320. stress strain Fig. 320. Concrete behavior with gap Bilinear Model In another study, Agrawal et al. (1965) claimed that the response of under reinforced concrete beams is governed by the steel. It was therefore proposed to use a simplification of the concrete model to simplify the analysis of doubly reinforced concrete beams. Idealized elastic plastic curves were drawn with a yield stress equal to the nominal strengthfc value for concrete, and an elastic modulus equal to the average stiffness of the initial portion of the actual stressstrain curve as shown in Fig. 321. Also, the tensile strength of concrete is neglected in this model. In FLPIER the stress strain curve for concrete is based on 13 points as shown in Fig. 322. Referring to Fig. 322, the value of the average elastic modulus is taken as Ec ( Eq. 3. 57 and the value of the yielding strain Ey is taken as a(2) Y Eq. 3. 58 E, Fig. 321. Bilinear model for concrete 1"oiue Fig. 322. Stressstrain curve for concrete in FLPIER Strain Rate Effect The strength of reinforced concrete sections is very dependent on the strain rate (Shkolnik (1996)), and degree of confinement of the section (Saatcioglu and Razvi (1992), and Soroushian et al. (1986)). At a fast strain rate, both the modulus of elasticity and the strength of concrete increase. The increase in strength can be as much as 17 % for a strain rate of 0.01/sec, as reported by Park and Paulay (1975). It is very difficult to test structural members under such conditions. So before using the results from static tests, it is important to consider the effect of strain rate in dynamic analysis. Important findings about the strain rate effects on reinfoced concrete members are summarized next (Otani, (1980)): 1) High strain rates increased the initial yield resistance, but caused small differences in either stiffness or resistance in subsequent cycles at the same displacement amplitude. 2) Strain rate effect on the resistance diminished with increased deformation in a strainhardening range. 3) Non substantial changes were observed in ductility and overall energy absorption capacity. Otani (1980) also suggests that the strain rate during an oscillation is highest at low stress levels, and gradually decreases toward a peak strain. Cracking, crushing and yielding will contribute to a reduction in the system's stiffness, elongating the period of oscillation. Such damage is caused by the lower modes of vibration having long periods. Therefore, the strain rate effect is small in earthquake analysis, and has a small effect on the response. So, the static hysteretic behavior can be used in a nonlinear dynamic analsyis of reinforced concrete structures. Confinement Effect The confinement of reinforced concrete columns is a good way to increase its strength and ductility, as shown in Fig. 323. Confined column Moment Curvatute Fig. 323. Confined and unconfined concrete models response Confinement basically decreases the slope of the descending branch of the loading curve of concrete, making the confined concrete member more ductile and less brittle. The reader is referred to Saatcioglu and Razvi (1992), and Soroushian et al. (1986) for a good discussion on the subject. Soroushian et al. (1986) proposed a very simplified model, that incorporates both, strainrate and confinement effects in the compression envelop curve for concrete. The following constitutive model (Soroushian et al. (1986)) was implemented in FLPIER: K 2E E Y KK 0.002K,K, 0.002K, K, f =8 0.002KK3 Eq. 3.59 K,K2 f[1z(e 0.002K,K,)] 2 0.002KK3 where f= concrete compressive stress E = concrete compressive strain K, =1+ V 4A p, = volumetric ratio of the hoop reinforcement to concrete core = E'" = _4A V sh' f ~= 28day compressive strength of concrete, adopted as 0.85fc fih = yield strength of transverse reinforcement 0.5 3 + .002f 5 (psi) +3 3p 0.002KK z = f 1000 +4 0.002KK3 Eq. 3.60 (MPa) 3 + 0.29M_ 3 h' +145 1000 p, 0.002K,K, 145f, 1000 4 sf h'= width of concrete core measured to outside of the transverse reinforcement, as shown in Fig. 3.24 for square, rectangular, and circular sections. h h' ' Fig. 324. Core width for different cross sections s = centerto center spacing of transverse reinforcement K2 = 1.48 + 0.1601ogo0 + 0.0127(log0 g)2 Eq. 3.61 K, = 1.08 +0.1601og0io + 0.0193(logi0 d)2 Eq. 3. 62 Note: For 9 < 10s / sec, K2 = K3 = 1.0. In this model, K2 represents the strain rate effect on the compressive strength of the concrete, and K3 takes care of the strain rate effect on the strain at maximum stress. It is assumed that the strain rate effect on the slope of the descending branch of the stressstrain diagram is similar to the strain rate effect on the compressive strength of concrete. This is supported by test results. The model representation can be seen in Fig. 325 below. There is no change in the unloading or reloading curves. Concrete Compressive Stress (f) KIKf'c 0.2KKf z 0.002KiK, Concrete Compressive Strain (e) Fig. 325. Confined concrete model The following modifications were also proposed by Soroushian et al. (1986) for the secant and tangent stiffness of concrete when subjected to dynamic loading: E = 1.241+ 0.111ogio + 0.127(logl0 g)2 Eq. 3.63 Ecs E = 1.061+ 0.464 log,1 + 0.00683(logo s)2 Eq. 3.64 E,, Where Ecd = dynamic secant modulus of elasticity, E,, = static secant modulus of elasticity, Eld = dynamic tangent modulus of elasticity, and E,, = static tangent modulus of elasticity. Note that when compared to the secant modulus, the tangent modulus seems to be less influenced by the rate of straining. These changes were also implemented in FLPIER. CHAPTER 4 MODAL ANALYSIS Because the most frequently used procedure for designing bridges for earthquakes is based on modal analysis, a brief explanation of this method of analysis is given next. It should be noted that modal analysis is an extensive subject and that what is presented here is only an introduction of the basic concepts necessary to understand the method. The interested reader is referred to Chopra (1995), or Paz (1985), for a more detailed discussion on modal analysis. Natural Vibration Frequencies and Modes Before getting to modal analysis it is opportune to introduce the eigenvalue problem whose solution gives the natural frequencies and vibration modes of a system. The free vibration of an undamped system in one of its natural vibration modes can be described mathematically by d(t) = q, (t) Eq. 4. 1 where the deflected shape n, does not vary with time. The time variation of the displacements is described by the simple harmonic function q,(t) = A, cosco,t + B, sinco,t Eq. 4.2 where A, and B, are constants of integration that can be determined from the initial conditions that initiate the motion. Combining Eqs. 4.1 and 4.2 gives 
Full Text 
xml version 1.0 encoding UTF8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd INGEST IEID EV79Y7PH8_ARQZD9 INGEST_TIME 20130122T13:17:33Z PACKAGE AA00012993_00001 AGREEMENT_INFO ACCOUNT UF PROJECT UFDC FILES 