THE DYNAMIC BEHAVIOUR OF THE CENTRE OF STIFFNESS OF R / C ECCENTRIC STRUCTURES UNDER SEISMIC EXCITATION Raymond K.W. Tong B. A. Sc. (Mechanical Engineering) University of British Columbia, 1982 B. A. Sc. (Civil Engineering) University of British Columbia, 1985 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F A P P L I E D S C I E N C E in T H E F A C U L T Y O F G R A D U A T E S T U D I E S D E P A R T M E N T O F C I V I L E N G I N E E R I N G We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A April 1988 © Raymond K.W. Tong , 1988 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of rJul f £noj/^gfK/'it, The University of British Columbia Vancouver, Canada DE-6 (2/88) a b s t r a c t Under seismic excitation, the centre of stiffness of an eccentric structure is stationary as long as the structure remains elastic. Once yielding occurs, the centre of stiffness will begin to move away from its original position, as the torsional forces induced by the eccentricity cause uneven distribution of yielding among members. This movement of the centre of stiffness very often increases the eccentricity of the structure causing further damage. The purpose of this thesis is to determine the significance of this magnification of the initial eccentricity. A procedure for locating the centre of stiffness was developed and incorporated into a time-step dynamic analysis program code named DRAINTABS. Two models were chosen to describe the moment-rotation relationship of reinforced concrete members; the elasto-plastic model and the Takeda model. The former is a bilinear model whereas the latter takes into account the strength degradation of reinforced concrete members under cyclic loading. A number of reinforced concrete buildings were studied. It was found that the centre of stiffness did not always move towards the side of the structure which was more heavily loaded due to the torque induced by the eccenticity. Excursions in the other direction were possible when the translational motion was not in phase with the torsional motion at the instant of maximum excursion. Moreover, when the strength degrading characteristic of R / C members was modelled, the eccentricity of the structure increased gradually with increasing length of excitation. However, this increase was found to be small and therefore insignificant. The procedure for locating the centre of stiffness was also incorporated into another ii analysis program code named PITSA which utilizes the modified substitute structure method. The results obtained were compared to those obtained using DRAINTABS. Although PITSA has been proven to be a relatively inexpensive yet reliable alternative to a time step analysis, it failed to predict the maximum displacement of the centre of stiffness with any acceptable degree of accuracy. iii Table of Contents abstract ii table of contents iv list of tables viii list of figures ix dedication xii acknowledgement xiii 1 I N T R O D U C T I O N 1 2 B R I E F R E V I E W O F P I T S A 3 2.1 Modified Substitute Structure Method 3 2.1.1 Damage Ratios 3 2.1.2 Substitute Stiffness and Damping Ratio 4 2.1.3 Computational Procedure 5 2.1.4 Convergence Scheme 6 2.2 Structural Idealization 7 2.3 Building Stiffness Matrix 8 2.4 Response Spectrum 8 3 B R I E F R E V I E W O F D R A I N T A B S 9 iv 3.1 Structural Idealization 9 3.2 Condensation of Frame Stiffness and Load Vector 10 3.3 The Building Stiffness Matrix and Load Vector 12 3.4 Step-by-Step Dynamic Analysis 12 3.5 Earthquake Excitation 13 3.6 Calculation of Damage Ratios 13 4 C O M P A R I S I O N S O F R E S U L T S 15 4.1 Structure 1: 3-Storey Framed Structure 16 4.2 Structure 2: 3-Storey Structure with Identical Frames 17 4.3 Structure 3: 5-Storey Coupled Wall Structure 18 4.4 Structure 4: 15-Storey Framed Building 19 4.5 Conclusions 19 5 L O C A T I N G T H E C E N T R E O F S T I F F N E S S 21 5.1 General Procedure for Locating the Centre of Stiffness 21 5.1.1 Extracting the Building Stiffness Matrix 22 5.1.2 Condensation of the Stiffness Matrix 23 5.1.3 Locating the Centre of Resistance 24 5.2 Additional Output from DRAINTABS and PITSA 25 5.3 Program Testing 26 5.3.1 Test Structure 26 5.3.2 Verification of the Building Stiffness Matrix 26 5.3.3 Discussions of Results 27 6 C A S E S S T U D I E D 29 v 6.1 Structure 1: 3-storey Structure Designed in Accordance with the Na-tional Building Code 29 6.1.1 DRAINTABS (Takeda Model) 29 6.1.2 DRAINTABS (Elasto-Plastic Model) 31 6.1.3 Comparisons of Results Obtained Using DRAINTABS and PITSA 32 6.1.4 Effect of Eccentricity on the Behaviour of the Centre of Stiffness 33 6.2 Structure 5: 3-Storey Structure with Weaker Columns 34 6.2.1 DRAINTABS (Takeda model) 34 6.2.2 DRAINTABS (Elasto-Plastic model) 35 6.2.3 PITSA 36 6.3 Structure 2: Three-Storey Structure with Identical Frames 36 6.3.1 DRAINTABS (Takeda-Model) 36 6.3.2 DRAINTABS (Elasto-Plastic Model) 37 6.3.3 PITSA 37 6.3.4 Effect of Eccentricity 37 6.4 Structure 6: Three-Storey Structure with Asymmetrical Stiffness Distri-bution 38 6.5 Conclusions 40 7 C O N V E R G E N C E P R O B L E M S IN U S I N G P I T S A 42 7.1 Convergence Speeding Routine 42 7.2 The Solution Searching Routine 45 7.3 Conclusions 47 8 C O N C L U S I O N S 49 vi B i b l i o g r a p h y A L i s t i n g o f S u b r o u t i n e s A d d e d t o D R A I N T A B S v i i List o f Tables 6.1 Comparison of maximum displacements of CR predicted by DRAINTABS and PITS A for structure 1 98 6.2 Comparison of maximum displacements of CR predicted by DRAINTABS and PITSA for structure 5 98 viii L i s t o f F i g u r e s 2.1 Definition of Damage Ratio 51 2.2 Plan view of floor showing frame and diaphragm horizontal displacement (PITSA) (from ref. 3) 52 2.3 Spectrum A (from ref. 3) 53 3.1 Diaphragm and frame horizontal displacements(DRAINTABS) (from ref.3) 54 4.1 (a) Member properties of Structures 1,2 and 3 55 4.1 (b) Layout and dimensions of structures 1, 2 and 3 56 4.2 (a) Member strengths of structure 1 57 4.2 (b) Member strengths of structure 1 58 4.3 (a) Member strengths of structure 2 59 4.3 (b) Member strengths of structure 2 60 4.4 Dimensions and properties of structure 3 61 4.5 (a) Member properties of structure 4 62 4.5 (b) Layout and dimensions of structure 4 63 4.6 El Centro NS earthquake record, 1940 64 4.7 (a) Damage ratios for structure 1 65 4.7 (b) Damage ratios for structure 1 66 4.7 (c) Damage ratios for structure 1 67 4.8 (a) Damage ratios for structure 2 68 4.8 (b) Damage ratios for structure 2 69 4.8 (c) Damage ratios for structure 2 70 ix 4.9 Damage ratios for structure 3 71 4.10 (a) Damage ratios for structure 4 72 4.10 (b) Damage ratios for structure 4 73 4.10 (c) Damage ratios for structure 4 74 5.1 Locating the centre of stiffness from the condensed stiffness matrix . . . 75 5.2 Layout and member data of test structure 76 5.3 Time history plot of CR for test structure 77 6.1 Hinge moment-rotation relationship for Takeda model (from ref. 7) . . . 78 6.2 Elasto-plastic moment-curvature relationship 79 6.3 (a) Time history plot of CR for structure 1 80 6.3 (b) Time history plot of CR for structure 1 81 6.4 Graph of maximum displacement of CR vs eccentricity for structure 1. . 82 6.5 (a) Time history plot of CR for structure 5 83 6.5 (b) Time history plot of CR for structure 5 84 6.6 (a) Time history plot of CR for structure 2 85 6.6 (b) Time history plot of CR for structure 2 86 6.7 Graph of maximum displacement of CR vs eccentricity for structure 2. . 87 6.8 Time history plot of CR for structure 6 88 6.9 Plots of displacement, rotation and CR of the first diaphragm of struc-ture 6 89 7.1 Flow-chart for convergence speeding routine (from ref. 2) 90 7.2 (a) Plots of damage ratio error and natural period for structure 5. . . . 91 7.3 (b) Plots of damage ratio error and natural period for structure 5. . . . 92 7.4 (c) Plots of damage ratio error and natural period for structure 5. . . . 93 7.5 (a) Damage ratios for structure 5 94 7.5 (b) Damage ratios for structure 5 95 x 7.5 (c) Damage ratios for structure 5 96 7.6 Flow-chart for solution searching routine STACHK (from ref. 6) . . . . 97 xi TO M Y MOTHER xii a c k n o w l e d g e m e n t The author would like to express his sincere gratitiude to his advisor, Dr. N.D. Nathan, for his valuable advice and guidance throughout this research. Thanks are also extended to Dr. R. Spencer for his conscientious proofreading. The help of my fellow graduate students during the preparation of this thesis is greatly appreciated. This research was funded by the National Research Council of Canada in the form of a Research Assistantship. xiii C h a p t e r 1 I N T R O D U C T I O N If a lateral force is applied at a point in a horizontal section of a structure, and the section moves in the same direction as the force, then, by definition, this point is the centre of stiffness (or the centre of resistance) of the section. The significance of the centre of stiffness is that any force that is not applied through it will create torsional motion of the section. It is usually assumed that the seismic design forces for each floor of a structure act through the centre of mass of the floor, and that the floor performs torsional oscillations about its centre of mass as well. For eccentric structures whose centres of mass do not coincide with their centres of resistance, additional torsional forces will be created. These forces are directly proportional to the amount of eccentricity, and are known as long as the structure remains elastic. In design, it is necessary to account for such torsional responses which may induce additional shear forces on the lateral resisting elements of the building. Most seismic codes for buildings recognize this necessity and have made provisions for torsional effects. The most common form of torsional provision is the requirement that an additional loading condition due to torsional moments or torques at each storey be considered. The torsional moment at each storey is obtained by multiplying the storey shear by a quantity termed 'design eccentricity'. In most seismic codes, the design eccentricity expression consists of two parts. The first part is termed dynamic eccentricity and is a function of the distance between the centre of mass and centre of stiffness. This dynamic eccentricity takes into account the difference in 1 Chapter 1. INTRODUCTION 2 stiffness and mass distributions of a building. The second part, commonly referred to as accidental eccentricity, accounts for other factors such as torsional ground motion, and the differences between actual and computed eccentricities. The accidental eccentricity is a function of the plan dimension of the building. Once yielding begins, structural elements that have yielded become softer, causing the centre of stiffness to move away from them. This shifting of the centre of stiffness may increase the eccentricity of the building, and thus create an extra amount of torque which in turn may cause further yielding of the structural elements that yielded initially with the whole cycle being repeated. It can be seen that the torsional motion of the structure can be substantially magnified by the movement of the center of resistance. It is the main objective of this thesis to study the dynamic behaviour of the centre of stiffness for a number of eccentric reinforced concrete structures in the non-linear range, and therefore to assess the extent to which the initial eccentricity is magnified by the movement of the centre of stiffness. It should be noted that buildings which are not eccentric during elastic response may become eccentric during nonlinear response if the structural elements on one side yield before the other . Chapter 2 BRIEF REVIEW OF PITSA The code name P I T S A stands for Pseudo Inelastic Torsional Seismic Analysis. It is an analysis program for the evaluation of inelastic responses and ductility demands of three dimensional buildings under seismic excitation. It was developed as a less expensive but reliable alternative to a time-step dynamic analysis. The analysis uses the modified substitute structure method which was developed by Yoshida 1 for the analysis of existing reinforced concrete buildings. It was later refined and extended by Metten 2 to the analysis of seismic resistant coupled structure walls. Tarn 3 further extended the method to three-dimensional structures. 2.1 Modified Substitute Structure Method To evaluate the inelastic responses of a structure, P I T S A replaces the real structure with a fictitious elastic structure with its stiffness and damping properties related to but differing from those of the real structure. Then the member forces are obtained from a modal spectral analysis of the substitute structure using a linear response spectrum. 2.1.1 Damage Ratios The damage ratio of a member is the ratio of its initial stiffness to the secant stiffness in its final configuration as shown in fig.2.1. A bilinear model with strain hardening equal to s is used to describe the moment-rotation relationship for the member. It should be noted that damage ratio is not equivalent to ductility ratio which is the ratio of yield 3 Chapter 2. BRIEF REVIEW OF PITSA 4 to maximum rotation. They are numerically the same only for the elastoplastic case where strain hardening is zero. The relationship between damage ratio and ductility ratio is given by where p, = damage ratio n = ductility ratio s = strain hardening ratio 2.1.2 S u b s t i t u t e St i f f n e s s a n d D a m p i n g R a t i o The flexural stiffness of substitute frame members is defined as, (EIU = ( 2 . 2 ) where [EI)si = cross-sectional flexural stiffness of the ith member in the substitute frame (EI)ai = cross-sectional flexural stiffness of the ith member in the actual frame fi — selected tolerable damage ratio for the ith member The term [EI)a is calculated using the fully cracked section. The substitute damping factor for the ith element is approximated by the following equation: Pi = 0.02 + 0.2(1 - -^=) (2.3) V M i This value of viscous damping represents the hysteric energy dissipated by the member. To obtain a composite damping ratio for the entire structure, a 'smeared damping Chapter 2. BRIEF REVIEW OF PITSA 5 ratio' is computed for each mode by assuming that each element contributes to modal damping in proportion to its relative flexural strain energy associated with the mode shapes. The flexural strain energy of the ith element in the mth mode is given by the formula, n,- = ^ T J ^ : (Ml + Ml - MaiMH) (2.4) where Mai, Mbi ' end moments of the ith member. The smeared damping ratio for the mth mode is then given by, Pm- E < n r (2.5) It can be seen from the above equation that members with large bending moments and/or large damping ratios are major contributors to the smeared damping ratio. 2.1.3 C o m p u t a t i o n a l P r o c e d u r e The modified substitute method involves an iterative procedure which can be summa-rized as follows: 1. A elastic modal analysis is performed in the first iteration with the damage ratio set to unity for all members. The damping ratio is arbitrarily assigned a value as it will not affect the final results. 2. Member moments for each mode are combined using the root-sum-square and complete-quadratic-combination methods. The damage ratios of those members whose CQC moments exceed their yield values are modified using the following formula: Mnpn . . = TTR— irr ^ 1 (2-6) My[l - s) + sfinMn Chapter 2. BRIEF REVIEW OF PITSA 6 where fin+i = damage ratio in the (n + l)th iteration Hn = damage ratio in the nth iteration Mn = CQC moment in the nth iteration My = yield moment 5 = strain hardening ratio 3. Starting with the second iteration, a smeared damping ratio is calculated for each mode and a new set of damage ratios are obtained. With the new damage ratios, a new stiffness matrix is calculated and another iteration is performed to further refine the damage ratios. The iteration is halted when the convergence criteria described in the next section are satisfied. 2.1.4 C o n v e r g e n c e Scheme To determine if the iterative process has indeed converged to the final solution, the following criteria are used: 1. The bending moment error which is the difference between the computed bending moment and the moment' capacity of a yielded member must be less than 5% of capacity for all members i.e. M n ~ M c a p < 0.05 if M > 1 (2.7) jyicap 2. The change in damage ratios between successive iterations must be less than 1 % of the damage ratio in the current iteration for members with damage ratios greater than five. In the case of a member with damage ratio less Chapter 2. BRIEF REVIEW OF PITSA 7 than five, the absolute difference of the damage ratios between successive iterations is to be limited to 0.1. The criterion can be summarized as follows: Hn-Hn-l Mr. < 0.01 if A* > 5 or (2.8) |j*n - A*n-l| < 0.1 if fJ, < 5 It is worthwhile to note that equation 2.8 is always the governing criterion for damage ratios greater than two. 2.2 S t r u c t u r a l I d e a l i z a t i o n The structural idealization in PITSA can be summarized as follows: 1. The structure is idealized as a series of arbitrarily oriented plane frames con-nected by horizontal rigid diaphragms. Isolated shear walls are considered as frames consisting of one column line having the associated wall proper-ties. Coupled shear walls are considered as frames consisting of two column lines and one bay of beams with rigid end section reaching from the wall centreline to the face of the wall. 2. The mass of each floor is lumped in the diaphragm at the centre of mass for horizontal inertia effects. 3. All torsional stiffnesses of components are ignored. 4. The compatibility of vertical displacements for columns common to two or more frames is enforced. However, no joint rotational compatibility is enforced for common columns and therefore some error is introduced when the frames are not mutually perpendicular. Chapter 2. BRIEF REVIEW OF PITSA 8 5. The moment-rotation relationship is assumed to be bilinear with no stiffness degradation. 6. No axial-fiexural interaction is assumed for columns. 2.3 Building Stiffness Matrix The procedure for assembling the building stiffness matrix is very similar to that in DRAINTABS which will be briefly described in chapter 3. Individual member stiffness matrices are combined to form frame matrices which are then condensed to eliminate those degrees of freedom which are not associated with inertia forces. The condensed stiffness matrices are then transformed from the local to the global coordinate sys-tem which consists of two translations and one rotation at the centre of mass of each diaphragm as shown in fig.2.2. These transformed matrices can then be combined together to form the building stiffness matrix. 2.4 Response Spectrum Six choices of response spectrum are provided in PITSA. Spectrum A which is shown in fig.2.3 will be used in this thesis. It is acquired by deriving averaged or smoothed spectra from four real earthquake records scaled to the same maximum acceleration. The four records are El Centro EW, El Centro NS, Taft N69W,and Taft S21W. C h a p t e r 3 B R I E F R E V I E W O F D R A I N T A B S DRAINTABS is a computer program for the inelastic dynamic response analysis of three-dimensional buildings under seismic excitation. It was developed at the Univer-sity of California Berkeley in 1976 by Rafael Guendelman and Graham H. Powell. This chapter describes briefly the computational procedure as well as the structure ideal-ization assumed in the program. A more detailed description of the program can be found in reference 4. 3.1 S t r u c t u r a l I d e a l i z a t i o n The structural idealization assumed in DRAINTABS can be summarized as follows: 1. Each structure consists of a series of independent substructures intercon-nected by rigid diaphragms. Each substructure can be of any arbitrary geometry and include structural elements of a variety of types. It is not necessary for all substructures to connect to all diaphragms. 2. The program does not enforce compatibility of vertical and rotational di-aplacements at joints common to two or more frames, and therefore the coupling of substructures through common columns cannot be fully ac-counted for. This type of idealization may not be accurate for structures such as framed tubes in which there is substantial coupling through common columns. To correct for common-column effects, the axial forces calculated 9 Chapter 3. BRIEF REVIEW OF DRAINTABS 10 for a column common to two or more frames are added together and the combined force is used to assess the P - M interaction effects. This is only an approximate correction procedure, and is valid only for structures whose behaviour is not significantly affected by axial deformations of the columns. 3. The mass of the structure may be lumped at the joints of the frames for vertical inertia effects, but should preferrably be lumped entirely in the diaphragms for horizontal inertia effects. 4. The displacement degrees of freedom for any frame are divided into two groups: internal degrees of freedom and connected degrees of freedom. The horizontal displacements at those joints which are connected to the di-aphragms are the connected degrees of freedom and all other displacements are internal degrees of freedom. 3.2 C o n d e n s a t i o n o f F r a m e Stiff n e s s a n d L o a d V e c t o r Frame stiffness matrices are assembled by combining individual element stiffness matri-ces. Six types of elements are available in DRAINTABS: truss element, beam-column element, infill panel element, semi-rigid connection element, beam element and beam element with degrading stiffness. The equilibrium equations for a single frame can be written in matrix form as: (3.1) where rj = internal degrees of freedom Chapter 3. BRIEF REVIEW OF DRAINTABS 11 rE = connected degrees of freedom Ri = load vector corresponding to rj RE = load vector corresponding to TE KH,KIE,KEI,KEE = frame stiffness submatrices The frame load vector is assembled by adding the loads applied directly at the nodes and the loads originating within elements as fixed-end forces. For each frame, a stiffness and a load vector in terms of the connected degrees of freedom only are obtained by condensing out the internal degrees of freedom, 77. The condensation process is as follows: Multipling out the submatrices in eqn.3.1 gives \Kn\in} + \KIE\{rE} = {Rj} (3.2) [KEi]{ri} + [KEE]{rB} = {RE} (3.3) Solving eqn.3.2 for 77 in} = [KT?]{Ri} - [KJ^][KIE]{rE} (3.4) Substituting into eqn.3.3 [KEI] {{KTr'KRj} - [Kj})[KIE\{rE}} + {KEE){rE} = {RE} (3.5) Hence [-Kcond]{»*E} = {Rcond} (3.6) where [Kcond] = \KEE)-\KEI)[K£]\KIE) (3.7) {Rcond} = {REy-lKEj^Kj^Rj} (3.8) Chapter 3. BRIEF REVIEW OF DRAINTABS 12 3.3 The Building Stiffness Matrix and Load Vector To form the stiffness matrix for the whole building, all the condensed frame matrices are transformed to a common displacement system in terms of the degrees of freedom of the diaphragms i.e. two displacements and one rotation at the centre of mass of each diaphragm as shown in fig.3.1. The centres of mass of the diaphragms need not lie on a vertical line. The building stiffness matrix can then be assembled by combining all the transformed frame stiffness matrices together. The building load vector consists of two components: loads applied to the individual frames and loads applied directly to the diaphragms. Loads applied to the frames are condensed during the static condensation process, and must then be transformed to the diaphragm system. 3.4 Step-by-Step Dynamic Analysis At any instant of time, the basic equation of dynamic equilibrium can be written in incremental form as { A x } , {Ai ;} , {Arc} = increments of acceleration,velocity and displacement, respectively [M]{Ax} + [CT]{Ax} + [KT]{Ax} = -\M\{Axa} (3.9) where { A * , } increment of ground acceleration [M] building mass matrix [CT] tangent damping matrix [KT] tangent stiffness matrix Chapter 3. BRIEF REVIEW OF DRAINTABS 13 The dynamic response is determined by a step-by-step integration of the equations of motion. The constant acceleration model which assumes a constant absolute ac-celeration within a time step is used. The type of viscous damping specified may be mass and/or stiffness dependent. Once the diaphragm displacements are known, they are transformed back to the individual frames and the member forces can then be calculated. 3.5 Earthquake Excitation Two independent horizontal ground motions plus a vertical ground motion may be specified as the earthquake excitation. In order to obtain the required ground mo-tions, input earthquake records can be scaled with respect to both time and ground acceleration. 3.6 Calculation of Damage Ratios To assess the extent of damage of a yielded structure using DRAINTABS, the damage ratios of all members have to be calculated using the output from the program. The damage ratio of a member is defined as the ratio of its initial stiffness to the secant stiffness in its final configuration as shown in fig.2.1. By dividing the moment cor-responding to maximum rotation by the maximum rotation, the secant stiffness and therefore the damage ratio can be calculated. For the elastoplastic case where there is no strain hardening, the damage ratio is numerically equal to the ductility ratio. Hence, with the assumption that the member undergoes reverse bending and that the point of inflection is near midspan, the damage ratio can be directly calculated from the output member rotations as follows: p, = Oy + Op By if s = 0 (3.10) Chapter 3. BRIEF REVIEW OF DRAINTABS 14 where By = rotation at yield ML 6EI 0p = plastic hinge rotation (3.11) DRAINTABS outputs the plastic hinge rotations of members and. the rotations at yield have to be calculated, using eqn.3.11. C h a p t e r 4 C O M P A R I S I O N S O F R E S U L T S Although the modified substitute structure method has been proven by Tarn3 to be a reliable alternative to time-step dymamic analysis in predicting the inelastic behaviour of multi-storey three dimensional structures, the accuracy of the method depends to a certain extent on a number of factors, such as the structure type and the response spectrum used. Moreover, results from DRAINTABS have also been shown to be dependent on the earthquake record used. Since both DRAINTABS and PITSA will be used in the next two chapters to study the inelastic behaviour of the centre of resistance, it is of interest to examine the agreement between the two programs on the ductility requirements of the structures that are to be used later on in this thesis. To describe the amount of inelastic deformation in flexural members, two indices can be used; the curvature ductility and the member ductility. The curvature ductility is the ratio of the curvature at ultimate to the curvature at yield moment. It depends on the stress-strain characteristics of the concrete and steel in the member. In using this approach, a plastic hinge length has to be assumed. It is usually taken as a fixed percentage of the depth of the member. A value of one-half of thie effective depth is said to be reasonable5. The member ductility is defined as the sum of the yield rotation and the maximum hinge rotation divided by the yield rotation. Since DRAINTABS outputs the hinge rotations at each end of every member, the member ductility based on end rotations was adopted in this study. Strain hardening was taken to be zero in both programs 15 Chapter 4. COMPARISIONS OF RESULTS 16 because only then is the member ductility numerically equal to the damage ratios output by PITSA. It should be noted that the equation used to calculate the end yield rotation of a member assumes that the member undergoes anti-symmetric bending at the ends as shown in fig.2.1. An examination of output from DRAINTABS showed that this was not always the case. The two ends of a member did not always yield simultaneously. Nevertheless, the moments at the two ends were close enough not to cause any significant error in the calculation. The design spectrum A was chosen for PITSA and the EL Centro NS earthquake record was used for DRAINTABS. Viscous damping was included at 2 % of the critical to represent the effect of non-structural elements, as this value is believed to be a reasonable value for damping in a concrete structure in the elastic range. Hysteretic damping of the structure itself is included automatically in the non-linear analysis. 4.1 S t r u c t u r e 1: 3-Storey F r a m e d S t r u c t u r e This was one of the structures used for program testing by Canisius5 in his master's thesis. The dimensions and member properties are shown in fig.4.1. It is a 2 bays by 3 bays three storey concrete framed building with an eccentricity of 54 inches in the x-direction. The member strength capacities were determined by sizing them ap-proximately according to the NBC code under a maximum acceleration of 0.3g. The strengths of the members are shown in fig.4.2. The member strengths under positive and negative bending were assumed to be the same. The El Centro NS earthquake record scaled to a maximum acceleration of 0.3g was used in DRAINTABS to excite the building in the y-direction. A plot of the acceleration record is shown in fig.4.6. The damage ratios obtained from the two programs are shown in fig.4.7. In most cases, the damage ratios predicted by PITSA were higher than those Chapter 4. COMPARISIONS OF RESULTS 17 predicted by DRAINTABS. The average discrepancy was rougly 20%. Nevertheless, PITSA was able to predict successfully the pattern of ductility requirements. The beams were found to have yielded much more extensively than the columns. How-ever, the agreement between the damage ratios obtained from the two programs was the same for the beams as for the columns. This indicates that where is no apparent correlation between the accuracy of PITSA and the ductility requirement of a member. 4.2 Structure 2: 3-Storey Structure with Identical Frames The layout of this structure is the same as structure 1. However, all the frames parallel to the y-axis(3,4,5and 6) were made identical such that they all had the same strengths as frame 3 of structure 1. This structure was chosen for two reasons. The first one is to determine if PITSA would give more accurate predictions for more uniform structures. Secondly, since all the frames were identical, the damage ratios calculated should rep-resent the loading pattern in the structure. The same earthquake record was used to excite the building in the y-direction. The results are shown in fig.4.8. Again, PITSA was able to predict correctly the damage pattern in the structure. As for structure 1, the damage ratios predicted by PITSA were on the average 20% higher than those predicted by DRAINTABS; same as for structure 1. The beams were found to have damaged much more severely than the columns with the amount of damage increasing with higher floors. Surprisingly, the maximum column yielding occured in column 4 of frame 6. It was expected that the corner column in frame 6 would be most severely damaged. Furthermore, since the building had an eccentricity in the x-direction, the frames on the same side as the centre of mass should yield more than the frames on the other side due to the torque induced by the eccentricity. This was confirmed by both programs; the damage ratios increased in the positive x-direction. Chapter 4. COMPARISIONS OF RESULTS 18 4.3 Structure 3: 5-Storey Coupled Wall Structure This structure was used by Tarn1 in his thesis for program testing. The dimensions and properties are shown in fig.4.4. The building had an eccentricity of 2 feet in the y-direction which represented a uniform eccentricity of 3 % throughout the height. The size of the core which consisted of two channel sections joined by two coupling beams, measured 30 feet by 18 feet. The diaphragms were 100 ft by 60 ft. Two earthquake records were used to excite the structure in the x-direction; the El Centro NS and the El Centro EW both scaled to a maximum acceleration of 0.3g. The results are shown in fig.4.9. As before, the predictions by PITSA were conservative. The agreement between the two programs was good for the columns. However, for the beams, the damage ratios predicted by PITSA were much higher. Nevertheless, PITSA was able to pick out the damage pattern very well. The damage ratios calculated using the El Centro EW record were higher than those obtained using the El Centro NS record. The reason is probably that for the El Centro EW earthquake, the acceleration peaked at a much later instant than that for the El Centro NS earthquake. Since the Takeda model was used in DRAINTABS to model the strength degradation of reinforced concrete members under cyclic loading, by the time the acceleration reached its peak, the strength of the building had already been weakened considerably. Consequently, extensive yielding occured in the members. On the other hand, the effect of strength degradation was not significant when the El Centro NS record was used as the peak acceleration occured very early at 2.12 sec. In fact, the design spectrum A used in PITSA was derived from six earthquake records. Therefore, better agreement would probably result if all six records were used and the average values calculated for comparisons. Chapter 4. COMPARISIONS OF RESULTS 19 4.4 Structure 4: 15-Storey Framed Building This structure will not be used later on in this thesis. It is studied here in order to determine the reliabilty of PITSA on tall and slender structures. The layout and member properties are shown in fig.4.5. It is a 15-storey framed structure with an eccentricity of 4 feet in the y-direction. The El Centro NS record was used to excite the building in the x-direction for a duration of 8 seconds. The results are shown in fig.4.10. The damage pattern was not predicted as well as in the previous cases. However, the agreement between the two programs was generally good. The damage ratios predicted by PITSA were not always higher than those predicted by DRAINTABS, as found in the previous cases. Nevertheless, for both the beams and the columns, the maximum damage ratio was predicted by PITSA. Also, PITSA correctly predicted that more yielding would occur in frame 4 which was closer to the centre of mass than frame 1. 4.5 Conclusions Based on the results obtained for the above structures, the following can be concluded: 1. PITSA is a reliable tool in predicting the damage pattern of framed as well as coupled wall structures within the limits of accuracy of practical seismic design. 2. In general, the predictions by PITSA are on the conservative side. However, the agreement between DRAINTABS and PITSA seems to be better for framed stuctures than for coupled wall structures. 3. The accuracy of PITSA does not seem to be dependent on the ductility requirements of the structure being analysed. Chapter 4. COMPARISIONS OF RESULTS 20 4. When using D R A I N T A B S , different earthquake records will produce slightly different results, even if all records are scaled to the same maximum accel-eration. Chapter 5 L O C A T I N G T H E C E N T R E OF S T I F F N E S S 5.1 General Procedure for Locating the Centre of Stiffness In general, the location of the centre of resistance during excitation differs from floor to floor. Output from DRAINTABS and PITSA do not give the location of the centre of stiffness explicitly. Therefore, subroutines have been incorporated in both programs to calculate the position of the centre of stiffness for each floor, using information provided by the programs. The procedure is outlined below: 1. Extract the building stiffness matrix [K] from the program. 2. Condense [K] into a 3 x 3 matrix [Kc] such that [ K c ] 3 x 3 { » " } 3 x l = {R}sxi condensed matrix column vector containing the three degrees of freedom of the i th floor for which the center of stiffness is 21 where [Kc] = w = Chapter 5. LOCATING THE CENTRE OF STIFFNESS 22 to be located > = column vector containing the three components of forces of the i th floor. 3. Reduce the first two columns of the condensed matrix which are the systems of forces required to move the floor in the x and y directions, respectively, to two resultant forces. 4. Locate the intersection of these forces, which will give the centre of resistance of the floor being investigated. 5.1.1 E x t r a c t i n g t h e B u i l d i n g S t i f f n e s s M a t r i x D R A I N T A B S Since a time-step analysis is employed in DRAINTABS, the building stiffness matrix is calculated for every time step, making it possible to study the behaviour of the center of resistance during the entire seismic excitation. Because the structure stiffness matrix is always symmetric, only half of it needs to be stored. In fact, as a storage saving measure, DRAINTABS stores the upper half of the stiffness matrix columnwise in compacted form ignoring all terms up to the first non-zero element in every column. To keep track of the addresses of the non-zero elements, the addresses of all the diagonal elements are stored. By knowing the addresses of all diagonal elements, the number of zeros in each column can be calculated. For instance, if the second diagonal element has an address of two, this means that the first element Chapter 5. LOCATING THE CENTRE OF STIFFNESS 23 in the second column must be zero (otherwise the address of the diagonal element would be three), and is therefore ignored. After the building stiffness matrix in compacted form has been extracted, it is expanded into a two dimensional array by subroutine TRANS. P I T S A Since the modified substitute structure method is used in PITSA, only the stiffness ma-trix of the substitute structure can be obtained from the program. As in DRAINTABS, only the upper half of the stiffness matrix is stored and is expanded into a two dimen-sional matrix by subroutine TRANS. The displacement of the centre of resistance from its original position found for the substitute structure is assumed to be the maximum displacement during the entire excitation period. 5.1.2 C o n d e n s a t i o n o f t h e Stif f n e s s M a t r i x Condensation of the centre of resistance is carried out by subroutine CONDEN. Before condensation is performed, rows and columns are interchanged so that the 3 x 3 sub-matrix in the left upper corner contains all elements belonging to the floor for which the centre of resistance is to be located. For instance, if the centre of stiffness of the third floor is to be found, rows 7 to 9 and columns 7 to 9 are interchanged with rows 1 to 3 and columns 1 to 3, respectively. The matrix is then partitioned as shown below: Koo Koi Kio K n where Koo is a 3 x 3 submatrix. Multiplying out the submatrices gives the following equations: [K0o]{ro> + [ K o ^ r J = {R0} (5.1) Chapter 5. LOCATING THE CENTRE OF STIFFNESS 24 [K10]{r0} + [KxiKrx} = {Rx} (5.2) Solving eqn.5.2 for {ri} = [Kn1] ({Ri} - [K10]{r0}) (5.3) Substituting into eqn.5.1, [Koo]{r0} + [Kcipr/] ({Ri} - [K10]{r0}) = {R0} or ([K00] - [KoiUK^JlKw]) {r0} = {R0} - M R " 1 ] ^ } or [Kc]3x3{>o}3xl = {R}sxi (5.4) where [Kc] = [Koo] - [KoiUK^pu] {R} = {R0} - [KoiUK^KRi} 5.1.3 Locat ing the Centre of Resistance By definition, any force that is applied through the centre of resistance in any direction will cause the structure to move in that direction without twisting. To find out the system of forces required to produce a unit displacement in the x-direction, {ro} can be set to {l 0 0} r in equation 5.4, and the resulting {R} should give the system of forces required. It should be obvious that the first two columns of the condensed stiffness Chapter 5. LOCATING THE CENTRE OF STIFFNESS 25 matrix contain the systems of forces necessary to produce a unit displacement in the x and y directions, respectively. Each of these systems consists of two forces in the x and y directions, and one moment; all of which are assumed to act through the centre of mass of the floor. In order to locate the centre of resistance, the two systems of forces contained in the condensed stiffness matrix must be reduced to two resultant forces. To do this, a resultant force is obtained for each system by finding the vector sum of the two translational forces. Then, to eliminate the moments, the resultant of each system is moved in a direction perpendicular to the force for a distance which is equal to the moment divided by the magnitude of the resultant force. This is illustrated in fig.5.1. Therefore, if a force is applied through point A in the direction shown, the floor will move in the x-direction. Exactly the same procedure can be applied to the second column of the condensed matrix, and the intersection of the two resultant forces obtained should be the centre of resistance of the floor, as shown in fig.5.1. The above procedure for locating the centre of resistance is performed by subroutine CR, and a listing of all the subroutines added to DRAINTABS and PITSA is given in Appendix A. 5.2 A d d i t i o n a l O u t p u t f r o m D R A I N T A B S a n d P I T S A In addition to the normal output from DRAINTABS, the locations of the centres of resistance of the floors specified by the user are given for every time step. Moreover, maximum displacement of the centre of stiffness from its original position, and the corresponding time step are also included. The initial position of the centre of stiffness is taken to be that corresponding to the first time step. In the case of PITSA, only the location of the centre of resistance of the substitute structure is given. As mentioned Chapter 5. LOCATING THE CENTRE OF STIFFNESS 26 earlier, the displacement of the the centre of stiffness of this structure should correspond to that of the real structure at the time when the centre of resistance is at its maximum excursion from its initial position. 5.3 Program Testing 5.3.1 Test Structure The layout and member data of the test structure are shown in fig.5.2. It is a two-storey framed concrete building with an eccentricity of two feet in the y-direction. The EL Centro NS earthquake record scaled to a maximum acceleration of .3g was used to excite the building in the x-direction for a duration of 5 seconds. 5.3.2 Verification of the Building Stiffness Matrix Both DRAINTABS and PITSA are very complicated programs and finding the building stiffness matrix is not as easy as one might expect. Therefore, the first logical step to verify the results is to check if the building stiffness matrix extracted is indeed the correct one. To calculate the stiffness matrix of a three-dimensional structure by hand could be an extremely tedious task. However, since the test building consists of identical beams and columns on all four sides, and is very simple in structure, calculating the stiffness matrix becomes manageable. In fact, this is the reason why such an unrealistically simple structure was chosen for testing. The hand-calculated stiffness matrix was compared with that extracted from both DRAINTABS and PITSA before yielding, and they were found to be almost identical. This verifies that the matrix taken out from both programs was indeed the building stiffness matrix. Chapter 5. LOCATING THE CENTRE OF STIFFNESS 27 5.3.3 Discussions of Results The behaviour of the centre of stiffness with time is plotted in fig.5.3, using the results obtained from DRAINTABS. The following observations can be made: 1. The program correctly predicted the centre of stiffness to be at the centre of the building before yielding began. 2. Since the building is symmetric with respect to the y-axis, the seismic forces exerted on both frame 2 and frame 3 are the same, and because these two frames are identical in strength, they will yield simultaneously. As a result no movement of the centre of stiffness is expected in the x-direction. This is verified by the results from DRAINTABS and PITSA. 3. Because the centre of mass of the building is on the plus y-side of the initial position of the centre of resistance, frame four will always receive more seismic force than frame one, regardless of the direction of the earthquake motion. Consequently, there will always be more yielding in frame four, causing the centre of stiffness to move towards frame one. This is verified by the results from DRAINTABS; all the excursions of the centre of resistance occured in the negative y-direction. 4. In DRAINTABS, the user can choose between two methods of modelling the stress-strain behaviour of reinforced concrete members; the Takeda model and the elasto-plastic model. Details of theses two models will be given in the next chapter. To analyse the test structure, the Takeda model was chosen to take into account the strength degradation of reinforced concrete members under cyclic loading. An inspection of the output file shows that with increasing length of excitation, the centre of stiffness did not oscillate Chapter 5. LOCATING THE CENTRE OF STIFFNESS 28 back to its initial position. Instead, the troughs of the oscillations gradually shifted towards frame one, indicating that the strength of frame four was indeed degrading under load reversals. However, this movement was too small to be detected on the plots. 5. Maximum displacement of the centre of resistance occured at 2.66 sec. for both floors, four-tenths of a second after peak excitation. 6. The behaviour of the centre of stiffness for both floors was quite similar, except that the second floor had slightly longer excursions. 7. An examination of the condensed matrices shows that there was at least one zero element in every column. To find out if the added subroutines would work for a full matrix, the structure was excited at 45° to the x-axis, and the centre of mass was relocated at a point which would give the building the same amount of eccentricity in both directions. By doing so, the condensed matrices will be made to contain no zero elements, and, furthermore, the centre of stiffness should move along a line at 45° to either axis. This was correctly predicted by both DRAINTABS and PITSA, proving that the subroutines will indeed work for this particular case as well. Based on the above observations, it is reasonable to conclude that the subroutines added to both programs can indeed correctly locate the centre of stiffness of a building. Chapter 6 C A S E S S T U D I E D 6.1 Structure 1: 3-storey Structure Designed in Accordance with the Na-tional Building Code 6.1.1 D R A I N T A B S ( Takeda Model) Under load reversals well into the inelastic range, the stiffness of a reinforced concrete member decreases due to concrete cracking and bond deterioration at the steel concrete interface. When the Takeda method is chosen in DRAINTABS, Takeda-type rules are used to describe the resulting hysteretic moment curvature relationships. Strain hard-ening and degrading flexural stiffness are approximated by assuming that the element consists of a linear elastic beam element with non-linear rotational springs at each end as shown in fig.6.1. Yielding may take place only in concentrated plastic hinges at element ends, and all plastic deformation effects, including the effects of degrading stiffness, are introduced by means of the moment-rotation relationships for the hinge springs. The moment-rotation relationship for each hinge is also illustrated in fig.6.1. The basic relationship is in the form of a bilinear curve, with an initial stiffness which is characteristic of monotonic loading condition. The degrading stiffness of the hinges is introduced when reversed loading is applied. The numbers on the legs of the diagram are yield codes which are part of the program output. The layout and member data of the structure are shown in fig.4.1. The structure is a 29 Chapter 6. CASES STUDIED 30 three-storey framed concrete building with an eccentricity of 54 ins. in the x-direction. The El Centro earthquake record scaled to a maximum acceleration of .3g was used to excite the building in the y-direction for a duration of 5 seconds. Because the centre of mass is to the right of the centre of resistance, under seismic excitation, the torque induced by the eccentricity will exert greater force on the frames on the right hand side of the centre of mass. In anticipation of the greater ductility demands on these frames, the building is designed such that the strengths of the frames (3,4,5 and 6) parallel to the y-axis increase in the positive x-direction. The behaviour of the centre of resistance with time for the three floors is plotted in fig.6.3(a). The following observations can be made: 1. Although the frames to the right of the centre of mass were made stronger, they still yielded more than those on the left initially. This is indicated by the fact that the first few excursions of the centre of resistance were to the left of its initial position. 2. After the first few excursions to the left, the centre of resistance started to oscillate about its initial position. Excursions to the right are possible because frames on the right hand side of the centre of mass have greater ductility demands only if the translational and torsional motion are in phase at the instant of maximum excursion. However, there are occasions when the translational and torsional motions become out of phase at the instant of maximum excursion and therefore make frames on the left subject to higher ductility demand. Currently, there is no clear cut criterion to predict whether such a condition will occur. Further study on this issue will be useful since very few seismic code provisions are directed to the design of elements on the same side as the centre of stiffness (measured relative to the Chapter 6. CASES STUDIED 31 centre of mass). 3. Maximum displacement of the centre of stiffness occured simultaneously at 2.4 seconds for all floors; two-tenths of a second after peak excitation. 4. The behaviour of the three floors was quite similar and the length of excur-sions increased with higher floors. 6.1.2 D R A I N T A B S ( E l a s t o - P l a s t i c M o d e l ) The moment-curvature diagram for this model is shown in fig.6.2. It is a much simpler model in that it assumes no strength degradation for reinforced concrete members under cyclic loading in the inelastic range. The time-history plots of the centre of resistance are shown in fig.6.3(b), and the following obsevations can be made: 1. In general, the behaviour of the centre of resistance was less oscillatory than that predicted by the Takeda model. This is probably due to the simplicity of the elasto-plastic moment-rotation relationship. The number of excursions is much smaller and they occurred predominantly on the left hand side of the centre of stiffness, suggesting that the frames to the right of the stiffness centre yielded more extensively than those on the other side, and that the translational and torsional motions were mostly in phase at instants of maximum excursion. However, the excursions predicted using the Elasto-Plastic model were generally longer than those obtained using the Takeda model. 2. For all three floors, the instant of maximum displacement occurred simulta-neously at 4.7 sec; much later(1.8 sec.) than that predicted by the Takeda model. The reason is probably due to the fact that the Takeda model takes Chapter 6. CASES STUDIED 32 into account the degradation in strength of reinforced concrete members un-der cyclic loading. As the frames on one side started to yield, the centre of stiffness began to move towards the other side. However, the strengths of the frames on the non-yielding side had already been weakened by previous load reversals. Since the position of the stiffness centre depends on the rel-ative stiffness between the frames on both sides of the centre of mass, the softening of the frames on the left reduced the maximum displacement of the stiffness centre to the left. This effect becomes more pronounced with time as the structure has undergone more load reversals. This explains why the lengths of the excursions predicted by the Takeda model were substantially shorter,and maximum displacement of the centre of stiffness occured at an earlier time. 3. All three floors exhibited very similar behaviour of the centre of stiffness. As found in the previous case, the length of excursions increased with higher floors. 6.1.3 C o m p a r i s o n s o f R e s u l t s O b t a i n e d U s i n g D R A I N T A B S a n d P I T S A Since PITSA employs the modified substitute structure method, it is not possible to study the behaviour of the stiffness centre with time. However, the displacement of the stiffness centre from its initial position for the substitute structure can be used to compare with the maximum displacement obtained using DRAINTABS. The results are tabulated in table 6.1. Because the Takeda model is a more realistic representation of reinforced concrete members, it was chosen in DRAINTABS for comparison purposes. As shown in table 6.1, the maximum displacements predicted by PITSA were much smaller than those predicted by DRAINTABS. Moreover, contrary to DRAINTABS, Chapter 6. CASES STUDIED 33 results from PITSA show that the length of maximum excursions decreased with higher floors. 6.1.4 E f f e c t o f E c c e n t r i c i t y o n t h e B e h a v i o u r o f t h e C e n t r e o f Sti f f n e s s To examine the effect of eccentricity on the behaviour of the centre of stiffness, the centre of mass of the building was placed at various distances from the initial position of the centre of stiffness and DRAINTABS was used to calculate the maximum dis-placement of the centre of stiffness corresponding to each eccentricity. The results are plotted in fig.6.4. At zero eccentricity i.e. when the centre of stiffness coincided with the centre of mass, maximum excursion occured to the right of the initial position of the centre of stiffness, indicating that there was more yielding in the frames on the left (frames 3 and 4). This was expected because the building was designed such that the frames on the right hand side were stronger in anticipation of the heavier seismic loading on them. However, as the centre of mass was moved to the right away from the centre of stiffness, the frames on that side became more heavily loaded as a result of the increasing eccentricity. Consequently, the centre of stiffness moved further and further to the left. As shown by the curve in fig.6.4, the variation of the maximum displacement with eccentricity is fairly linear; there are no drastic changes of slope along the curve. If the centre of stiffness were to move away from the centre of mass with increasing eccentricity due to the degradation of members on that side, the curve of fig.6.4 would show increasing downward slope. Chapter 6. CASES STUDIED 34 6.2 Structure 5: 3-Storey Structure with Weaker Columns According to the present design philosophy, buildings should be designed such that the ductility demands on the beams are much greater than those on the columns. By doing so, yielding will occur in the beams rather than in the columns. This is desirable because a failure in the beams is less catastrophic than a failure in the columns. In the previous cases, the columns may have been excessively strong. To examine the behaviour of structures with weaker columns, a 3-storey structure was designed based on an elastic dynamic analysis, and the member strengths were obtained by dividing the resulting bending moments by three and four for columns and beams, respectively. The design strengths of the frames are shown in fig.4.3. The layout of the structure is identical to that of the structure studied previously. Again, the El Centro earthquake scaled to a maximum acceleration of 0.3 g was used to excite the building in the y-direction. 6.2.1 D R A I N T A B S ( Takeda model) The results obtained from DRAINTABS using the Takeda model are shown in fig.6.5(a). The following observations can be made: 1. The frequency of oscillations was very high, and there were as many excur-sions to the left as there were to the right, implying that the torsional and the translational motions were often out of phase. 2. For all floors, the maximum excursion was to the right of the building, and it occured at 1.8 seconds; four-tenths of a second before peak excitation. The third floor had the maximum displacement followed by the first and then the second. Chapter 6. CASES STUDIED 35 3. The lengths of the excursions after the first 2.5 seconds were considerably shorter , possibly due to the fact that the strengths of the members were substantially reduced by the cyclic loading. 4. An inspection of the member rotations shows that two beams and two columns in frame 3 had actually collapsed. This explains why the centre of oscillations shifted to the right of the building approximately after the first two seconds. The amount of shifting was much greater for the second and third floors. 6.2.2 D R A I N T A B S ( E l a s t o - P l a s t ic m o d e l ) Results obtained using the Elasto-Plastic model are shown in fig.6.5(b). The following observations can be made: 1. In general, the behaviour of the centre of mass was less oscillatory than that predicted by the Takeda model. There were much fewer excursions, but they were generally longer than those predicted by the Takeda model. 2. Maximum displacement of the centre of stiffness occurred at 4.6 sec. for the second and third floors, and they were to the left of the building. However, for the first floor, maximum excursion occured at a much earlier instant, at 2.2 sec, and it was to the right of the building. 3. During the first second into the nonlinear range, the excursions were mostly to the right of the building, indicating that there was more yielding in the frames on the left hand side. However, with increasing length of excitation, the centre of stiffness began to make more excursions to the left. In fact, for the third floor, the centre of oscillations shifted to the left of the initial position indicating the some of the members on the right hand side had Chapter 6. CASES STUDIED 36 actually collapsed. This is contrary to expectation; one would expect the members on the same side as the centre of mass to suffer major damage. 6.2.3 P I T S A When attempting to use PITSA to analyse the structure, a convergence problem was encountered, and the program failed to converge after twenty one iterations. This problem will be further investigated in the next chapter. 6.3 Structure 2: Three-Storey Structure with Identical Frames The layout of this structure is the same as the previous one. However, the frames which are parallel to the y-axis are made identical in strength such that they all have the same strengths as frame 2 of structure 1. The same acceleration record was used to excite the building in the y-direction. 6.3.1 D R A I N T A B S (Takeda-Model) The results for this structure are shown in fig.6.6(a), and we can conclude the following: 1. The excursions of the centre of stiffness were predominantly to the left of the building. This was expected because, due to the eccentricity of the building, the frames on the right hand side received more seismic loading, and as they are identical in strength, there should be more yielding in these frames, causing the centre of stiffness to move away from them. 2. Maximum excursion occurred soon after peak excitation at 2.4 sec. for the second and the first floors. On the other hand, maximum excursion took place at a much later instant at 4.9 sec. for the first floor. Chapter 6. CASES STUDIED 37 3. The second floor had the maximum displacement of the centre of stiffness, followed by the first and the third. 6.3.2 DRAINTABS(Elasto-Plastic Model) The results are plotted in fig 6.6(b), and the following observations can be made: 1. Except for a few small peaks, the excursions of the centre of stiffness were almost exclusively to the left, suggesting the extensive yielding in the frames on the right hand side. 2. As before, the number of excursions was much smaller than that predicted by the Takeda model, but they were generally longer than those obtained by the Takeda model. 3. Maximum excursion occurred simultaneously at 4.6 sec. for all floors; 2.4 seconds after the occurrence of peak excitation. 4. The length of maximum excursion decreased with higher floors. 6.3.3 P I T S A Results from both PITSA and DRAINTABS are tabulated in table 6.2. Again the Takeda model was chosen in DRAINTABS. As found before, the maximum displace-ments obtained from PITSA were substantially smaller than those from DRAINTABS. Moreover, DRAINTABS predicted the second floor to have the maximum displacment, whereas PITSA predicted the maximum displacement to occur on the first floor. 6.3.4 Effect of Eccentricity Again DRAINTABS was used to obtain a plot of maximum displacement of the centre of stiffness vs. eccentricity. The curve is shown in fig.6.7. As expected, since the building Chapter 6. CASES STUDIED 38 is symmetric in strength, the centre of stiffness remained stationary throughout the entire excitation period at zero eccentricity. With increasing eccentricity, the centre of stiffness gradually moved further towards the left in a rather smooth and linear fashion; but it tended to concave upward, which was unexpected. 6.4 Structure 6: Three-Storey Structure with Asymmetrical Stiffness Dis-tribution The layout and member strengths of this structure is the same as structure 2 (3-storey framed structure with identical frames). However, in this case, the stiffness of frames 3 and 4 was increased by 15 % so that the initial position of the centre of stiffness was to the left of the centre of mass which was placed at the centre of the building. The reason for chosing this structure was that by making the structure unsymmetrical w.r.t. stiffness, it was hoped that the translational and rotational motions would be coupled, and thus preventing the centre of stiffness from moving to the right where the frames were softer. The time history plot of the centre of stiffness is shown in fig.6.8. Since the centre of mass was to the right of the initial position of the centre of stiffness, the frames on the right should be more heavily loaded, and it was therefore expected that the centre of stiffness would move predominantly to the left. However, contrary to what was expected, the excursions were mostly to the right, although the maximum displacement of the centre of stiffness occurred on the left hand side at 4.58 seconds for all floors. To investigate further the behaviour of the centre of stiffness, the time history for the first floor, together with the displacement and rotation of the first floor diaphragm were all plotted on the same graph as shown in fig.6.9. The initial position of the centre of stiffness was to the left of the centre of mass Chapter 6. CASES STUDIED 39 at 343 inches. There were three major excursions in the plus x-direction, occurring approximately at 1.76, 2.13 and 2.75 seconds, respectively. As can be seen in fig.6.9, each of these three instances occurred when the displacement and the rotation of the diaphragm became out of phase at the instant of maximum excursion. An examination of the yield codes of the members in frames 3 and 6 shows that during the descending portion of the rotation curve, the rate of increase of the deflection of frame 6 decreased while that of frame 3 increased due to the rotation of the diaphragm. When the decrease in deflection of frame 6 due to rotation of the diaphragm could no longer be offset by the increase due to the diaphragm displacement, frame 6 would actually start to unload. On the other hand, the deflection of frame 3 kept on increasing causing more yielding of its members. Since the stiffness of a member during yielding is always smaller than the stiffness during unloading, the left hand side of the structure became softer relative to the other side causing the centre of stiffness to move to the right. When the displacement of the diaphragm reached its peak, it began to move back towards the centre, and frame 3 would eventually begin to unloaded. As all the frames were unloading together, the centre of stiffness began to move back to its original position because, according to the Takeda model, the unloading stiffness of a member is almost the same as its initial elastic stiffness. It is these characteristics of the Takeda model that govern the entire history of the motion of the centre of stiffness. The maximum displacement of the centre of stiffness occurred at 4.58 seconds on the left hand side of the structure. From fig.6.8, it can be seen that at that instant, the displacement and the rotation of the diaphragm peaked simultaneously, thus confirming the hypothesis that the direction of the excursion of the centre of stiffness depends on the relative phase between the translational and rotational motions, and hence the positions of members in these hysteretic cycles and thus their stiffnesses. As the excursions in this case were predominantly to the right of the structure, the Chapter 6. CASES STUDIED 40 strength of the frames on the left must have degraded more than those on the left, and therefore it was expected that the centre of oscillations would gradually move away from its original position towards the centre of mass. This was confirmed by the many plateaus on the time history plot of the centre of stiffness. 6.5 Conclusions From all the cases studied, the following conclusions can be drawn: • The movement of the centre of stiffness is dependent on the hysteretic model used. • The behaviour of the centre of stiffness predicted by the elasto-plastic model is always less oscillatory than that predicted by the Takeda model. However, the length of the excursions predicted by the former is usually longer. • The maximum displacement of the centre of stiffness predicted by the elasto-plastic model usually occurs at a much later instant than that predicted by the Takeda model. • The relationship between the maximum displacement of the centre of stiffness and the eccentricity seems to be fairly linear. • Although PITSA has been proven to be a reliable tool in predicting the damage ratios of members under seismic excitation, it does not seem to be able to predict the maximum displacement of the centre of stiffness with any acceptable accuracy. • The direction of the excursion of the centre of stiffness depends on whether or not the translational motion is in phase with the rotational motion. Chapter 6. CASES STUDIED 41 • W i t h increas ing length of e x c i t a t i o n , the centre of stiffness w i l l g r a d u a l l y m o v e away f r o m the s ide tha t has y i e l d e d m o r e extensively . However , th is m o v e m e n t was f o u n d not to be s igni f icant . C h a p t e r 7 C O N V E R G E N C E P R O B L E M S I N U S I N G P I T S A As noted in chapter 6, a convergence problem was encountered when PITSA was used to analyse structure 5. The main purpose of this chapter is to investigate this problem, and to examine the effectiveness of two routines which are supposed to expedite convergence. 7.1 C o n v e r g e n c e S p e e d i n g R o u t i n e While the modified substitute structure method has been shown to be a reliable tool of analysing the inelastic behaviour of 3-dimensional structures, its convergence rate and therefore its effectiveness depends very much on the type of structure being analysed. Early computer runs using this method showed that for some structures, the damage ratios converged very slowly, or oscillated about the final answer. The original modified substitute structure method program contained a routine which proved effective in arriving at the final answer for those cases in which the changing damage ratios were either decreasing or increasing steadily. This routine operated on the basis of adding to the damage ratios that were to be returned to the main program a factor multiplied by the change in the damage ratios over the last iteration. In this manner, changing damage ratios were moved faster in a direction which hopefully was towards the true answer. This routine is effective in cutting down the number of iterations required for cases where the damage ratios move in one direction. However, in those cases where the damage ratios oscillate, the routine would actually act as a deterrent to convergence. Metten5 in his master's thesis modified the routine such that it could better establish 42 Chapter 7. CONVERGENCE PROBLEMS IN USING PITSA 43 the damage ratio trends by keeping track of damage ratios from more than just the last iteration. At each iteration, three damage ratios are stored, the current one plus those from the last two iterations. The nine possible trends in the ratio are then determined and those ratios that seem to be consistently decreasing or increasing will be modified by appropriately adding or subtracting a factor multiplied by the difference of the last two values as in the original problem. On the other hand, in the case of oscillating damage ratios, the oscillations are damped into producing an answer lying between the last two values. In cases where the damage ratios do not change significantly over two oscillations, no modification will be made. The above scheme is summarized in fig.7.1. the factor f3 is a positive number less than unity. A value of zero effectively shuts off the routine, while a value greater than one will cause divergence. During the first few wildly changing iterations, the convergence speeding routine is shut off to let the modified substitute structure method program naturally home closer to the final answer. The value of /? is usually chosen arbitrarily. There is no optimum value for all structures as different structures converge to the right solution in very different ways. Metten5 found in general that the convergence rate increases with the value of beta. The convergence speeding routine has been incorporated into PITSA. It is activated after the first 8 iterations, and beta is assigned a value of 0.8. The results for structure 5 obtained from PITSA with the routine added are shown in fig.7.2. It can be seen that the damage ratio errors and the 1 st mode natural periods oscillated wildly during the first 25 iterations, and there was no indication of the program converging to the final solution. To examine the effect of the convergence speeding routine, a second run was made ) on the same structure with the routine shut off. The results are shown in fig.7.3. Surprisingly, without the convergence speeding routine, the amplitudes of oscillations of both the damage ratio errors and the periods decreased steadily. However, after 30 Chapter 7. CONVERGENCE PROBLEMS IN USING PITSA 44 iterations, the damage ratio error still failed to satisfy the tolerance limit set to be 0.01. Nevertheless, it was expected that the convergence criterion would eventually be met within the next 10 iterations. Therefore, the convergence speeding routine actually acted as a deterrent to convergence in this case. As shown in fig.7.3, the program was very close to convergence after 23 iterations. To see if the damage ratio error tolerance is too stringent, it was increased from 0.01 to 0.05. As expected, the program converged after 19 iterations, and the damage ratios were compared to those obtained from DRAINTABS. As can be seen in fig.7.5, with the exception of frame 3, the agreement between the two programs was surprisingly good. PITSA was able to predict the distribution of damage correctly. For the beams, the damage ratios predicted by PITSA were higher in most cases, whereas for the columns, the predictions were less conservative. The average discrepancy was approximately 15 %. For frame 3, results from DRAINTABS showed that columns 5 and 6, and beams 2 and 3 underwent very large rotations which implied that they had actually failed. On the other hand, the damage ratios predicted by PITSA for those members were relatively small. Although PITSA failed to detect the failure of a number of members, it was probably caused by the fact that the structure was poorly designed, and not by the relaxation of the convergence criterion. It was initially thought that the convergence problem was caused by those members in frame 3 which had collapsed. However, an investigation of the source of the damage ratio error in each iteration showed that this was not the case. Since PITSA predicted the ductility requirements very well for the other frames, it is reasonable to suggest that the damage ratio error tolerance limit can be reduced for structures which require a great number of iterations to converge. The extent to which it should be reduced depends on the accuracy of results desired. Chapter 7. CONVERGENCE PROBLEMS IN USING PITSA 45 7.2 The Solution Searching Routine This routine was developed by Hui 6 in his thesis. It was claimed to be a better procedure for expediting the convergence of the modified substitute structure when the natural period of the substitute structure lies within the rapidly descending portion of the response spectrum. This routine is especially useful when the natural period is in the vicinity of a corner in the spectrum. Details of this routine are given in Hui's thesis and will not be repeated here. The procedure is based on the premise that the solution of any problem using the modified substitute structure method is actually bounded by the values of spectral acceleration corresponding to the early stages of the iterative procedure. Therefore, these values can be used as the upper and lower bounds that will limit the range of the capacity curve where the solution has to be found. The procedure is outlined below:6 1. To detect a convergence problem: after the first five iterations, start to compare the changes in first mode response period with the previous two iterations and determine whether reversal has occurred in the last three iterations. 2. Store the values of the first mode spectral acceleration for the previous iteration. This will be used as an upper or lower bound for the search once a convergence problem is detected. 3. When a convergence problem is detected, use a control flag to turn on the search routine in place of parts of the normal procedure. 4. Use the binary search method to search for the solution. During each trial, a constant first mode spectral acceleration is used until the program converges in the normal way. Chapter 7. CONVERGENCE PROBLEMS IN USING PITSA 46 5. Then check this acceleration against the acceleration from the spectrum, for the same period and the same level of damping. If the values do not agree then repeat step 4 with another trial value of constant acceleration. If the values do agree within a reasonable limit of tolerance, then the program has actually converged. The above procedure is performed by a subroutine named STACHK, a listing of which is given in the appendix. It should be noted that STACHK only works for the first mode of the analysis. It will detect any reversal of the response period after the fifth iteration. If the reversal is larger than a predefined limit, the search routine will be turned on and a trial value of spectral acceleration fed into the program. Once convergence is attained, this acceleration is checked against that of the spectrum, at the same level of damping, using the following criterion: So>cai ~ Sasp SADIF = where SCbgp < 0.015 (7.1) Sacai = the trial value of spectral acceleration at the current level of damping SaBp = the spectral acceleration determined from the spectrum at the same level of damping If this criterion is not met, then the upper or the lower bound has to be adjusted according to these rules: 5 < = San-! for Sacal > Sasp (7.2) Saln = So„_! for Sacal < SaBp (7.3) Chapter 7. CONVERGENCE PROBLEMS IN USING PITSA 47 where San-! the trial value of spectral acceleration at 2 % damping from the last iteration the new upper bound at 2 % damping Sa\ 'n the new lower bound at 2 % damping (7.4) The next trial will again be the mean of these new bounds. This will continue on until eqn.7.1 is satisfied. A flow chart of the subroutine STACHK is shown in fig.7.6. To find out if the solution searching routine can improve the convergence rate for the present problem, STACHK was incorporated into PITSA, and the same structure analysed. The results are shown in fig.7.4. A convergence problem was detected at iteration 11, and the solution searching routine was activated. When compared to fig.7.2, it can be seen that with STACHK in effect, the magnitude of oscillations was substantially reduced, and the program seemed to be converging to the correct answer, although a vast number of iterations would be required. In order to examine the effectiveness of the convergence speeding routine, beta was set to zero and the analysis repeated. Surprisingly, the results obtained were exactly the same as before. 7.3 C o n c l u s i o n s • When using the modified substitute structure method, the convergence speed-ing routine seems to be an effective tool in reducing the number of iterations for structures whose damage ratios converge to the final answer monotonically. However, for cases where the damage ratios tend to oscillate, the convergence Chapter 7. CONVERGENCE PROBLEMS IN USING PITSA 48 speeding routine will actually act as a deterrent to convergence. • The solution searching routine was found to be effective in reducing the magnitude of oscillations for problem cases, and therefore decreasing the number of iterations required for convergence. • The damage ratio error tolerance limit set in PITSA was found to be stringent and could be relaxed for poorly designed structures which require a large number of iterations to converge. The extent to which it should be reduced depends on the degree of accuracy desired. Chapter 8 CONCLUSIONS Based on the results obtained from this study, the following conclusions can be made: • The behaviour of the centre of stiffness of eccentric structures in the inelastic range is generally very oscillatory. Excursions to both sides of the initial position are possible because the displacement and the rotation of the structure become out of phase occasionally at instant of maximum excursion. When this happens, the side which is more heavily loaded due to torque created by the eccentricity will start to unload while the other side continues yielding. If the unloading stiffness is very close to the elastic stiffness in the hysteretic model, the unloading side becomes stiffer causing the centre of stiffness to move towards it. When both sides are unloading, the centre of stiffness will start moving back to its original position. • With increasing length of excitation, the eccentricity of the structure will increase due to the gradual shifting of the centre of stiffness away from its original position towards the side whose strength has been degraded by reverse loading. However, the amount of increase was found to be small. Nevertheless, it is worth noting that this finding was based on the assumption that reinforced concrete members under cyclic loading will behave in accordance with the Takeda model. Very different behaviour of the centre of stiffness is expected when a different model is used. Unless the real behaviour of reinforced concrete members in the inelastic 49 Chapter 8. CONCLUSIONS 5 0 range is known, it cannot be concluded with certainty that the magnification of the initial eccenticity of a structure during an earthquake is insignificant. • Although PITSA is a reliable method of predicting the damage patterns of build-ings under seismic excitation, it fails to predict the maximum displacement of the centre of stiffness accurately. Figure 2 .1 : Definition of Damage Ratio 52 Figure 2.2: Plan view of floor showing frame and diaphragm horizontal displacement (PITSA) (from ref. 3). 53 o . o « 1 1 1 J 1 t . , ; • 0.5 1.0 1.5 2.0 2.5 3.0 1 2 '0 8 6 4 2 Period, sec Frequency, hertz S p e c t r a l A c c e l , for = 8 S p e c t r a l A c c e l , tor 0=0.02 6+100-/3 Figure 2.3: Spectrum A (from ref. 3) 54 Figure 3.1: Diaphragm and frame horizontal displacements(DRAINTABS)(from ref.3) Member data : Columns Beams Area (in2) 576 495 Cracked Moment of Inertia (in4) 27648 44921 E (ksi) 3000 3000 Strain Hardening Ratio 0 0 Diaphragm Data : Weight in x-direction = 450 kips Weight in y-direction = 450 kips Rotary Moment of Inertia = 228750 kip.ft2 Seismic Excitation : El Centro NS scaled to a maximum acceleration of 0. Eccentricity, e : 54 ins. Figure 4.1: (a) Member properties of Structures 1,2 and 3. 56 -®-C M 3 @ 20' 8 12 16 4 8 12 7 11 15 3 7 11 i—1 6 10 14 ~~-ce O FLOOR AREA - 100'x60 CORE AREA - 30'X18' Figure 4.4: Dimensions and properties of structure 3 I 62 Member data : Columns Beams Area (in2) 1000 700 Cracked Moment of Inertia (in4) 42000 26000 E (ksi) 3000 3000 Strain Hardening Ratio 0.0001 JD0Q1 Yield Moment (kip.in) 18000 3600 Diaphragm Data : Weight in x-direction = 360 kips Weight in y-direction = 360 kips Rotary Moment of Inertia = 156250 kip.ft2 Seismic Excitation : El Centro NS scaled to a maximum acceleration of 0.3 g. Eccentricity, e : 48 ins. Figure 4.5: (a) Member properties of structure 4. 63 4 CM-1 3 @ 2() = 60' •* Figure 4.5: (b) Layout and dimensions of structure 4 0.4^ 0.3H 0.2-\ cn c o "o Q) CD O O < 0.0 - 0 . H - 0 . 2 H -0.3-f 20 30 Time(sec) Figure 4.6: El Centro NS earthquake record, 1940. 65 Figure 4.7: (a) Damage ratios for structure 1. Figure 4.7: (b) Damage ratios for structure 1. 67 Figure 4.7: (c) Damage ratios for structure 1. 68 Figure 4.8: (a) Damage ratios for structure 2. Figure 4.8: (b) Damage ratios for structure 2. 70 Figure 4.8: (c) Damage ratios for structure 2. Legend U2 U CtHTRt NS B. COfTSOCW COLUMNS COLUMNS Figure 4.9: Damage ratios for structure 3. 72 O . S F R A M E 1 Figure 4.10: (a) Damage ratios for structure 4. 73 F R A M E 1 B E A M S 5 F R A M E 4-B E A M S Figure 4.10: (b) Damage ratios for structure 4. F R A M E 1 cr> O.S -I I i 2 3 -* a e 7 e 9 i o n 12 13 I A i s C O L U M N S V > 0.3 H *3 L e g e n d CZ3 O R T B S 1 2 3 C O L U M N S I O n 1 2 1 3 n Figure 4.10: (a) Damage ratios for structure 4. 75 Figure 5.1: Locating the centre of stiffness from the condensed stiffness matrix 20 ' >-Plan @ I N Elevation Elevation weight of diaphragm = 60 Kips rotary moment of inertia = 4000 Kips.ft2 Beams Columns EI (K.ft2) 0.125 x 106 0.2083 x 106 EA (K) 0.900 x 106 0.1500 x 107 Figure 5.2: Layout and member data of test structure 13-12- 1 st FLOOR or ii-o ~o i o -o o a > 9 -Max. Displ.= 2.11 ft. @ 2.66 sec. 13-100 2 0 0 3 0 0 4 0 0 12- 2 nd FLOOR CH 11-o tn X) 10-o o o I > 9 -Mox. Displ.= 2.72 ft. @ 2.66 SEC. 100 2 0 0 3 0 0 Time Steps 4 0 0 5 0 0 Figure 5.3: Time history plot of CR for test structure 78 Figure 6.1: Hinge moment-rotation relationship for Takeda model (from ref. 7) 79 Figure 6.2: Elasto-plastic moment-curvature relationship 7 0 0 -6 0 0 -5 0 0 -TAKEDA Model 1 st FLOOR Max. Displ. = 106" @ 2.36 sec . 100 2 0 0 3 0 0 4 0 0 2 nd FLOOR Max. Displ. = 166" @ 2.36 sec . —i— 100 2 0 0 3 0 0 4 0 0 3 rd FLOOR 4 0 0 -3 0 0 -Max. Displ. = 189" [ ' @ 2.37 sec . T T 100 2 0 0 3 0 0 Time Steps 4 0 0 5 0 0 Fig.6.3(a)-Time History Plot of CR for Structure 1 7 0 0 -6 0 0 -5 0 0 -ELASTO-PLASTIC 1 st FLOOR 4 0 0-1 3 0 0 -2 0 0 -Max. Displ. = 2 2 0 " @ 4 .56 sec. 100 2 0 0 3 0 0 400 2 nd FLOOR r~v~ Max. Displ. = 2 2 3 " @ 4 .56 sec. 100 2 0 0 3 0 0 4 0 0 3 rd FLOOR Max. Displ.= 274" @ 4 .56 sec. 1 i 1 i 2 0 0 3 0 0 Time Steps 100 4 0 0 500 Fig.6.3(b)-Time History Ploi of CR for Structure 1 82 Figure 6.4: Graph of maximum displacement of CR vs eccentricity for structure 1. 83 Max. Displ. = 269" @ 1.85 sec. TAKEDA Model 1 st FLOOR 100 200 Max. Displ. = 233" @ 1.85 sec. 300 400 2 nd FLOOR 100 Max. Displ.= 312" @ 1.81 sec. 200 300 400 3 rd FLOOR — r -100 200 300 Time Steps 400 500 Fig.6.5(a)-Time History Plot of CR for Structure 5 Max. Displ. = 303" @ 2.19 sec. T P lr L / I M J 1 st FLOOR ELASTO-PLASTIC 100 200 300 400 V 2 nd FLOOR U L_i Max. Displ. = 327" @ 4.57 sec. 100 200 300 400 3 rd FLOOR Max. Displ.= 386" @ 4.57 sec. 100 200 300 Time Steps 400 50 Fig.6.5(b)-Time History Plot of CR for Structure 5 85 7 0 0 -«2 6 0 0 -o 5 0 ° -TAKEDA Model 1 st FLOOR O c i_ o o O I X 4 0 0 -3 0 0 -2 0 0 -1 0 0 -Max. Displ. = 163" <§> 4.86 sec. 8: i 100 70 « 6 0 0 -2 0 0 3 0 0 4 0 0 2 nd FLOOR o 5 0 0 ' in a> "5 c o o (J I X 4 0 0 -3 0 0 -2 0 0 -1 0 0 -Max. Displ. = 192" @ 2.36 sec. 8=-7 0 «2 6 0 0 -100 2 0 0 3 0 0 4 0 0 3 rd FLOOR g 5 0° -in "o _C T» o o o I X 4 0 0 -3 0 0 -2 0 0 -100 - Max. Displ.= 124" @ 2.37 sec. 100 i • i 2 0 0 3 0 0 Time Steps 4 0 0 5 0 0 Fig.6.6(a)-Time History Plot of CR for Structure 2 Figure 6.7: Graph of maximum displacement of CR vs eccentricity for structure 2. TAKEDA-Model i 100 I 2 0 0 1 st FLOOR -v-Max. Displ. = 64" @ 4.56 sec . 3 0 0 4 0 0 2 nd FLOOR Max. Displ. = 112" @ 4.56 sec . 100 2 0 0 3 0 0 4 0 0 3 rd FLOOR Max. Displ.= 151" @ 4.56 sec . —i 5 0 0 100 2 0 0 3 0 0 Time Steps 4 0 0 Fig.6.8-Time History Plot of CR for Structure 6 Fig.6.9 Plots of the Displacement, Rotation and CR of the First Diaphragm of Structure 6 oo CO l<0 / DAMOLD = Damage Ratio for i-2 Iteration DAMB = Damage Ratio for i-1 Iteration DAMRAT = Damage Ratio for i Iteration DAMDIF = DAMRAT-DAKB DR = Damage Ratio returned to program, Figure 7.1: Flow-chart for convergence speeding routine (from ref. 2) NO STACHK BETA=0.0(1-8) BETA=0.8(9- ) Fig.7.2 Plots of Damage Ratio Error and Natural Period for Structure 5. NO STACHK Beta=0.0 1 s i M o d e ' i i i i i i 2 6 10 14 18 22 26 Iteration Number FIG.7.3 Plots of Damage Ratio Error and Natural Period for Structure 5 30 CO ts5 FIG.7.4 Plots of Damage Ratio Error and Natural Period for Structure 5 CD w 94 Figure 7.5: (a) Damage ratios for structure 5. 95 Figure 7.5: (b) Damage ratios for structure 5. 96 Figure 7.5: (c) Damage ratios for structure 5. 97 s e t u p p e r a n d 1ower b o u n d Sa ( i s 2% damp i n g ) s e t I F L A G = 0 c a l c u l a t e mean Sa a s g u e s s (? 2% damp i n g ) m o d i f y Sa w i t h a p p r o p r i a t e d a m p i n g RETURN c h e c k c o n v e r g e n c e : c o m p a r e Sa w i t h Sa f r o m s p e c t r u m y e s PRINT t r i a l n o . i n b i n a r y s e a r c h & SADIF c a l c u l a t e t h e d i f f e r e n c e b e t w e e n two S a ' s ( S A D I F ) LOCK : 0 = f o l l o w n o r m a l r o u t i n e 1 • f o l l o w s e a r c h r o u t i n e 2 - s e a r c h c o m p l e t e d I FLAG : 0 - n o r m a l i t e r a t i o n 1 = l a s t i t e r a t i o n ICOUNT : i t e r a t i o n n o . a d j u s t l o w e r o r u p p e r b o u n d S a ( «• 2'/. d a m p t n g ) PRINT c a l c u l a t e Sa a c t u a l Sa & SA D I F s e t I FLAG = 0 s e t LOCK ' 2 RETURN Figure 7.6: Flow-chart for solution searching routine STACHK (from ref. 6) Chapter 8. CONCLUSIONS 98 Floor PITSA DRTBS (Takeda) 1 st 76" 106" 2 nd 63" 166" 3 rd 56" 189" Table 6.1: Comparison of maximum displacements of CR predicted by DRAINTABS and PITSA for structure 1. Floor PITSA DRTBS (Takeda) 1 st -33" -163" 2 nd -27" -192" 3 rd -25" -124" Table 6.2: Comparison of maximum displacements of CR predicted by DRAINTABS and PITSA for structure 5. Bibliography [l] Yoshida, S,"Modified Substitute Structure Method for the Analysis of Existing R/C Structures," Master's Thesis U.B.C., March 1979. [2] Metten, A.W.F.,"The Modified Substitute Structure Method As A Design Aid for Seismic Resistant Coupled Structural Walls," Master's Thesis U.B.C., March 1981. [3] Tarn, K.S.,"Pseudo Inelastic Torsional Seismic Analysis Utilizing the Modified Substitute Structure Method," Master's Thesis U.B.C., April 1985. [4] Guendelman-Israel, R., and Powell, G.H., "DRAIN-TABS: A Computer Program for Inelastic Earthquake Response of Three Dimensional Buildings," Report No.77-8, University of California, Berkeley, March 1977. [5] Cansius, T.D.G., "Inelastic Torsional Seismic Analysis of FRame Buildings", Mas-ter's Thesis U.B.C., Jan 1986. [6] Hui, L., "Pseudo Non-linear Seismic Anaysis," Master's Thesis U.B.C.,October 1984. [7] Kanaan, A., and Powell, G.H., "General Propose Computer Program for Inelastic Dynamic Response of Plane Structures", Report No. EERC 73-6, University of California, Berkeley,1973. [8] Tso, W.K., and Sadek, A.W., "Inelastic Response of Eccentric Structures", 4 th CCEE Vancouver, Canada, June 1983, pp.261-270. 99 Appendix A Listing of Subroutines Added to D R A I N T A B S 100 101 4320 SUBROUTINE TRANS(M,A,ABF.NEOB,NLEVEL.I S T E P ) 4320 5 C • • 4320 7 C * THIS SUBROUTINE EXPANDS THE ONE-DI MENS IONAL MATRIX * 4320 8 C • CONTAINING THE BUILDING S T I F F N E S S MATRIX IN COMPACTED * J 3 2 0 9 c * FORM INTO A FULL TWO-DIMENSIONAL ARRAY. 4320 95 c * 4320 97 c * INPUT: -AOORESSES OF ALL DIAGONAL ELEMENTS STORED IN * 4320 98 c M ( I ) . * 4320 99 c -ALL NON-ZERO ELEMENTS IN THE UPPER HALF OF THE * 4320 995 c * S T I F F N E S S MATRIX. STORED IN A ( I ) . * 4320 997 c * - S I Z E OF THE BUI L D I N G S T I F F N E S S MATRIX. NEOB. * 4 3 2 0 998 c * OUTPUT:-FULL B U I L D I N G MATRIX. A B F ( I . d ) . * 4320 999 c 4321 DIMENSION M ( 1 ) . A ( 1 ) , A B ( 1 0 0 . 1 0 0 ) . K D I A ( 1 0 0 ) . A B S ( 1 0 0 . 1 0 0 ) . 4322 C A B F ( 1 0 0 . 1 0 0 ) 4323 K D I A ( 1 ) = 1 4324 N=1 4325 DO 8 0 1 0 1 = 1 .NEOB 4326 DO 8 0 2 0 u =1.NEOB 4327 AB( I , J ) = 0 . 0 4328 8 0 2 0 CONTINUE 4329 8 0 1 0 CONTINUE 4330 DO 8 0 4 0 1=1,NEOB 4331 IF ( M ( I ) . N E . K O I A ( I ) ) GOTO 8 0 3 0 4332 DO 8 0 5 0 0=1.1 4333 A B ( I . J ) = A ( N ) 4334 A B ( d . l ) = A ( N ) 4335 N=N+1 4336 8 0 5 0 CONTINUE 4337 KK=I+1 4338 K D I A ( K K ) = M ( I ) + K K 4339 GOTO 8 0 4 0 4340 8 0 3 0 CONTINUE 4341 N Z = K D I A ( I ) - M ( I ) 4342 K=NZ+1 4343 IF ( I . E Q . N E O B ) GOTO 8 1 1 0 4344 KK=I+1 4345 K D I A ( K K ) = M ( I ) + K K 4346 81 10 DO 8 0 7 0 d=K.I 4347 A B ( I . J ) = A ( N ) 4348 A B ( d . I ) = A ( N ) 4349 N=N+1 4350 8 0 7 0 CONTINUE 4351 8 0 4 0 CONTINUE 4356 IR=1 4357 DO 8 2 0 0 MM=1,NLEVEL 4358 DO 8 2 1 0 K=1 ,3 4359 L=(K-1)*NLEVEL+MM 43 6 0 DO 8 2 2 0 1=1.NEOB 4361 A B S ( I R . I ) = A B ( L . I ) 4362 8 2 2 0 CONTINUE 4363 1R=IR+1 4364 8 2 1 0 CONTINUE 4365 8 2 0 0 CONTINUE 4366 !C=1 4367 00 830O MM=1,NLEVEL 4368 0 0 8 3 1 0 K=1.3 4369 L=(K-1)*NLEVEL+MM 102 4 3 7 0 DO 8 3 2 0 1=1,NEOB 437 1 A B F ( I , I C ) = A B S ( I , L ) 4372 8 3 2 0 CONTINUE 4373 IC=IC+1 4374 8 3 1 0 CONTINUE 4375 8 3 0 0 CONTINUE 4 380 RETURN 438 1 END 4382 SUBROUTINE C 0 N D E N ( A B . N E O B . N L E V E L . I F L . A K , I S T E P ) 4 3 8 2 . 5 C * 4 3 8 2 . 7 C * THIS SUBROUTINE CONDENSES THE B U I L D I N G S T I F F N E S S MATRIX * 4382 . 8 C * INTO A 3X3 MATRIX CONTAINING ELEMENTS WHICH BELONG TO THE * 4382 . 9 C * FLOOR FOR WHICH THE CENTRE OF R E S I S T A N C E IS TO BE LOCATED. * 4382 . 95 C r + 4382 . 97 C • INPUT : -BUILDING S T I F F N E S S MATRIX, A B ( I . d ) * 4 3 8 2 . 98 C -NUMBER OF THE FLOOR OF I N T E R E S T , I F L * 4382 . 99 C * OUTPUT: -CONDENSED FLOOR S T I F F N E S S MATRIX. A K ( I . d ) * 4382 . .995 C A * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 4383 OIMENSION AB( 100, 1 0 0 ) , A K 1 ( 1 0 0 , 1 0 0 ) , A K 2 ( 100, 1 0 0 ) , A K 3 ( 1 0 0 , 1 0 0 ) , 4384 C A K 4 ( 1 0 0 . 1 0 0 ) . A A ( 2 0 0 0 ) . A I ( 1 0 0 , 1 0 0 ) . P 1 ( 1 0 0 , 1 0 0 ) , P 2 ( 1 0 0 , 1 0 0 ) , 4385 C A K ( 3 , 3 ) 4385 . 2 W R I T E ( 1 2 , 9 2 1 0 ) I F L 4385 . . 3 92 10 FORMAT(///, '* » * * CONDENSED S T I F F N E S S MATRIX FOR FLOOR' , 1 4 , 3 X , 4385 . . 4 C '****'./) 4 3 8 5 . 5 C * * • * NO CONDENSATION REQUIRED IF THE STRUCTURE CONSISTS OF 4 3 8 5 . , 7 C ONE FLOOR ONLY * * * « 4386 IF ( N L E V E L . E O . 1 ) GOTO 9 0 0 0 4 3 8 6 . .5 C * • * * INTERCHANGE ROWS AND COLUMNS SO THAT THE 3X3 SUBMATRIX IN 4 3 8 6 . .7 C THE UPPER L E F T CORNER CONTAINS ALL ELEMENTS BELONGING TO 4 3 8 6 . .8 C THE ( I F L ) t h f l o o r * * * * 4387 I I = ( I F L - 1 ) * 3 + 1 4388 DO 9 0 1 0 1=1,3 4389 DO 9 0 2 0 d =1,NEOB 4 3 9 0 A = A B U ,d) 4391 AB ( I , d)=AB( I I ,d) 4392 A B ( 1 1 . d ) = A 4 3 9 3 9 0 2 0 CONTINUE 4394 1 1 = 1 1+1 4395 9 0 1 0 CONTINUE 4396 I I = ( I F L - 1 ) * 3 + 1 4397 DO 9 0 3 0 d= 1 . 3 4398 DO 9 0 4 0 1=1,NEOB 4399 A = A B ( I . d ) 4 4 0 0 A B ( I , d ) = A B ( I , I I ) 4401 A B ( I . I I ) = A 4402 9 0 4 0 CONTINUE 4403 1 1 = 1 1 + 1 4404 9 0 3 0 CONTINUE 4404 .5 C * * • * START PA R T I T I O N I N G OF THE B U I L D I N G MATRIX * • • * 4 409 00 9 0 5 0 1=1.3 44 10 DO 9 0 6 0 J=1.3 441 1 A K 1 ( I . d ) = A B ( I , d ) 4412 9 0 6 0 CONTINUE 44 15 9 0 5 0 CONTINUE 44 16 DO 9 0 7 0 1=1.3 44 17 M=1 M 1 8 DO 9 0 8 0 d=4,NEQB A K 2 ( I , M ) = A B ( I , d ) 103 4 4 2 0 M=M+1 4421 908O CONTINUE 4424 9 0 7 0 CONTINUE 4425 L=1 4426 0 0 9 1 0 0 1=4,NEOB 4427 DO 9 1 1 0 d= 1 . 3 4428 AK3 ( L , J ) = A B ( I , d ) 4429 9 1 1 0 CONTINUE 4 4 3 0 L = L+1 4431 9 1 0 0 CONTINUE 4436 L=1 4437 DO 9 1 2 0 1=4.NEOB 4438 M=1 4439 DO 9 1 3 0 d=4.NEQB 4 4 4 0 AK4( L , M)=AB(I , d) 4441 M=M+1 4442 9 1 3 0 CONTINUE 4443 L = L+1 4444 9 1 2 0 CONTINUE 4449 N=NE0B-3 4 4 5 0 DO 9 1 4 0 1=1,N 4451 DO 9 1 5 0 J = 1.I 4452 I D = ( I * ( I - 1 ) ) / 2 + J 4 4 5 3 A A ( I D ) = A K 4 ( 1 , 0 ) 4456 9 1 5 0 CONTINUE 445 7 9 1 4 0 CONTINUE 4458 C A L L S Y M I ( A A . N , I E . N 1 ) 4 4 5 9 DO 9 1 6 0 1=1.N 4 4 6 0 DO 9 1 7 0 d= 1 ,N 4461 I F ( d . G T . I ) GOTO 9 1 6 0 4462 I D = ( I * ( I - 1 ) ) / 2 + d 4 4 6 3 A I ( I , d ) = A A ( I D ) 4464 A I ( J , I ) = A A ( I D ) 4 4 6 5 9 1 7 0 CONTINUE 4466 9 1 6 0 CONTINUE 44 7 3 C A L L G M U L T ( A K 2 , A I , P 1 , 3 , N . N , 1 0 0 , 1 0 0 , 1 0 0 ) 4474 C A L L G M U L T ( P 1 , A K 3 , P 2 , 3 , N , 3 , 1 0 0 , 1 0 0 , 1 0 0 ) 4475 DO 9 1 8 0 1=1,3 4476 DO 9 1 9 0 d=1 .3 4477 A K ( I , J ) = A K 1 ( I , J ) - P 2 ( I , d) 4478 9 1 9 0 CONTINUE 44 7 9 W R I T E ( 1 2 . 9 2 0 0 ) ( A K ( I . L ) , L = 1 . 3 ) 4 4 8 0 9 2 0 0 F 0 R M A T ( 3 E 1 2 . 4 ) 4481 9 1 8 0 CONTINUE 4 4 8 2 GOTO 9 4 0 0 4 4 8 3 9 0 0 0 CONTINUE 4484 DO 9 3 0 0 1=1.3 4 4 8 5 DO 9 3 1 0 d=1,3 44 8 6 A K ( I , J ) = A B ( I . d ) 44 8 7 9 3 1 0 CONTINUE 4 4 8 8 W R I T E ( 1 2 , 9 2 0 0 ) ( A K ( I . L ) . L = 1 , 3 ) 44 8 9 9 3 0 0 CONTINUE 4 4 9 0 9 4 0 0 CONTINUE 4491 RETURN 4492 END 44 9 3 SUBROUTINE C R ( A , C M X , C M Y . I F L , I S T E P , C R X M , C R Y M . N S T E P S . C R X F . C R Y F ) 4 4 9 3 . 5 C 4493 . 7 C * T H I S SUBROUTINE LOCATES THE CENTRE OF R E S I S T A N C E OF THE * 104 4493 . 4 4 9 3 . 4493 . 4493 . 4493 . 4 4 9 3 . 4493 . 4493 . 4493 . 4494 4494 . 4494 . 4495 4495 . 449G 4 4 9 6 . 4496 . 4 4 9 6 . 4 4 9 6 . 4497 4498 4498 . 4 4 9 8 . 4498 , 4498. 4499 4502 4503 4 504 4505 4506 4507 4508 4509 4510 451 1 4512 4513 4514 4515 4 5 1 5 . 4516 45 1 6 . 4 5 1 6 . 4 5 1 6 . 4517 4518 45 1 8 . 4 5 1 8 . 4 5 1 8 . 4518 . 4519 4522 4523 4524 4525 4526 4527 8 9 95 97 98 99 995 997 998 5 7 04 5 7 8 I F L t h FLOOR. * INPUT : -CONDENSED FLOOR S T I F F N E S S MATRIX * -NUMBER OF THE FLOOR OF I N T E R E S T . I F L * -X AND Y COORDINATES OF THE CENTRE OF MASS; CRX.CRY* OUTPUT:-X AND Y COORDINATES OF THE CENTRE OF R E S I S T A N C E ; * CRXF.CRYF • -MAXIMUM DISPLACEMENT OF CR AND THE CORRESPONDING * TIME STEP DIMENSION CMX(1),CMY( 1 ) . A ( 3 . 3 ) , C R X F ( 2 0 ) , C R Y F ( 2 0 ) . C R X M ( 2 0 ) . C C R Y M ( 2 0 ) , C R Y X ( 2 0 ) . C R X Y ( 2 0 ) . I S T E P X ( 2 0 ) . I S T E P Y ( 2 0 ) , C C R X 0 ( 2 0 ) , C R Y 0 ( 2 0 ) IND1=0 IF ( ( A ( 1 , 1 ) . E O . O . ) . O R . ( A ( 2 . 1 ) . E O . O . ) ) GOTO 500 0 R A T I 0 = A B S ( A ( 2 . 1 ) / A ( 1, 1 )) IF ( ( R A T I O . L T . 1 0 0 . ) . A N D . ( R A T I O . G T . 0 . 0 1 ) ) GOTO 5001 IF ( R A T I O . G T . 1 . ) A ( 1 . 1 ) = 0 . 0 IF ( R A T I O . L T . 1 . ) A ( 2 . 1 ) = 0 . 0 GOTO 5 0 0 0 5001 D = S Q R T ( A ( 1 , 1 ) « * 2 + A ( 2 . 1 ) * * 2 ) A N G 1 = A B S ( A T A N ( A ( 2 , 1 ) / A ( 1 , 1 ) ) ) IF ( ( A ( 1 , 1) . G T . O . ) . A N D . ( A ( 2 , 1 ) . G T . 0 . ) ) ANG1=ANG1 IF ( ( A ( 1 , 1 ) . L T . 0 . ) . A N D . ( A ( 2 , 1 ) . G T . 0 . ) ) ANG1=3. 1416-ANG1 IF ( ( A ( 1 . 1 ) . G T . 0 . ) . A N D . ( A ( 2 , 1 ) . L T . 0 . ) ) ANG1=2*3. 1416-ANG1 IF ( ( A ( 1 , 1 ) . L T . O . ) . A N D . ( A ( 2 , 1 ) . L T . 0 . ) \ ANG1=3. 1416 +ANG1 DM=A(3.1)/D DX1=-0M/SIN(ANG1) GOTO 5 0 2 0 5000 CONTINUE IF ( A ( 1 . 1 ) .EO.0.0) GOTO 503 0 D M = A ( 3 . 1 ) / A ( 1 , 1 ) CRY=CMY(IFL)+DM IND1=2 GOTO 5 0 2 0 5 0 3 0 CONTINUE DM = A ( 3 . 1 ) / A ( 2 , 1) CRX = CMX( I F D - D M IN01=1 5 0 2 0 CONTINUE IND2=0 IF ( ( A ( 1 .2) . E O . O . ) . O R . ( A ( 2 . 2 ) . E O . O . ) ) GOTO 6 0 0 0 R A T I 0 = A B S ( A ( 2 . 2 ) / A ( 1 . 2 ) ) IF ( R A T I O . G T . 1 . ) A ( 1 , 2 ) = 0 . 0 IF ( R A T I O . L T . 1 . ) A ( 2 , 2 ) = 0 . 0 GOTO 6 0 0 0 6001 D = S 0 R T ( A ( 1 , 2 ) * * 2 + A ( 2 . 2 ) * * 2 ) A N G 2 = A B S ( A T A N ( A ( 2 , 2 ) / A ( 1 . 2 ) ) ) IF ( ( A ( 1 . 2 ) . G T . O . ) . A N D . ( A ( 2 . 2 ) . G T . O . ) ) ANG2=ANG2 IF ( ( A ( 1 .2) . L T . O . ) . A N D . ( A ( 2 , 2 ) . G T . O . ) ) ANG2=3.1416-ANG2 6000 IF ( ( A ( 1 ,2).GT.O. ) . A N D . ( A ( 2 , 2 ) . L T IF ( ( A ( 1 . 2 ) . L T . O . ) . A N D . ( A ( 2 , 2 ) . L T , DM=A(3.2)/D 0X2=-DM/SIN(ANG2) • GOTO 6 0 2 0 CONTINUE IF ( A ( 1.2).EO.0.0) GOTO 6 0 3 0 D M = A ( 3 . 2 ) / A ( 1 , 2 ) CRY=CMY( I F D + O M ) ) ANG2-2*3.1416-ANG2 ) ) ANG2=3.1416+ANG2 4528 IND2=2 4529 GOTO 60 2 0 4 5 3 0 6 0 3 0 CONTINUE 4531 D M = A ( 3 . 2 ) / A ( 2 . 2 ) 4532 CRX=CMX(IFL)-DM 4533 IND2=1 4534 6 0 2 0 CONTINUE 4535 IF ( ( I N D 1 . E O . O ) . A N D . ( I N D 2 . E 0 . 0 ) ) GOTO 5040 4536 IF ( ( I N D 1 . G T . 0 ) . A N D . ( I N D 2 . G T . 0 ) ) GOTO 5050 4537 IF ( I N D 1 . G T . 0 ) GOTO 5 0 6 0 4538 IF ( I N D 2 . E 0 . 1 ) GOTO 5 0 7 0 4539 D C R X = ( C R Y - C M Y ( I F L ) ) / T A N ( A N G 1 ) 4 5 4 0 CRX=DX1+DCRX + C M X ( I F L ) 4543 GOTO 5080 4544 5 0 7 0 CONTINUE 454 5 DCRY = (CRX-CMX( I F D - D X 1 ) *TAN( ANG1 ) 4546 CRY=CMY(IFL)+DCRY 4548 GOTO 5080 4549 5 0 6 0 CONTINUE 4 5 5 0 IF ( I N D 1 . E Q . 1 ) GOTO 7 0 7 0 4 551 D C R X = ( C R Y - C M Y ( I F L ) ) / T A N ( A N G 2 ) 4552 CRX=DX2+0CRX+CMX(IFL) 4554 GOTO 5080 4555 7 0 7 0 CONTINUE 4556 D C R Y = ( C R X - C M X ( I F L ) - D X 2 ) * T A N ( A N G 2 ) 4557 CRY=CMY(IFL)+OCRY 4 5 5 9 GOTO 5080 4 5 6 0 5 0 4 0 CONTINUE 4561 AA=0X2-DX1 4562 X = A A * T A N ( A N G 2 ) / ( T A N ( A N G 2 ) - T A N ( A N G 1 ) ) 4563 CRX=DX1*X+CMX(IFL) 4564 CRY =CMY(IFL)+X *TAN(ANG1) 4566 5051 GOTO 5080 4567 5 0 5 0 CONTINUE 4569 5 0 8 0 CONTINUE 4569.2 C R X F ( I F L ) = C R X 4569.3 C R Y F ( I F L ) = C R Y 4569.4 IF (I S T E P . N E . 1) GOTO 8002 456 9 . 4 5 C R X 0 ( I F L ) = C R X F ( I F L ) 4569.47 C R Y 0 ( I F L ) = C R Y F ( I F L ) 4 569.5 8 0 0 2 C R X L = C R X - C R X O ( I F L ) 4569.7 ' CRYL=CRY-CRYO ( I F L ) 4 5 7 0 IF ( A B S ( C R X M ( I F L ) ) . G T . A B S f C R X L ) ) GOTO 80 0 0 4571 C R X M ( I F L ) = C R X L 4572 C R Y X ( I F L ) = C R Y L 4573 I S T E P X ( I F L ) = I S T E P 4574 80O0 CONTINUE 4575 IF ( A B S ( C R Y M ( I F L ) ) . G T . A B S ( C R Y L ) ) GOTO 8 0 1 0 4576 C R Y M ( I F L ) = C R Y L 4577 C R X Y ( I F L ) = C R X L 4578 I S T E P Y ( I F L ) = I S T E P 4 5 7 8 . 0 9 8 0 1 0 CONTINUE 4578.18 IF ( I S T E P . L T . N S T E P S ) GOTO 80 2 0 4578.22 C R X M ( I F L ) = C R X M ( I F L ) + C R X O ( I F L ) 4578.24 C R Y X ( I F L ) = C R Y X ( I F L ) + C R Y O ( I F L ) 4 5 7 8 . 2 5 C R X Y ( I F L ) = C R X Y ( I F L ) + C R X O ( I F L ) 4578.26 C R Y M ( I F L ) = C R Y M ( I F L ) + C R Y O ( I F L ) 4578.27 W R I T E ( 1 5 . 8 0 3 0 ) I F L , C R X M ( I F L ) . I S T E P X ( I F L ) L i s t i n g o f D R T B S - 2 ( 4 3 2 0 . 4 5 8 0 ) a t 13:51:52 on APR 6. 1988 f o r CC1d=T0NG o n 4 5 7 8 . 3 6 8 0 3 0 FORMAT('FLOOR' .I 4 , 4X . 'MAX IMUM DISPLACEMENT OF CR IN X-DIR 4 5 7 8 . 3 9 C F 1 2 . 4 , 4 X . ' A T TIME STEP '.16) 4 5 7 8 . 4 5 W R I T E ( 1 5 . 8 0 4 0 ) I F L , C R Y M ( I F L ) . I S T E P Y ( I F L ) 4 5 7 8 . 5 4 8 0 4 0 FORMAT('FLOOR*,I 4,4X. 'MAX IMUM DISPLACEMENT OF CR IN Y-DIR 4 5 7 8 . 5 8 C F 1 2 . 4 . 4 X . ' A T TIME STEP '.16) 4 5 7 8 . 6 3 8 0 2 0 CONTINUE 4 5 7 9 RETURN 4 5 8 0 END ~~