THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING THE ANALYSIS OF STRUCTURAL RESPONSE TO EARTHQUAKE FORCES Glen V. Berg A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the University of Michigan 1958 May, 1958 IP-291

Doctoral Committee: Professor Bruce G. Johnston, Chairman Professor George E. Hay Professor Leo M. Legatski Associate Professor Ernest F. Masur Professor Lawrence C. Maugh Professor Robert M. Thrall

ACKNOWLEDGMENT The author wishes to acknowledge the guidance and assistance received from Professor Bruce G. Johnston, chairman of his doctoral committee, and from the other committee members. Special thanks are due to Professor E. F. Masur, who reviewed the manuscript during its preparation and after completion of the draft and offered numerous valuable suggestions. Professor R. M. Thrall also offered many helpful suggestions and, with his colleagues, undertook the proof of an essential theorem. Professor Thrall generously gave permission to include the proof of the theorem as an appendix to this dissertation. Mr. D. A. DaDeppo evaluated certain matrixes used in the calculations. The Statistical Research Laboratory made available the high-speed digital computer, without which the computations could not have been made. The typing, drafting, and reproduction were done by the staff of the Industry Program of the College of Engineering. To each of the individuals who has contributed to the production of this dissertation, the author is sincerely grateful. ii

TABLE OF CONTENTS Page ACKNOWLEDGMENT................................................. ii LIST OF TABLES................................................. iv LIST OF FIGURES.......................................... v NOMENCLATURE................................................... vii FOREWORD....................................................... x I INTRODUCTION................1................... 1 II ANALYTICAL SOLUTION OF THE DIFFERENTIAL EQUATIONS OF MOTION FOR DAMPED LINEAR DYNAMIC SYSTEMS 5 III EVALUATING THE ANALYTICAL SOLUTION................. 37 IV NUMERICAL SOLUTION OF THE DIFFERENTIAL EQUATIONS OF MOTION FOR NON-LINEAR DYNAMIC SYSTEMS......... 47 V THE EQUATIONS OF MOTION FOR A MULTI-STORY BUILDING SUBJECTED TO EARTHQUAKE...................... 53 VI CALCULATED RESPONSE OF MULTI-STORY BUILDINGS TO EARTHQUAKE GROUND MOTION......................... 71 VII ENERGY RELATIONS................................ 97 VIII SUMMARY AND CONCLUSIONS......................111 APPENDIX A PARTITION THEOREM FOR EUCLIDEAN N-SPACE, by Hans Samelson, R. M. Thrall, and Oscar Wesler.... 115 REFERENCES................................................ 125 iii

LIST OF TABLES Table No. Page I ITERATIVE SOLUTION FOR REAL EIGENVALUE............. 39 II ITERATIVE SOLUTION FOR COMPLEX EIGENVALUE........... 41 III MATRIXES, FOUR-STORY STRUCTURE..................... 77,78 IV MATRIXES, EIGHT-STORY STRUCTURE.................... 79 V MAXIMUM MAGNITUDES OF RESPONSE PARAMETERS, FOURSTORY STRUCTURE, HELENA GROUND MOTION............ 91 VI MAXIMUM MAGNITUDES OF RESPONSE PARAMETERS, FOURSTORY STRUCTURE, FIRST TEN SECONDS OF EL CENTRO GROUND MOTION................................... 95 VII MAXIMUM MAGNITUDES OF RESPONSE PARAMETERS, EIGHTSTORY STRUCTURE, HELENA GROUND MOTION... 96 iv

LIST OF FIGURES Figure Page 1 Damped Linear Dynamic System....................... 6 2 Undamped System................................... 15 3 Damped System.................................. 22 4 Damped System, Damping Linearly Independent of Inertia and Stiffness............................ 30 5 Idealized Moment-Rotation Diagram.................. 56 6 Correction for Plastic Hinge Rotation.............. 58 7 Framing Diagram, Eight-Story Bent................ 72 8 Framing Diagram, Four-Story Bent.7.............. 73 9 Interfloor Viscous Damping........................ 76 10 Index Identification for Matrixes.................. 80 11 East-West Component of Earthquake at Helena, Montana, October 31, 1935....................... 87 12 Response of Four-Story Undamped Elastic Frame, Helena Earthquake................................ 88 13 Response of Four-Story Viscous Damped Frame, Helena Earthquake.............................. 89 14 First Mode Response of Four-Story Viscous Damped Frame, Helena Earthquake......................... 90 15 Response of Four-Story Undamped Frame with Shear Walls, Helena Earthquake......................... 91 16 Response of Four-Story Undamped Elastic Frame, El Centro Earthquake........................... 92 17 Response of Four-Story Elastic-Plastic Frame, El Centro Earthquake............................ 93 18 Energy Per Unit Mass, Four-Story Undamped Elastic Frame, Helena Earthquake................. 106 19 Energy Per Unit Mass, Four-Story Viscous Damped Frame, Helena Earthquake.......................... 107 v

LIST OF FIGURES (Cont'd) Figure Page 20 Energy Per Unit Mass, Four-Story Undamped Frame with Shear Walls, Helena Earthquake........... 108 21 Energy Per Unit Mass, Eight-Story Undamped Elastic Frame, Helena Earthquake...........109.............. log 22 Energy Per Unit Mass, Eight-Story Viscous Damped Frame, Helena Earthquake......................... 110 vi

NOMENCLATURE Symbols are defined where they first appear in the text. Those which appear frequently are listed below for reference. [A] inertia matrix [B damping matrix [c 3 stiffness matrix [D], [E] matrixes in reduced system E strain energy (F) force vector in reduced system G shear rigidity H hinge rotation moment [I] identity matrix M moment P plastic hinge moment Q resistance of frame R resistance function T kinetic energy U shear strength {Y) iteration vector a, ai mass, inertia coefficient b, bij damping coefficient c, cij stiffness coefficient d, dil force-rotation coefficient (f) force vector g generalized force h generalized force, damped system h time interval vii

iV- 1 i,j indices, floor or story k,a indices, moment location 2,m,n,p, q,r constants in Runge-Kutta procedure n index, order of system s index, time step t time w ground displacement x displacement coordinate (x) displacement vector (Y) displacement vector in reduced system Pfi fraction of critical damping 5ij Kronecker delta 71 modal displacement coordinate, damped system X eigenvalue, damped system,L4kJ moment-deflection coefficient Vk2 moment-rotation coefficient story deflection modal displacement coordinate { ~}) modal displacement vector T time parameter of integration plastic hinge rotation {( ) modal column { 4') modal column, damped system undamped natural frequency [~] modal matrix Ah forward difference viii

)( T transpose of a vector [T transpose of a matrix [ i-1 inverse of a matrix Re real part of a complex quantity absolute value determinant Summation convention (Chapters IV to VIII): A repeated index in any term implies summation over the full range of values that the repeated index can assume. ix

FOREWORD In the early morning hours of July 28, 1957, Mexico City was Jolted by the most severe earthquake ever recorded there. More than thirty modern buildings in the heart of the city suffered structural damage. One four-story building and two six-story buildings collapsed completely. Several other tall buildings, under construction at the time of the shock, also collapsed. Tall buildings suffered appreciably more damage than low buildings. Yet the forty-three-story Torre Latino Americana, an office building more than twice as tall as any other building in the city, and located in the area of most severe shock, survived with no structural damage whatever. This was not mere chance, for earthquake forces had been given detailed attention in the design of that building. The catastrophic San Francisco earthquake of 1906 focused a great deal of attention on the problem of minimizing the loss of life and property in earthquakes. Until that time, the notion of designing structures to withstand earthquake shocks had received little attention. True, it was recognized that proper execution of construction details would yield a structure less likely to suffer severe damage in an earthquake; yet the engineer who might have wanted specifically to design a structure adequate to withstand an earthquake without undue sacrifice of economy could have found little to guide him. Soon after the San Francisco earthquake, the Seismological Society of America was founded for the purpose of promoting scientific investigation of earthquake phenomena. Thereafter, steady advances were made toward understanding the origin, propagation, and effects of earthquake disturbances. x

For many years a major obstacle was the lack of accurate information on ground movement in earthquakes. To overcome this, in 1931 the U. S. Coast and Geodetic Survey began a program of measuring and recording ground accelerations during strong-motion earthquakes by the use of automatic recording seismographs. More than sixty seismograph stations are now in operation in the western United States. More recently, in 1950, the Earthquake Engineering Research Institute was founded to encourage advancement in the scientific design of structures in earthquake areas. The United States has not been alone in its search for answers to earthquake problems. Many other countries throughout the world have regions of seismic activity, and most of these countries are active in seismic research. Notably, Japan has contributed a great deal to the advancement of seismology and structural engineering through its Earthquake Research Institute. Earthquake problems related to engineering might be divided into three principal fields: First, there are the problems of ground motion. This field covers measuring, recording and analyzing ground motion in earthquakes, studying the origin of earthquake shocks and their propagation through the earth, and the yet unsolved problem of predicting the locations and intensities of future earthquakes. Second, there are the problems of structural mechanics. In this field, the primary concern is the response of structures to ground motion, methods of calculating response, and finding the effects of structural properties on the response. Third, there is the problem of taking the findings in the first two fields, along with a knowledge of engineering materials, observations of earthquake damage, and a generous amount of practical experience, and formuxi

lating a rational set of design rules to guide the practicing engineer in the design of earthquake-resistant structures. This dissertation is intended primarily as a contribution to the second of these fields,that of structural mechanics, with a brief look at the relation of this field to current design practice. xii

CHAPTER I INTRODUCTION The basic tool of the theory of linear vibrations was forged by Daniel Bernoulli in 1753, when he enunciated the principle of resolution of vibrations into independent modes. The general theory of vibration of an undamped linear dynamic system with a finite number of degrees of freedom was developed by Lagrange in 1762-1765. The work of these two great mathematicians has served as a basis for most of the developments in the theory of vibrations over the past two centuries. Today, the formal solution of the problem of linear vibrations, both damped and undamped, is complete and has been extended to systems with infinitely many degrees of freedom(l'2). The theory of vibrations of non-linear systems has not fared as well as its linear counterpart. Bernoulli's principle of resolution of vibrations into independent modes, which plays such an important role in linear theory, is inapplicable to non-linear systems. Exact solutions have been obtained only for a few special cases of single degree of freedom systems(3'4). For steady-state forced vibrations of single degree of freedom systems with non-linear damping, methods have been devised for finding approximately equivalent linear systems(5). For more general vibrations it may sometimes be possible to approximate the system by a lines system over a limited range, or by one of the non-linear systems for which the solution is known. If strong non-linearities are present, this approach is usually unsatisfactory; in such cases, for single degree of freedom systems, the solution can often be obtained by the use of an electronic analog computer(6) Alternatively, numerical methods can be employed(6'7'8 -1

-2For a non-linear system having many degrees of freedom, the complex circuitry required for the analog computer becomes prohibitive, and numerical methods provide the only suitable approach(9). Except for the simplest cases, a high-speed digital computer is essential. In nearly all earthquake response analyses accomplished to date, the procedure has been to represent the structure by a linear. dynamic system. An undamped system is usually used to transform the differential equations of motion into uncoupled form, and linear damping is then introduced into the uncoupled equations as a fraction of critical damping in each mode(6'10). A method has been developed(ll) for applying this procedure to multi-story elastic-plastic shear buildings, using a technique mentioned previously; namely, replacing the non-linear elasticplastic system by a linear elastic system over a limited range. To this writing, so far as the author has been able to determine, no actual calculation of the response of a multi-story elastic-plastic structure to earthquake ground motion has been reported. This dissertation presents, first, the classical method of normal modes for solving the differential equations of motion of undamped linear dynamic systems having finitely many degrees of freedom. Conditions are then established for which a damped linear system can be solved by first obtaining the uncoupled equations of motion for the corresponding undamped system, and then modifying the uncoupled equations to account for damping. This is valid only in special cases. For the general case, a method is presented for deriving uncoupled equations of motion for the damped linear system. This method is mathematically exact, and does not rely upon "small damping" for its validity.

-3Second, methods are given for evaluating the analytical solutions. The methods employ iteration and numerical integration, and are adaptable to the high-speed digital computer. Third, a numerical method is presented for solving the differential equations of motion for non-linear systems. This is a general method, and imposes no restrictions on the non-linearities other than piecewise continuity of the functions. Fourth, the differential equations of motion are set up for the response of a multi-story building to earthquake ground motion, taking into account elastic-plastic deformation in the frame members, a modified Coulomb damping due to inelastic deformation of shear walls and partitions, and viscous damping. Uniqueness of the solution for an elasticplastic frame is proved. Fifth, a number of solutions are presented showing the computed response of four- and eight-story buildings to recorded earthquake ground motions. Various systems, both linear and non-linear, are used to represent the buildings, illustrating qualitatively the effects of variations in the systems. Although it is not intended that the results be interpreted quantitatively, a comparison is made between computed results and building code provisions. Finally, energy equations are given and energy curves are plotted for several of the response solutions. A comparison of the energy and response curves shows clearly the important role of dissipative forces in reducing the amplitude of vibration.

CHAPTER II ANALYTICAL SOLUTION OF THE DIFFERENTIAL EQUATIONS OF MOTION FOR DAMPED LINEAR DYNAMIC SYSTEMS To write the differential equations of motion for a damped linear dynamic system, consider a system of n discrete bodies, each of which is connected to every other body and to the base by a system of linear springs and viscous dampers. Such a system, of order two, is shown symbolically in Figure 1. Assume that the springs and dampers are weightless, and that all motion takes place in the x-direction, which is the direction in which the springs and dampers act. Let the motion of the base in the x-direction be any prescribed function of time; and, further, let each body be acted upon by an external force in the x-direction, each force being any prescribed function of time. We then have a linear dynamic system of n degrees of freedom. Let ai = mass of the ith body, (FL-1T2)* Xi = displacement of the ith body, relative to a coordinate system fixed in the base, (L) w(t) = displacement of the base, relative to a "fixed" frame of reference, (L) Pi(t) = external force applied to the ith body, (F) and let dots denote differentiation with respect to time. Define the stiffness coefficient cij as the force exerted on the ith body by the springs when the configuration of the system is * Dimensions are shown in units of Force, Length, and Time. Other units are admissible, as long as consistency is maintained. -5

=w(t) X xI — x2 xl and x2 are displacements relative to base w is displacement of base relative to "fixed" reference axis. Figure 1. Damped Linear Dynamic System

-7xj = 1 unit length; xk = O, k A j; and xk = 0 for all k. The dimension of cij is FL-1. Define the damping coefficient bij as the force exerted on the ith body by the dampers when the configuration of the system is x = 1 unit velocity; xk = 0, k ~ j; and xk 0 O for all k. The dimension of bij is FL-1T. Because all the bodies are interconnected by the spring and damper system, a displacement or velocity in any one body will cause spring or damper forces to be exerted on all bodies in the system. For any configuration, the force exerted on the ith body by the spring and damper system is The Newtonian equation of motion for the ith body is then n n ~t( ~t@ iij C ~ ~= P,(t) (2.1) Thus, for the entire system, Xta, ~ ~ ~ J jx(t) QiZ(t). i,.n) (2.2) Let fi(t) P(t)- a= Pi(t). n) (23) Equations (2.2) then become n Q~LC, j:, C~Xi fj (t — (L\ ) (2.4) C) C)ZbJX jz J $

-8which are the general differential equations of motion for a damped linear dynamic system. In the earthquake response problem, Pi(t) = 0, and Equations (2.2) become a special case of Equations (2.4). To preserve generality, the general form will be used in the derivations which follow. A system of n homogeneous second order differential equations with constant coefficients can always be transformed into a single differential equation of order 2n, or, alternatively, can be reduced to a system of 2n differential equations of first order(l2) Also, by means of suitable change of variable, a system of n simultaneous homogeneous differential equations can be transformed into a system of n equations each containing only one variable. This is an application of the Smith canonical form for polynomial matrixes(l3). The latter technique will be employed in this chapter, applied to the system of second order equations for the undamped case, and to the reduced system of first order equations for the damped case. To facilitate the mathematical derivations, matrix notation will be employed(l). Let the inertia matrix be a1 O.. 0 0 a2. 0 [A]..... 0 0.. an

-9let the damping matrix be bll bl2 ~. bin b21 b22 ~ * b2n [B]=. bnl bn2. bnn and let the stiffness matrix be C11 C12 * Cln C21 C22 ~ * C2n [Ic] Cnl Cn2 * * Cnn Let the displacement vector be (x) =nJ and the driving force vector be fl(t) f2(t) (f(t)) =. fn(t)

-10 - Equations (2.4) can then be written [A](X) + [ e{k) + [c]{X) = (f(t)). (2.5) We shall impose the conditions that [A] be a diagonal matrix with positive diagonal elements; CB] and [C] be symmetric; and [A], [B], [C], and {f(t)} be real. These conditions are always satisfied for a structural system. We shall consider first the classical solution for the undamped system. The equation for the undamped system is [Al{X} + [C]{)x) = (ft). (2.6) To solve this equation, one must first solve the homogeneous equation [Al()' [C](x) 0. (2.7) The solution takes the form {(x ) ke iw {Vl (2.8) in which k is a scalar of dimension L,w is a scalar of dimension T-1 and ({) is a dimensionless vector. Then {(ta =-) -k ev t( l}: (2.9) and Equation (2.7) becomes _ k e; wt 4(4 [A]{- [C[]{#)} = Q (2.10) [ILr]- k [A] [C]]{[) - 0, (2.11) in which [I] is the identity matrix. The conditions imposed on [A] guarantee the existence of its inverse [A]-1. Equation (2.11) can have a non-trivial solution only if 1'[1- - [A]I [C]I 0. (2.12) Equation (2.12) is called the characteristic equation of [[A]-l[C]] and its roots 3 are called the latent roots or eigenvalues

-11 - of [[A]l[CI]. In any stable system, [A] and [C] are both positive definite; that is to say, they are symmetric and their principal minors are all positive. These conditions guarantee that the eigenvalues of [[A]-l[C]] are real and positive. For each eigenvalue LW, there exists a non-zero eigenvector {(j) which satisfies Equation (2.11). Since Equation (2.11) is homogeneous, any scalar multiple of {iJ} will also satisfy the equation. In view of this, it is convenient to assign unit value to the last element of ()} and determine the rest of the elements accordingly. The eigenvector (}j) is also called the modal column associated with a?. To establish the orthogonality of the modal columns, we have, from Equation (2.10), X [A~](j}- [[C](~J) 0 (2.13) and w (ck})T[A] _ - (,kT[C1]T.o (2.14) Premultiplying Equation (2.13) by {pk}T and postmultiplying Equation (2.14) by ({j}, one gets C +) [Al{i - 9k}T[C]({ 0i) 0 (2.15) and {k)Yi [Al]7} (_ 7 [C] ( } = o 0 (2.16) By hypothesis, [A] - [A]T and [C] = [C]T. Therefore, (Cj CO) (OkT[A](5i) o (2.17) from which ({6) [A] (J) = 0. (j # k)* (2.18) Equation (2.15) then gives { ) C]{<j 1 0. (C j k)* (2.19) Equations (2.18) and (2.19) are the orthogonality properties of the modal columns for an undamped system. * The case of repeated eigenvalues is excluded.

-12If the w's are all distinct, the n modal columns {fj} form a linearly independent set. One can then represent any displacement vector{ x} as a linear combination of the ()'ts; that is, one can write (x) = Z(sXYj), (2.20) j=l in which tj is a scalar of dimension L. Let us define the modal matrix [] as the matrix which has for its columns the n modal columns {~), 2~}. (on). Also, let us define the modal displacement vector { } as the vector comprising the elements tj, S2,...'n. Equation (2.20) can then be written in the form (> = [d(' } (2.21) or, equivalently, {t) = CH] {X), (2.22) The columns of [i] are linearly independent; therefore, [i] is nonsingular, and [~]-l exists. Consider the case of free vibration of an undamped system, with (x) and (k) prescribed at t = O. Combining Equations (2.7) and (2.20), one obtains nIn Z[A]{(#) t) + L[C]{(i)(tg) 0 O. (2.23) j=l J j:l Premultiplying each term by (Ok)T gives n n kif{)'A]Ej [A f [C]{~j)A rt = ~- (2.2o4) j=' j Because of the orthogonality properties of the modal columns, all terms vanish except the ones for which j = k, leaving (ek)T[A]() k} + (I k)rC](tk)ik = (k- 2z n) (2.25)

-13From Equation (2.10), [c] {) = C [A ( k), (2.26) which reduces Equation (2.25) to (ok}T[A]k} 4 0,k T (2.27) ( k)T[AW ]{pk) ~k + Ok k> [A] }k 0 or simply.. (k. ) (2.28) Sk - k Sk. (k The solution is (2.29) ik(t)= 6k~) 5 in Okt + Sk(O) COS Wkt Note that reversing the sign of a does not change Sk(t). Whether we use the positive or negative value of co is immaterial. Combining this with Equation (2.20), the complete solution for free vibration with prescribed initial conditions is (x(t)) =(t(n}( si-n et+ f(.) COS't) (2.30) From Equation (2.30) it can be observed that if some Sj(0) = gj(0) = 0, the mode corresponding with that w) does not appear in the motion. Specifically, if gj(O) = [j(O) = 0 for j h p, and p(0) and =p(O) are not both zero, (x(t)) will be proportional to {I}? at all times. Consider now the case of forced vibration of an undamped system. Again one can use the modal series expansion of (x) given in Equation (2.20). Combining this with Equation (2.6) yields j-I \ j=l:{ *[A <} 4+ =_([]t)}. (2.31) Premultiplying each term by {,k}T one gets At~?(~ ) A { ]) {~, k}T(f(t)) (2.32) j i jzI

-14In view of the orthogonality properties, this reduces to (k}T1[A]{} Vi) + (k)T[C{l k}) = (k)T f(t)). (2.33) With the use of Equation (2.26) we obtain (k) [A]{ k)k + O(k L[A]I(k)k = (k){f(t)) (2.34) or simply (2.35) 4- W kk(t) in which gk ( _) = k}T~ r(f (t)} (2.36) The function gk(t) is called the generalized force. It has dimension LT-2 The solution of Equation (2.35), for zero initial conditions, is given by Duhamel's integral -k g() s nlk (t S) d r (2.37) Here again it makes no difference whether we choose the positive or negative value of Wck. Example —Consider the undamped elastic frame shown in Figure 2, which is equivalent to the spring and mass system shown. For either system, [A] kip sec2 in-1,

-15 - I= 19,160 in4 772 kips o I -c rO Cmm'I cu 4 I=7500 in4 1158 kips ) 0M o,, cur L: 30' = E= 30000 kip in2 Elastic Frame -40 kip in 260nkip in 240 kip'in Equivalent Spring and Mass System Figure 2. lUndamped System

and F500 -2401 [c] g = Ikip in-. -240 200j Equation (2.12) then becomes (2 _ 500) 80 =0, (2.38) 120 (M2 - 100) which expands into the polynomial co)4 - 830 2 21, == 0. (2.39) The roots are 2o = 29.8388 sec-2 and (2.40) 2 = 236.8278 sec2. Equation (2.11) is ('2 500) 80 120 (W2 _ 100)1 (, (2.41) from which 8- 8o0 j 1 ( 100 {02 500 -120 in which the unit value for the second element is arbitrarily assigned. Evaluating Equation (2.42) for the two values of a2 gives ({L}) ={5847} and (2. 43 ) 12} {-.1402}

-17or.] 5847 -1.14o2 (2. [+]= 01 1 j (2.44) The inverse of the modal matrix is _-I *.57974.66104 []= [.57974.33896 Let the system be initially displaced to the position (X(0)} = 1 in., and let it be released from that position with zero velocity at t = 0. From Equations (2.22) and (2.45), the modal displacement and velocity vectors at t = 0 are ~O~o~= [Q~1 r 1.24081 {(0)} K x(O) = in. (2.46) -.2408 and {(o = [j]-l {x(o)} = (2.47) Putting the above results into Equation (2.30), one obtains.for the motion of the system.7255.2745) {(t)} -. } cos 5.4625t + J4cot {5.3892t in. (2.48) tl.24o8- l-2408J Let the undamped system be subjected to a constant driving force equal to half its weight. That is, {f(t)} = 38J kips, t> 0 Also, let the system be at rest in the position of equilibrium at t = 0.

Then, from Equation (2.36), {.5847 1} f579 gl(t) 5847 1 [3 ] 58 47) = +239.47 in.sec-2 (2.49) 0 2 1 and (-1.1402 1} 579 g2(t) = -1.1402 1} [3 3 ] -1.1402= -46.47 in.sec2 (2.50) Equation (2.37) then gives (t) = 5.4625 239.47 sin 5.4625 (t - T) dr (2.51) = 8.O254 (1 - cos 5.4625t) in. and t E2(t) = 15f3892 I -46.47 sin 15.3892 (t - T) d(2.52) = -0.1962 (1 - cos 15.3892t) in. Combining these results with EQuations (2.43) and (2.20), one gets r4.6923 (1 - cos 5.4625t) +.2237 (1 - cos 15.3892t)) } 8.0254 (1 - cos 5.4625t) -.1962 (1 - cos 15.3892t) (4.9160o F-4.6923 F-.2237 = j 9 + Jcos 5.4625t + P cos 15.3892t in. (2.53) 17.8292 L-8.0254 +.1962J In analyzing the motion of a damped linear system, an expedient frequently employed is to drop the damping term temporarily and find the eigenvalues, modal columns, and uncoupled equations of motion (Equations (2.35) ) for the undamped system. A damping coefficient, expressed as a "fraction of critical damping," is then assigned for each mnode and inserted into Equations (2.35), converting them to ik~< ~ 2 Z _ (t), (k=, 2,.. n) (2.54)

in which Pk is the fraction of critical damping for the kth mode. Equations (2.54) are then solved for k and the displacements are determined from Equation (2.20). For the case of forced vibration with zero initial conditions, the solution of Equation (2.54) is _k = ___ P _p'K_(t-__) _ (t-C) dT. (2.55) For free vibration, the solution is = IkLinkk -t5 t( (2.56) + () coS k - t In using this procedure, one implicitly assumes that the modal columns for the undamped system remain valid for the damped system as well. In general, this is not so. However, it turns out that for the very special case in which the damping matrix is a linear combination of the inertia and stiffness matrixes, the undamped modal columns do remain valid and lead to completely uncoupled equations. For this case, let [B] = o[A] + [c], (2.57) where a is a constant of dimension T and 7 is a constant of dimension T. Equation (2.5) can then be written [AJ(x) + o [A](x) + 7 [C](x) + [C](x) = (f(t)) (2.58) It is claimed that the undamped modal columns are valid for this case.. If this is so, Equation (2.20) is also valid, and we can combine this

-20with Equation (2.58) to get n n n j_ \ j.i Jt (2.59) n _C_': ( f (t) Premultiplying each term by (k)T and taking into account the orthogonality properties (Equations (2.18) and (2.19) ), one gets { kT([A]{I:) (', I ~) + {ok)C]{ ~ ({.( tk) (2.60) But, by virtue of Equation (2.13), [C]{k) = 2o[A]{ k), (2.61) which reduces Equation (2.60) to - (' k) k + 0okZ k gk(t), (k=,2.) (2.62) in which gk(t) is defined by Equation (2.36). If we let 4 t +0k ~k k (2.63) 2 CC~k Equation (2.62) becomes identical with Equation (2.54). Not only are the equations completely uncoupled, but they are in terms of the undamped frequencies as well. Equations (2.55) and (2.56) give the solutions for forced and free vibration. One should note particularly from Equation (2.63) that if [B] is proportional to [C], (a = 0), Pk is directly proportional to wC so that the fraction of critical damping is higher in the modes having higher undamped frequencies. On the other hand, if [1B is proportional to LAl,

-21(y = o), Pk is inversely proportional to wk so that the fraction of critical damping is lower in the modes having higher undamped frequencies. For example, let the undamped system of Figure 2 be modified by the addition of damping, as shown in Figure 3. The damping matrix is 14 -6 [B] = ["6 1kip sec in-l, (2.64) which corresponds to Equation (2.57y) with z = 0.5 sec-l (2.65) 7 = 0.025 sec. The values of cu and [L] were found- in the previous example. They are 1 = 29.8388 sec-2, (2.40) 1 = 236.8278 sec-2.5847 -1.1402 ["] = (2.44) Evaluating Equation (2.63), one gets P1 =.1140, (2.66) ~2 =.2086. Let the system be displaced to ({x in., and let it be released with zero velocity at t = 0. The values of {((0)} and {*(0)} are exactly the same as for the previous example, namely {g(O)} = 4 in.,'(2.46) - 2408(2.47) {~}={}. (2.47)

-22I = 19160 in4 772 kips to o 6 kip sec in Eqiv I=e7500 in4 F 1158 kips t/ o 8 kip sec. in E=30000 kip in-:30' Damped Elastic Frame -40 kip in260 kip in 1 240 kip in' -I 8 kip sec in 6 kip sec in Equivalent Damped Spring and Mass System Figure 3. Damped System

-23 - Using these numerical results to evaluate Equation (2.56), one obtains 1 = e-'6230t(.1424 sin 5.4268t + 1.2408 cos 5.4268t) (2.67) 2 = e-32103t(-.0514 sin 15.0506t -.2408 cos 15.0506t) (2.68) The above results, with Equation (2.21), lead to'7302t 6230t x) 7302 e 6230t cos(5.4268t -.1143) 1.2489j + e'7 -3.2103t cos(15.0506t -.2102) in. (2.69) -.2462 Consider now a system in which the matrix [B] may be linearly independent of [A] and [C]. What follows is equally applicable if [B] is not linearly independent, and in that special case it leads to precisely the same result as the method above, although by a more laborious route. The equation of motion [A]"{) + [B]j(} + [C](x} = {f(t) (2.5) can also be written in reduced form(l,2) [D]() t [E](y) = (F) (2.70) in which [D ] [] (2.71) [A] [0] (2.72) [ x0 (2.73)[o ]c X~,

-24and F{ ) {(f(t)}} (2.74) Note that [D], [E] and (y) are not homogeneous dimensionally. To solve Equation (2.70), one must first solve the homogeneous system [D]({} + [E](y) = 0, (2.75) for which the solution takes the form ({y) = ke t(y}), (2.76) Putting this and its derivative into Equation (2.75), one gets k eXt{M[D( E) + (2.77) or [xLII1 + [D]'[E]1(,} - (2.78) That [D]-1 exists can be seen from the relation [D-I ([A]-[7B][Al [A] (2.79) an inverse; hence [D] also has an inverse. The values of X for which Equation (2.78) has non-trivial solutions are the eigenvalues of [-[D]- [EJ]. They occur as complex conjugate pairs*, and for a stable system they have negative real parts. The eigenvectors (y) also occur as complex conjugates; and since they * FDr the undamped system, we were c9ncerned with the eigenvalues of LAj-lc Jj. ~Because both [A] and IC] are positive definite, the eigenyaiues were real and positive. In the present case, neither [D] nor LEjis positive definite; hence the eigenvalues may be complex.

-25can be determined only to within an arbitrary scalar multiplier, it will be convenient to assign unit value to the last element of each eigenvector and adjust the remaining elements accordingly. Define a modal column {J) associated with the eigenvalue Xj such that the solution of Equation (2.75), expressed in terms of the original coordinates, is of the form (X3= k xJi t(j) (2.80) Differentiating this and putting the results into Equations (2.73) and (2.76) gives iy)= kei = e Njtf (2.81) or ~or r{ (A }} (2.82)The modal columns ({) occur as complex conjugate pairs. They are not the same as the modal columns for the undamped system. To establish orthogonality of the eigenvectors, we have, from Equation (2.77), [DlI } + [E]) - (2.83) and Xk( [D] {I}[T [ 0] E (2.84) Premultiplying Equation (2.83) by~ek}T and postmultiplying Equation (2.84) by (&}, one obtains hAyi{8T[D]{8i) t (i }T[E1i } = 0 (2.85)

-26 - and X (kki}[D]TY} + (k}T[E]lkh } O. (2.86) By hypothesis, [A], [B] and [C] are symmetric. From Equations (2.71) and (2.72) it follows that [D] and [E] are symmetric; that is, [D] = [D]T and [El = [E]T. Therefore, subtracting Equation (2.86) from (2.85) gives (X; -kh)(\i [lD]I} OQ (2.87) from which ({k) L[DiC]{\ 0 (j k)* (2.88) and this with Equation (2.85) gives ik)T[E]tIJ= O. (j k) (2.89) Equations (2.88) and (2.89) are the orthogonality properties of the eigenvectors. Expanding Equation (2.88) with the aid of Equations (2.71) and (2.82), one gets [X{pk} ] [o] [Aj x J}1 { }J i[A] [ { (} (2.90) which yields ( A kX){'k}T[Al{'} +t ai)['{} ~ ( # )(2.91) Similarly, Equations (2.72), (2.82), and (2.89) lead to FA(S(;7}1T l[A]1 [o] xs{j} }| ~(2.92) kj [0 [c] {} xJ * The case of repeated eigenvalues is excluded.

-27or T- >\{k}T[Ai{J} (#kf[C1{"j) 0. (j ) (2.93) -- x~ ~{+~}~+[Al{J} [C](,J} o. 6 k) If LBI is a linear combination of [A] and [C], the Orthogonality Equations (2.91) and (2.92) are equivalent to Orthogonality Equations (2.18) and (2.19). If the 2n eigenvalues are all distinct, the 2n eigenvectors form a linearly independent set, and any configuration of the system can be expressed as a linear combination of the eigenvectors. That is, we can write Zn (,y~t)} ~ g (J} 5 (t) (2.94) where each nj is a scalar function which may be complex. Equation (2..75) then becomes Zn Zn [o]{(,}5 j [E]{ i} j - {F} (2.95) Premultiplying by {k}T gives 2\v~.? ]{8j..j + yY} []j{y}j {4j F} (2.96) j-I In view of the orthogonality properties, all terms on the left side vanish except those for which j = k,* giving (k}T[oD] {k} + (} -[E]{k} S (}(F (2.97) From Equation (2.83), [E j] |} =- k [] (2.98) * Again we exclude the case of repeated eigenvalues.

-28which reduces Equation (2.97) to S - kk -- h (t) (k- t,. n) (2.99) where (Yw Y{ F} (,,' (f (t)) hk(t) -y{})r[]{tkt ) Ak(( )T[A].io + (}'[f] ) (2.100) The solution of Equation (2.99) is Ckk ~t; ko h(2.101) rZ k ( t ): Yq(0)ct + e k(t- )h ( — d'r, which, with Equation (2.94), describes completely the motion of the system. Since the X's, (I)'s, h's, and r's may all be complex, it might seem that iy(t)} could also be complex. However, it turns out that (y(t) remains real at all times, as will be shown. First of all, [A], [B] and [C] are all real, and it follows from Equations (2.72) and (2.79) that -[D1-1 [E] is real. Therefore, if there is an eigenvalue of -[D]-1 [El which is not real, its conjugate is also an eigenvalue. Let Xp and Xp+l be a conjugate pair of eigenvalues. From Equation (2.83) [Xp[D + [E]]{Y}= 0, (2.102) and [p+lILD] E]]IyP Y* 0. (2.103) Replacing each element in Equation (2.102) by its conjugate gives [ [x ]~1+ [E]] (fP} = o. (2.104)

-29But since [DI and [E] are real, and X = Xp+1, this becomes opts[D]+ t[E]] (} =o. (2.105) In view of Equations (2.103) and (2.105), {1p+l} must be a scalar multiple of {P}; and since a unit value has been assigned for the last element of each, {y'I}= {P } (2.106) Also, since [F] and [D] are real, Equations (2.100) and (2.106) imply that hp+,= -, (2.107) and it then follows from Equation (2.99) that p+l = -p (2.108) In view of the above, the terms of {y(t)} in Equation (2.94) corresponding to the conjugate pair of roots p and Xp+l are {\GAP}y + ~Pl pr( = {Y} + {i}p = 2 Pe({(Y P}) (2.109) For example, consider the system shown in Figure 4, which is the same as that used in the earlier examples except for the damping elements. The matrixes are [A]= 3 j kip sec2 in-1 14 -10 [B] = kip sec in-1 -10 10 [a~] = - kip in-1 -240 200

-30I- 19160 in4' 772 kips n qDm -I Fa 10 kip sec in H /= 7500 in 1 58 kips._ I1'.o p 4 kip sec inps n 260 kiEp in 240 Dip inS Vat X40 dip sec ~i 4kip sec in 0 kipin sc Equivalent Damped Spring and Mass System Figure 4. Damped System, Damping Linearly Independent of Inertia and Stiffness

-31From Equations (2.72) and (2.79), [D] [A ]A-'C]] -[I-] [] (2.110) Evaluating this for the present case, and putting the result into Equation (2.78), one gets 14 10 500 (x+y) - 3 -5 (x+5) -120 100 -1 0 X 0 I{Yj} = o (2.111) 0 -1 0 X The eigenvalues are the values of X for which the determinant of the matrix is zero. Expanding the determinant and setting it equal to zero gives the polynomial 4 + 29 t3 + 820 2 21200.11.3 3 3 The roots are X1 = -0.5106 + 5.4675 i ~2 = -0.5106 - 5.4675 i (2.113) X3 = -4.3228 + 14.6855 i X4 = -4.3228 - 14.6855 i From Equations (2.82) and (2.111), one can solve for the modal columns in terms of X, obtaining 13 +80 +5A 100..~}3.5 = 2<l} 14 - 500 5 X> + IZ 0 (2.114)

-32Here again, the unit value for the last element of the modal column is arbitrarily assigned. Putting the above values of X into Equation (2.114) yields { _ { 6.5886 +.0L4822i1 and'p (2.115) } -10350 +.2283i The complete eigenvectors are -.5642 + 3.1937i lJoyl z1 =3Y -.5106 + 5.4675i l I ~l -.5886 +.04822i and (2.116) 1.1215 - 16.1863i fJ+3l _ JMAl _ J-4.3228 + 14.6855i -1.0350 +.2283i _ 1 It may be well to point out that finding the eigenvalues and eigenvectors by constructing the characteristic equation of the matrix and and finding its roots, as was done above, is entirely impractical except for systems of order three or less. For higher order systems it is easier to find the eigenvalues and eigenvectors by matrix iteration. In fact, a method sometimes used to solve a high order polynomial is just the reverse of the above, namely, to construct a matrix having the polynomial for its characteristic equation and then find the eigenvalues by matrix iteration(14). A suitable iterative procedure for finding the dominant eigenvalue and eigenvector is presented in Chapter III.

-33 - Let the system of Figure 4 be displaced to {x) = { and released in free vibration with zero initial velocity. Then, from Equation (2.109), o-.5642 + 3.1937i 1.1215- 16.1863i (O 0 2 Re -.5106 + 5.4675i -4.3228 + 14.6855i {y(0) |.5886 +.04822i} r(0) + -1.0350 +.2283ij 3(0)J Ll1 1 (2.117) which leads to rql(O) =.6295 -.02576i (2.118)'13(0) = -.1295 +.02584i Using Equation (2.101), with hk = 0 for free vibration, one obtains rll(t) = (.6295 -.02576i) e(- 5106 + 5.4675i)t 713(t) = (-.1295 +.02584i) e(-4.3228 + 14.6855i)t 9) and Equations (2.73), (2.94), and (2.109) lead to.5886 +.048221 (-1.0350 +.2283i) {x(t)} =2 Re 1 rl(t) + 2 Re 13(t (2.120) Combining Equations (2.119) and (2.120) and simplifying gives -. 5106t {.7442 cos(5.4675t +.0409) =x3 e | 1.2602 cos(5.4675t -.0409) + e 3228t {.2800 cos(14.6855t -.4138) in. (2.121) -.2642 cos(14.6855t -.1969) ( This can be compared with Equation (2.69) to see the effects of the change in damping. For an example of forced vibration, let the same system be at rest in the equilibrium position at t = 0, and let it then be subjected to a constant driving force equal to half its weight; that is,

-34From Equation (2.100),.586 +.Q48E2j2 i2 1 579} 2(-,5lo5.,4675i)L56+'o4A~i I~ L1 aJ 1f + (2.122),588d36;O482 1Z 114 -lol { 56&66t.O4Z2i1 =-.4389 - 2 I?2 3 L and {-O5 0Q~5 1 t.2283 L 1} 79 2(-4 32?ki~4M4855L) -lI O35o+263i I[3 o])-1v.0350 t2283i + {-t0550~Z2a3i 1}(4 -O] ({.0350+. 2283 ( =.4 389 t+.5512 Inserting these values in Equation (2.101) and carrying out the integration, one gets l = (-4.0037 +.4541i) (e(''5106 + 5.4675i)t - 1 and (2.124) T3 = (.08911 -.05611i) (e(-4'3228 + 14.6855i)t _ 1i and finally this, with Equations (2.94) and (2.109), leads to 4(t).9160 + e-'5106t f-4.7595 cos(5.4675t -.0312) xt 7.82921) t -8.0589 cos(5.4675t -.1129) + e-4 3228t -.2232 cos(14.6855t -.7791) + e -'l+.2106 cos(l4.6855t -.5619)f i. (2.125) The following observations are worthy of note. If the damping matrix is linearly dependent on the inertia and stiffness matrixes, the

-35 - modal columns are real; in any one mode of vibration, all bodies in the system oscillate in phase, and the natural frequency of the damped system in any mode is lower than the undamped natural frequency in that mode. Also, in any one mode of free vibration, the system passes through the equilibrium configuration twice in each cycle of oscillation. On the other hand, if the damping matrix is linearly independent of the inertia and stiffness matrixes, the modal columns are in general complex; in any one mode of vibration, all bodies in the system oscillate at the same natural frequency, but not necessarily in phase; and the natural frequency of the damped system in any one mode is not necessarily lower than the undamped natural frequency in that mode. In the last example, the damped frequency in the fundamental mode is in fact higher than the undamped frequency in the fundamental mode. Also, the system will not necessarily pass through the equilibrium configuration in its oscillations in a single mode of free vibration. For a stable damped system the X's are complex, with negative real part. The natural frequency of the system in any mode, in radians per second, is given by the imaginary part of X, and the exponential term is indicated by the real part of X. In other words, if X = L+io, (2.126) the elements of the displacement vector in the corresponding mode of free vibration will be of the form keLt cos(pt - o), (2.127) in which k and 0 are real constants. This in itself tells a great deal about the vibrational characteristics of the system.

CHAPTER III EVALUATING THE ANALYTICAL SOLUTION In Chapter II it was shown that the motion of a linear dynamic system can be resolved into modes of vibration, and that the motion in each mode can be expressed as the product of a constant vector, the modal column, and a scalar variable in the form of Duhamel's integral. The eigenvalues and generalized forces which appear in Duhamel's integral, and the modal columns as well, may be complex. In carrying out the analytical solution, the major problems are to find the eigenvalues and modal columns, and, having these, to evaluate Duhamel's integral. One iterative method for finding the dominant eigenvalue of a matrix arises from the properties of high powers of a matrix(l). Let [U] be a real matrix having a dominant eigenvalue v, and a corresponding eigenvector (I). The "dominant" eigenvalue is the eigenvalue of greatest modulus, whether real or complex. Let (YO) be an arbitrary non-zero real vector. Denote {Y1} = [U]{YO, (31) and in general {Ym+l) = [U]{Ym). (3.2) Let V0 be any non-vanishing element of {Yo), and let Vm be the corresponding element of Ym) If the dominant eigenvalue v is real and distinct, it is known that Vm - as — m c (3.3) Vm -3 -37

-38and an VY} a s n-. (3.4) /m For example, consider the undamped system of Figure 2. The equation of interest is [o [I] - [Al'[cj]{ = o. (2.11) To express this in terms of the above notation, let [U] = [A]-l[c] v 2 (3.5) and ({) = ()} Then [vI]- Lu] ]{)=o. (3.6) Inserting numerical values in Equation (3.2), one gets for the iteration 500 -80 { m)= {Ym-1) (3.7) -120 100 If one starts with and takes Vm to be the last element of {m the iteration gives the results results shown in Table I.

TABLE I. ITERATIVE SOLUTION- FOR REAL EIGENVALUE m 1 2 3 4 5 6 -8o -2.1333 x 104 -5.1236 x lo6 -1.2155 x 109 -2.8793 x 1011 -6.8193 x 1013 YMf 100 1.96 x 1 4.5200 x 106 i.o668 x lo9 2.5254 x 1011 5.9807 x 1013 m- 100 196 230.61 236.02 236.72 236.82 rm-i L'J

-40From the last column of Table I, V6 V = V = 02 = 236.82, (3.8) and {Y6},,.>={Y)=t<}= 0 - ) ~ -(3.9) which compare well with the results found in Equations (2.40) and (2.43). Usually the mode having the lowest frequency is the most important mode in the response. This can be found by the same iterative scheme after rewriting Equation (2.11) in the form [[c]'[A] -' ]i= ) o (3.10) Then let [U] = [C]-l[A] v, (3-11) and The equation again takes the form of Equation (3.6), but this time the dominant eigenvalue v gives the smallest value of w. The matrix [C] is positive definite for any stable system; hence [C]-1[A] exists. If eigenvalues other than the dominant one are needed, a new matrix can be constructed which has the same eigenvalues as the original matrix except for the dominant one, which is replaced by zero(l) Repetition of the foregoing iteration then gives the second eigenvalue and eigenvector.

-41For the case where the dominant eigenvalues are complex conjugates and distinct, let a dominant eigenvalue be V = U + iv, and let Wm = VmVm-2 - V2m-1 (3.12) and Zm = VmVm3 - Vm1Vm_2~ In this case Wm 2 2 m -u + v as m -, (3.13) Wml Zm 2Wml Mu as m — o o, (3.14) and (Ym) -Re (km(})) as m -oa, (3.15) where km is a complex scalar. For example, consider the damped system of Figure 4. The equation here is [X;[I + [D'[LE]{v) 0O. (2.78) To find the fundamental mode, which is the mode of lowest frequency, rewrite this as [[El[D] [D] + I] 0. (3.16) Then let [U] = -[E] [D] and (3 * 17) —.

-42 - and the equation takes the desired form of Equation (3.6). Moreover, [0] [' 1 4E]'[D] 0 L~c]-1A] [C]-1[B] and the fact that [C] is positive definite guarantees the existence of -[E]-l[D1. When one puts numerical values into the matrixes, Equation (3.2) becomes 0 0 1 0 0 0 0 1 {Ym Ym} (3.18) -.014151 -.011321 -.009434 -.009434 -.016981 -.023585 +.038679 -.061321 Starting with o= 0 and carrying out the iteration of Equation (3.18), one gets the results shown in Table II. Here again, Vm has been taken as the last element ofYm) From the tabulated values and relations (3.13) and (3.14), we get the dominant eigenvalue v = u + iv = -.01689 +.18131i. (3.19) This gives for the fundamental mode 1 = -.5094 - 5.4677i, (3.20) v which is close to the result found in Equation (2.113).

TABLE II. ITERATIVE SOLUTION FOR COMPLEX EIGENVALUE m 1 2 3 4 5 6 7 8 9 0 -.009434 -.010653 +1.1187x10-3 +3.4581x10-4 -5.1829x10-5 -9.7314x10-6 +2.0619x10-6 +2.5247x10-7 1 -.061321 -.020190 +2.4324x10-3 +5.5118x10-4 -9.6788x10-5 -1.4942x10-5 +3.7027x10-6 +3.7035x10-7 -.009434 -.010653 +.0011187 +3.4581x10-4 -5.1829x10-5 -9.7314x10-6 +2.0619x10-6 +2.5247x10-7 -7.6971x10-8 -.061321 -.020190 +.0024324 +5.5118x10-4 -9.6788x10-5 -1.4942x10-5 +3.7027x10-6 +3.7035x10-7 -1.3529x10-7 Wm - - -.0005568 -1.7045x10-5 -5.3923x10-7 -1.7603x10-8 -5.8163x10-10 -1.9243xl10o- -6.3808x10-13 CWml _ _ -.03061.03163.03265.03304.03309.03316 Wm-l Zm - - 6.501ox10-13

-44The modal column for the fundamental mode is determined from relation (3.15) and Equation (2.82). From these, for large m, {Ym} QIe (km(\f'} = Pe k, i{})(3.21) Let X = (S} = {}I> +;t{P (3.22) and km= m + b Substituting these in Equation (3.21) gives {Ym} (= f - b (P)P(2)} (3.23). (pP) All elements in Equation (3.23) are real, and u and p have been determined, leaving a, b, (0) and {P) to be found. Since the last element of the modal column is arbitrarily assigned unit value, the last elements of (a) and (f) are 1 and 0, respectively. Finding the remaining elements and the values of a and b from Equation (3.23) is straightforward. In the present example, from column 9 of Table II and Equations (3.20) to (3.23), 2. 27 -.50o94 1+ 5.4671. -5094, -5.-677 g 3.7035 }o. [{ -1{. -.7697

-45from which {0}= (X)+ iLp) 587 - 04842i1 (3.25) This agrees closely with the earlier solution, given in Equation (2.115). The value of km merely serves to give unit value to the last element of the modal column. It has no further significance. Unless the driving force is extremely well-behaved, the solution of Duhamel's integral in closed form is likely to be unobtainable. For a piecewise linear driving force, the integral could be evaluated piece by piece, but the number of terms in the solution would increase with each piece added. For any sectionally continuous driving force it would be theoretically possible to expand the driving force in a Fourier series and then integrate term by term to any desired degree of accuracy. However, in many problems encountered in practice, and certainly in earthquake problems, the driving force is too irregular to permit these procedures. Instead, one must resort to approximate methods of integration; and to do so, it is advantageous to work directly with the differential equation instead of Duhamel's integral. For the damped system, the differential equation of interest is h5 -Aq h(t), (3.26) which is the same as Equation (2.99) except that the subscripts have been dropped. This equation in complex variable can be replaced by two equations in real variable. Let k = y + iZ X = l + ip (3.27) and h(t) = f(t) + ig(t).

-46 - Substituting these in Equation (3.26), one gets (yr + iz) - (4 + ip)(y + iz) = f(t) + ig(t). (3.28) Consider the real and imaginary parts of this equation separately. Then y - 4y + pz = f(t) (3.29) and z - py - ~z = g(t). This pair of differential equations in real variable is easily solved numerically on the high speed digital computer. One method of solution is given in Chapter IV. Alternatively, the electronic analog computer could be used. For the undamped system, the differential equation for each mode is a second order equation in real variable, which can be resolved into two first order equations for numerical solution. This is a special case of a more general problem treated in Chapter IV, and will not be considered here. All of the methods given in this chapter have been programmed for the IBM Type 650 high speed digital computer.

CHAPTER IV NUMERICAL SOLUTION OF THE DIFFERENTIAL EQUATIONS OF MOTION FOR NON-LINEAR DYNAMIC SYSTEMS In the analytical methods presented in Chapter II, the differential equations of motion for a linear dynamic system are first transformed into uncoupled differential equations, which are then solved one at a time, either analytically or numerically. A second approach to the problem is to solve the original differential equations of motion numerically without transforming them into uncoupled equations. With the aid of high-speed computing equipment, this can be accomplished for moderately large systems. With proper choice of procedure, irregularities in the driving force cause little difficulty. Also, some of the restrictions on the system can be removed. Non-linearities, such as those which might arise from inelastic deformation or Coulomb damping, can be accommodated in fact as well as in theory. One need not require that the stiffness matrix and viscous damping matrix be symmetric, although these are inherently symmetric for any structural system. A method for solving the differential equations of motion numerically is developed below. The method is aimed at the problem of structural vibrations due to earthquake, but it is applicable to other dynamic problems as well. It will be convenient to adopt the usual subscript notation and a summation convention in which a repeated index in any term implies summation over the full range of values which the index can take. For example, n a..x. means j- ax.. -47

and n n aijbjkXk means IL aijbjkxk An index appearing in parentheses implies that the summation convention is suspended for that particular index in that particular term. For example, aij x(j) means ailxl, ai2x2... ainXn (n quantities) n not aijxj. j=1 The differential equations of motion in their original form were a(i) ei + bij xj + cij xj = fi(t) (i=,2...n) (2.4) where ai is the ith diagonal element of [A]. There are many methods of approximating the solution to such a set of equations numerically. While these differ in detail, they all employ the same basic idea. Let h be a small finite increment of time. One takes the known discrete values of xi, xji, and "xi (i = 1,2,.....n) at times t, t-h, t-2h, etc., as required by the particular numerical process, and the values of fi, which are explicit functions of time, and projects these to evaluate xi, xi, and xi at time t + h. In this mariner, one advances step by step through the solution. No one method of numerical integration can be claimed to be the best for all problems. For the particular application considered in this study, it is advantageous to use a single-step method —that is, one which projects to time t + h from the values of xi, xi, and xi at time t and the functions fi, without using the values of the variables at earlier times. This is advantageous first, because it permits changing the time interval h at any step; and second, because it requires no special starting procedure for the initial steps of the solution. The Runge-Kutta

-49procedure is a single-step method well suited to high-speed computers. Its derivation is available in texts on numerical analysis, and will not be given in detail here. In brief, if we have the equations Yi = 0i(Yl'Y2')''Yn't) (i,2,4.1) the formulas for the Runge-Kutta third order procedure are kio' h0i(Y1y 2, ~Yn, t) -kil= hi(Yl + Pklo' Y2 + Pk20'' n + Pkno, t + ph) (4.2) ki2 - h0i(Yl + (q-r)k1O + rkll'" Yn + (q-r)kn0 + rknl, t + qh) ri(t+h) = Yi(t) + tkio + mk il+ ki(43) The relations governing the constants X, m, n, p, q, r are derived by expanding both sides of Equation (4.3) in power series in h, and equating coefficients of the powers up through h3. The relations are +m+ n =l mp + nq = 1/2 (4.4) mp2 + nq2 = 1/3 npr = 1/6, to which there are infinitely many solutions. One of the more useful solutions is 2 = 1/4 m =0 n = 3/4 p = 1/3 q = 2/3 r = 2/3.

-50This makes m = (q - r) = 0 and thus reduces the number of arithmetic operations required. A second solution, due to Conte and Reeves(5), is 2 =.62653829327 m =.85614352807 n = -.48268182134 p =.62653829327 q =.07542588774 r = -. 55111240553. In this solution, 2 = p = (q - r), which leads to a reduced storage requirement for the digital computer. To adapt Equations (2.4) to the Runge-Kutta process, one can write them in the form i a(i) [fi(t)- Cijxj biZ (45) i Z The solution is then straightforward. For small systems, one can carry out the solution with an ordinary desk calculator. However, there are numerical processes other than Runge-Kutta which are better adapted to manual computation. Using a numerical method of solution, one can solve more general dynamic systems than by analytical procedures. It is only necessary to write the equations as i = 1 [fi(t) - Ri (x. x~ x xn, Z1, Z2 ~ ~ n, t)] (4.6) xi = Zi

-51in which Ri is the resistance function. The functions fi and Ri need not be well-behaved functions —they need only be defined everywhere and piecewise continuous. Equation (4.6) is, of course, merely Newton's second law.

CHAPTER V THE EQUATIONS OF MOTION FOR A MULTI-STORY BUILDING SUBJECTED TO EARTHQUAKE A multi-story building in vibration is an extremely complicated dynamic system. It has infinitely many degrees of freedom and is far from linear. Many of its components are erratic in their behavior, even in carefully controlled laboratory specimens. Finding precisely the motion of such a system except by actual measurement is beyond accomplishment. One can simplify the system by means of reasonable approximations to reduce it to a solvable mathematical model; however, in doing this, one must bear in mind that the solution, no matter how precise it may be mathematically, is only the response of a model. In an earthquake, the ground shakes irregularly and sometimes violently in all directions. This induces structural response in six components of motion —three translational and three rotational() This study is limited to a single translational component of ground motion and a single translational component of structural response. One should note that a single translational component of ground motion may induce all six components of structural motion. Even if "rocking" of the structure is prevented or ignored, rotation about a vertical axis will occur unless the structure is symmetric about a plane parallel to the ground motion or about two perpendicular planes. If the structure has a plane of symmetry parallel to the direction of ground motion, and if "rocking" is prevented, the response will have a single translational component. This is the case considered here. The mass of an actual structure is distributed over the entire structure. In most cases, one can approximate the actual system reasonably -53 -

-54well by assuming that its entire mass is concentrated at several discrete points. Thus one obtains a "lumped mass" system. If the system consists of n point masses, it has 3n degrees of freedom in a general motion, and n degrees of freedom in the one-directional translation considered here. A multi-story building does in fact have most of its mass situated at the floor levels. One can consider it as a lumped mass system, taking the mass of each floor and half the mass of the walls and columns in the stories immediately above and below as a point mass located at floor level. The error introduced in this approximation can be expected to be small in the lower modes of vibration 17) which are the modes of primary importance in structural response to earthquake. Three components of the resistance function R (Equation (4.6) ) are considered in this chapter, namely, flexural resistance of the frame, shear resistance of the walls, and viscous damping resistance. It is well known that the flexural resistance of an elastic frame can be expressed in the form of a stiffness matrix c such that Qi cijxj. (5.1) In Equation (5.1), the forces Qi(i=1,2,... n) are the forces necessary to hold the frame in static equilibrium in the configuration defined by the lateral deflections xj(J=1,2,... n). The j-th column of the matrix c represents the set of static lateral forces that would be required to hold the frame in equilibrium with unit lateral deflection at the j-th floor and zero lateral deflection at all other floors. A conventional static analysis is sufficient to evaluate the stiffness matrix. The effect of dead load on the stiffness (the "overturning" effect) is usually negligible, but can be taken into account in evaluating the matrix if desired.

-55If the lateral deflections become sufficiently large, the elastic limit may be exceeded at one or more locations in the frame. This occurrence will, of course, affect the flexural resistance. For this study, it is assumed that the frame is rigid-jointed; that is, there is no relative rotation between the ends of the members at at a joint. In reality, this assumption is not strictly true except perhaps in a frame in which each member at a joint is welded to develop the full moment capacity of the member. It is further assumed that all members in the frame, columns and girders alike, have ideal elastic-plastic moment-rotation characteristics. The typical moment-rotation curve follows a linear elastic branch until the moment reaches the plastic hinge moment. Upon further increase in strain, the moment remains constant and a plastic hinge forms, making a "kink" in the member. Upon reversal of strain after reaching plasticity, the moment-rotation curve follows a path parallel to the original elastic branch. It remains on the second elastic branch until the plastic hinge moment, either positive or negative, is reached, and then follows a plastic branch as before. Extreme reversals cause the moment-rotation curve to follow a hysteresis loop. This idealized behavior is illustrated in Figure 5. As soon as plasticity is reached at any location, the frame behaves just as though the member were hinged at that location, with a constant moment applied to the hinge. This condition prevails until the next hinge forms or until a decrease in strain occurs at an existing plastic hinge location, causing the hinge to disappear. Cohen, Levy and Smollen(ll) developed a procedure for adapting the method of normal modes to an elastic-plastic frame with infinitely rigid girders, and Schenker(6) indicated a method of extending this to

-56 - z 0 iI R OTATION Figure 5. Idealized Moment-Rotation Diagram

-57cover the frame with flexible girders. In both procedures, an elastic solution is carried out by the method of normal modes to the point where plasticity is reached at some location. The frame is then modified by inserting a real hinge at the plastic hinge location, the new normal modes and frequencies are computed, and the solution is continued as a superposition of the conditions at the instant of transition and the changes found by solving the modified frame. Each time a new hinge forms or an old one disappears, the frame is modified accordingly. In this way a continuous solution can be found. For a step by step numerical solution, a comparable procedure can be devised. The conditions at the end of a time step can be found as a superposition of the conditions at the beginning of the step and the changes during the step. To find the changes during the time step, one first computes the "elastic increments"; that is, the changes which would be caused by the motion during the time step if the frame were rigidjointed and elastic. The resulting moments may exceed the plastic hinge moments at one or more locations in the frame. To correct for this, "corrector" solutions are superimposed in which the frame has real hinges at these locations, and hinge rotation moments are applied at these hinges in such a way that the total moments, obtained from superposition of -the elastic and corrector solutions, nowhere exceed the plastic hinge moments. The superposition of elastic and corrector solutions is illustrated in Figure 6. A hinge rotation at one location affects the forces and moments throughout the frame, and one must find the hinge rotation moments which meet all the requirements of the elastic-plastic system. An iterative imethod of solution is developed below.

-58H2 (a) ELASTIC (c) CORRECTOR NO. 2 KINKS (b) CORRECTOR NO.I (d) SUPERPOSITION OF a, b, AND c. Figure 6. Correction for Plastic Hinge Rotation

-59The elastic increments of force and moment are linear functions of the deflections; and the increments of force and moment due to hinge rotation are linear functions of the hinge rotations. Consider now the variation of forces and moments with time. Let the subscripts k and 2 denote locations in the frame, such that each value of the subscript denotes the end of a specific member at a specific joint. Let Mk be the end moment at location k, taken as positive when it tends to rotate the end of the member clockwise. Let Oi be the hinge rotation at location i, taken as positive when it corresponds to clockwise rotation of the end of the member; and let H2 be the hinge rotation moment at location X, defined as that moment which would cause a real hinge at 2 to rotate through the angle O, if the rest of the frame remained rigid-jointed and elastic. H, is positive when it tends to induce positive hinge rotation. Hi is proportional to ). In this notation, Qi = Qi(x, ) Qi=(x,H), (5.2) and Qi = E xj + I Hii (5-3) Because the Q's are linear functions of the x's and H's, EQ and EQ are constants. In fact, ih = cist (5.4) which is the elastic stiffness matrix. Define =Qi di,, ( )

-60which will be called the force-rotation matrix. Then Qi z cijij + dishie (5.6) Similarly, Mk = Mk(x,H) (5.7) and =Mk. k( ~' a;; Xj s+ +. (5H8) aMk aMk Again, because of linearity of the functions, vx- and. - are constants. J 2 Define 6Mk E~ = [Lk;} (5.9) which will be called the moment-deflection matrix, and aMk Hk = Vkgt (5.10) the moment-rotation matrix. Then Mk = HkjXj + VkH. (5l11) The M's, M's, and H's at every location must satisfy the constraints of the elastic-plastic system, namely, 1. The moment M cannot exceed the plastic hinge moment P in magnitude. 2. If a positive plastic hinge exists, the moment must be either constant or decreasing. If the moment is constant, the hinge may rotate in a direction consistent with the moment. If the moment is decreasing, the -hinge cannot rotate. 3. Equivalent conditions exist for a negative plastic hinge. 4. If the moment is less than the plastic hinge moment in magnitude, the hinge cannot rotate.

These conditions can be expressed thus: 1. -Pk ~ Mk Pk. 2. If Mk Pk' Mk 0, Hk 0, and kk = 0.(5.12) 3. If Mk = -Pk' Mk> 0, k k 0, and MkHk =. 4. If -Pk < Mk < Pk' Hk =. We will have use for the following uniqueness theorem: Given any rigid-jointed rectangular plane framework of ideal elastic-plastic flexural members, and given any x's and any M's not exceeding the plastic hinge moments, there exists one and only one set of M's and H's which satisfies Equation (5.11) and Constraints (5.12). The proof follows. First, observe that if, for one or more k, -Pk < Mk < Pk' then Hk = 0 and there is no constraint on Mk. Thus the solution for Hk is unique, it does not affect other locations, and the effect of conditions at other locations upon Mk is immaterial. Hence it is sufficient to consider only those locations for which IMkl = Pk' Define Yk = - (Sgn Mk)Mk bk = (Sgn Mk)kjxj (513) (5.13) Ski = (Sgn Mk)(Sgn M,)vkY ZY = -(Sgn Me)H2. Equation (5.11) then can be written SkiZ' - kYl = bk, (5.14) where 6k, is the Kronecker delta, 5kR = 1, k = 2, 6kQ= O, k 7.

-62Constraints (5.12) become Yk > 0, Zk )> 0, (5.15) Ykzk 0=. The vector b is unrestricted. To establish uniqueness, it must be shown that the matrix i has properties which guarantee that corresponding to any vector b there exists one and only one pair of vectors y and z that satisfies Equation (5.14) and Constraints (5.15) This amounts to a partition problem in n-dimensional Euclidean space. Consider the columns of the matrixes -6 and i as vectors 7!2''' I~n; ln 2' n. That is, let ({i} - {-ji } and () = ~(ji' Then let al,,2... An be a set of vectors such that every ai is either Qi or fi. There are 2 such sets. If the r's and I's partition the space, every vector b in the space can be formed by linear combination of one andonly one set of ai's with non-negative coefficients.* In other words, for any given vector b there is one and only one set of ails for which b} qi {ci}, (5.16) where every qi is non-negative. If this is so, for those i for which ~i = yi i = qi and zi = O; and for those i for which i = i' i = 0 and zi = qi. Thus the vectors y and z corresponding to any given vector b will be unique and will satisfy the constraints. By virtue of a theorem due to Samelson, Thrall and Wesler, given in the Appendix, necessary and sufficient conditions for the * If any coefficient is zero, say q = O, obviously we can have either am = ~m or A = Cm' However, in this case Ym = Zm = 0 and uniqueness is preserved.

-63columns of 5 and -6 to partition the space are that every principal minor of the matrix i be positive. Consider now the frame with real hinges at one or more locations, in the equilibrium configuration. This is the frame used for the "corrector" solutions, above. Let the hinged joints in this frame undergo hinge rotations A, with lateral deflection of the floors prevented. The resulting strain energy in the frame is E =1/2 OkMk (5 17) The moments are linear functions of the hinge rotations, so we can write MkPkRe (5.18) in which p is a moment-rotation matrix (not identical with v). Combining Equations (5.17) and (5.18) yields E = 1/2 OkPk2eQ* (5.19) We have defined Hk as the moment that would be required to rotate the hinge at k through the angle Ok with all other joints rigid and the frame elastic. Hk and Ok are directly proportional. Thus Ok = (k)Hk (5.20) and Mk =() kQQ' (5.21) where each PQ is a positive constant. Also, since we started from the equilibrium configuration, Equation (5.10) gives us for this case Mk R VkRHi, (5.22) and this with Equation (5.21) makes VkR = P()Pk. (5.23) Because Equation (5.19) is a strain energy equation, it is a positive definite quadratic form. The principal minors of p are therefore positive. Since v can be obtained from p by multiplying each row of p by

-64a positive constant, the principal minors of v must also be positive. Finally, since ~ was obtained from v simply by changing the signs of the off-diagonal elements of certain rows and the corresponding columns (see (5.13) ), the principal minors of i are positive. It follows from the partition theorem that a unique solution exists. The solution can be found by the following iterative scheme. Let zP be the p-th approximation to the vector z, with all elements of zP non-negative. (Zero is a convenient and satisfactory first approximation.) For convenience in notation, define kO = ~kn+l = Zo = Zn+l = 0. Starting with k = 1, and abandoning the summation convention temporarily, let k-l n+l +~=- bk - k Z p+l (5.24) k=O ki =k+l Then, pl P+l p+l uk if +1> O0 zl (not summed over k) kk (5.25) and if ul 4.o, zkp+l _ Repeat this for k=2, k=3,... k=n. Then let bk, (k=1 (5.26)+l kkaZ - bk, (k=l,2,... n). (5.26) The (p+l)st approximation to the vectors y and z is then complete. It satisfies Equation (5.14) and the constraint zk > 0, but may fail to satisfy the other constraints in (5.15). The process is repeated until the amount by which the approximation fails to satisfy all the constraints is insignificant.

-65 - This is, in essence, an adaptation of the Gauss-Seidel iterative method of solving linear equations (18) It can be seen that Equations (5.14) and (5.15) reduce finally to a system of linear equations, since at least one of zkyk must be zero for each value of k. The non-zero elements form a system of n (or fewer) linear equations in the same number of unknowns. In this application, the Gauss-Seidel procedure is used to set up the equations as well as to solve them. The procedure lends itself well to machine computation, and it turns out that convergence is quite rapid. In developing the above procedure, we have dealt with rates of change of deflection, hinge rotation, and moment; that is to say, we have worked with differentials. As soon as one introduces finite differences in place of differentials, one loses the claim to uniqueness, since it cannot be claimed that reversals in strain do not occur within the finite time interval. However, by taking the time intervals sufficiently small and making the further assumption that strain rates are everywhere monotonic within any one time interval, the exact, unique solution of the differential system can be approximated to any desired degree of accuracy by the finite difference system. To set up the finite difference system, let the subscripts s and s+l denote a condition at the end of the s-th and (s+l)st time step; e.g., xi s+l is the value of the variable xi at the end of the (s+l)st time step. Also, let xi xi,=x - xi i i Xi,s+l Xi,s' Q~k = ks- Mis' (5.27) SMk = Mk,s+l ks Hk,s+l k,s'

-66 - Then the finite difference equations are ~Q i = cijAxj + diAH- (5.28) and AMk = IkjAXj + Vk/AH2, (5.29) which correspond to Equations (5.6) and (5.11).. The elastic-plastic constraints are: if 0,S Pk either AMk = Pk - Mk,s and 0Hk 0 (5.30) or -Pk - Mk,s < Mk< Pk Mk,s and AHk = and if Pk Mk,s ~ either'Mk = - Pk - Mk,s and 0Ak > 0 (5.31) or - Pk - Mk,s < AMk < Pk - Mks and k = O In the finite difference system, by taking Axj sufficiently large, one can construct mathematically a situation for which no solution exists. For example, if, for some value of k, Mkjs were negative, one could choose the Ax. so that ukjAxj would be a very large positive quantity, so large that a solution of Equation (5.29) would require either that LAMk exceed Pk - Mk,s or that AHk be negative, both of which are in violation of Constraints (5.31) for negative Mk,s. The occurrence of such a situation in the course of an actual numerical integration would mean that very large changes in moment were taking place within a single time step — a condition that would rapidly destroy the accuracy of the numerical results. Should such a situation arise, it simply means that the time interval was taken too large in the first place. It is well known that the walls in a multi-story building may have a sizable influence on the response of the building to an earthquake.

-67However, the exact behavior of walls is very difficult to evaluate. In experimental work on determining the periods of vibration of the fifteenstory Alexander Building in San Francisco, Blume found that the tile partitions offered little or no contribution to the rigidity of the (19) building at small amplitudes of vibration 9) The brickwork participated in the rigidity, but not to the extent that its modulus of elasticity would imply. A greater participation might be expected at the greater amplitudes which would be experienced in a strong motion earthquake. Housner points out that a large share of the effect of masonry walls in large amplitude vibrations is due to the energy dissipated by the masonry units cracking and grinding together (20). This is, in effect, a modified form of Coulomb friction damping. Even plaster partitions may be able to dissipate appreciable amounts of energy this way. Experiments conducted at Stanford University on brick and concrete shear walls with and without openings(21'22) suggest that the loaddeflection characteristics of such walls can be reasonably described by a bi-linear curve of the elastic-plastic type. Tests involving strain reversal in masonry shear walls have not been reported, to the author's knowledge. For reasons of mathematical convenience, and in the absence of better information, it will be assumed that an elastic-plastic type load-deflection curve, with hysteresis upon reversal of strain, is a reasonable representation of the overall behavior of a shear wall. It is emphasized that no claim is made that masonry has elastic-plastic properties in shear.

-68With the assumption of elastic-plastic type behavior, the shear resistance of the walls is not difficult to evaluate. Let Vi = shear resistance of the walls in the i-th story (the story below the i-th floor), Gi = shear rigidity of the walls in the i-th story, Ui = shear strength of the walls in the i-th story, and ~i = xi - xi-1 = deflection within the i-th story. Considering the variation of V with time, and the constraints of elasticplastic behavior, one gets the following system:.If -Ui = Vi and i; 0, Vi = 0. If -Ui = Vi and ii > 0, or if -Ui < Vi < Ut J Vi G(i),ti (5.32) or if +Ui = Vi and ti < O, If +Ui = Vi and ti 0/ O, Vi = 0. As before, we can approximate this by a finite difference system, taking the time interval sufficiently small and making the assumption that strain reversals do not occur within the time interval. In finite difference form, let X~i = hi,s+l - hisThe system then is: If Vi + G )~(i > i U V i U. If -Ui < Vi,s +G(i)Si Ui V,+l V5 + G(i)i Ui Vis+ Vi, + G (5.33) If Vis + G (i) -Ui Vi,s+l = -Ui The resistance force at the i-th floor due to shear walls is Vi - Vi+1. Note that if the frame is assumed to have infinitely rigid girders, and if the columns in any story are elastic-plastic and all yield simultaneously (the usual elastic-plastic shear building

-69assumptions), the flexural resistance of the frame can be found by the procedure given above for shear walls. This is much simpler than the method outlined earlier for a general elastic-plastic rectangular frame. Not a great deal is known about the damping forces in actual structures. Customarily, it is assumed that damping is viscous; that is, the damping forces are linear functions of the velocities. There are two reasons for this. First, viscous damping is mathematically convenient. Until the advent of the high speed computer, it could have been said for all practical purposes that the assumption of viscous damping was a mathematical necessity, for any other assumption led to non-linear equations which could be solved only in a restricted sense(3"4). Second, in the course of small amplitude shaking machine studies it has been observed that the shaker force and the building amplitude tend to be nearly proportional, which is characteristic of a linear system with viscous damping. This has been taken not as a sign that damping forces are truly viscous, but rather that on the whole, a building behaves similarly to a system with linear damping at small amplitudes of vibration(23). The observations were made at amplitudes of the order of one per cent of the amplitudes that might be expected in a strong motion earthquake. Whether they would hold for much larger amplitudes has been neither established nor refuted. Hisada and Nakagawa have conducted a number of large amplitude vibration tests in Japan with the use of a shaking machine (24) They report that the effective viscous damping coefficient tends to increase with an increase in amplitude, due to progressive structural damage. In building vibration tests, the driving force produced by the shaking machine is sinusoidal. Jacobsen has shown that for a single

-70degree of freedom system subjected to a sinusoidal driving force, the assumption of viscous danmping can be made without appreciable error, provided the coefficient of viscous damping is chosen so that the energy dissipated in each cycle is the same in the viscous damped system as in the actual system(5). The extension of this concept to a system having many degrees of freedom and subjected to a random driving force is not immediately clear. If viscous damping is assumed to exist, the damping forces are linear functions of the velocities, and can be expressed as iD =t 4. (5.34) i ijxj' in which Di is the damping force at the i-th floor, b is the viscous damping matrix, and x is the velocity of the j-th floor relative to the ground. The matrix element bij is the damping force that would exist at the i-th floor if the structure were passing through the equilibrium configuration with unit velocity at the j-th floor and zero velocity at all other floors.

CHAPTER VI CALCULATED RESPONSE OF MULTI-STORY BUILDINGS TO EARTHQUAKE GROUND MOTION To illustrate qualitatively the effects of variations in the dynamic system upon the response, and also to demonstrate that the methods developed in preceding chapters can in fact be carried out, response calculations have been made for buildings of two different heights, and for two different earthquakes. The buildings chosen for the calculations represent steel framed structures proportioned for the lateral wind loads specified in typical United States building codes. The general framing diagrams are shown in Figures 7 and 8. The frames are not designed for seismic forces, and are more flexible than might be expected of seismically designed buildings of comparable size. The mathematical models used f6r the response calculations are based on frames one bay wide, "equivalent" to the three-bay frames in the sense that the total flexural stiffness and ultimate flexural strength of all the columns in each story are the same in the "equivalent" frame as in the original frame. The same equivalence applies also to the beams in each floor. The reason for this substitution is that in the elasticplastic case, computing the moments and hinge rotations for the original three-bay frames would have required more storage space for matrix elements and computed results than was available in the 2000-word memory of the computer. Even though the original three-bay frames could have been used for all cases except the elastic-plastic case, the "equivalent" frames were used for all analyses in order to avoid introducing unwanted differences which might invalidate comparison of the results for the -71

72 Symm. about t 12VWF27 12 UF27 I=708.3 in4 * C 14VWF 30 14 30 =964.8in4 oo -D04 - N d _______ 14,V 3O 14 W 30 do 4 F 3O 1 14 F 34 | I I I do oNN O- w C, 14W:'30 14W 30 do i1 16 W 36 16 36. I=1434.9in4*. *? I 3 @. 20'-0" = 60'-0" 3 20'- 0" Elevation of typical bent Elevation of typical bent Original Frame Equivalent Frame * Includes Iof floor slob Figure 7. Framing Diagram, Eight-Story Bent.

Symm. about _ 12 VF 27 12 VF 27 I=708.3 in4*( 0 __ __14 VF 30 14 VF 30 C I=96 doin.... 14_W' 30O 14 W 30 N do OD) C 0 Elevation of typical bent Elevation of typical bent Original Frame Equivalent Frame * Includes I of floor slob Figure 8. Framing Diagram, Four-Story Bent.

-74different systems. It turns out that the equivalent frames are slightly stiffer than the original frames, the difference being of the order of one per cent. For the four-story building, calculations were made using four different mathematical models to represent the structure. The first of these was an undamped elastic system in which the resistance was taken as that due to the elastic frame alone,\with no allowance for damping forces, inelastic frame deformation, or the resistance of walls and partitions. Second was a damped elastic model in which a resistance due to interfloor viscous damping was added to the elastic resistance. In the third model, the system was taken as the undamped elastic resistance of the bare frame, plus an elastic-plastic type component of resistance due to shear walls. The shear wall component was purposely kept small to avoid its dominating the response. For the fourth model, the system was taken as an elastic-plastic bare frame, without shear walls or viscous damping. For each of these systems, the elastic resistance of the bare frame was the dominant resistance element, and in the last three cases it was modified by one additional component of resistance. For the eight-story building, the analyses were limited to the elastic undamped case and the elastic frame with inter-floor viscous damping. For all analyses it was necessary to have the stiffness matrix c (Equation (5.4) ); and for the elastic-plastic case, the matrixes d, X and v (Equations (5.5), (5.9) and (5.10) ) were also needed. These matrixes were evaluated on the high speed computer by Mr. D. A. DaDeppo, Research Assistant, Engineering Research Institute, University of Michigan. The computer program was written so that the instructions, member stiffnesses, story weights, and story heights were first loaded into the

-75 - memory. The computer then wrote the slope deflection equations, solved them by iteration, used the results to evaluate the matrix elements, and printed out the elements on punched cards. For the viscous damping matrix, b, for each building, it was assumed that inter-floor viscous damping forces of a magnitude of 3 kip sec. in 1 were effective in each story. This type of damping is shown symbolically in Figure 9. The matrixes for each of the buildings are given in Tables III and IV; and the index numbers are identified in Figure 10. For the four-story building with shear walls, the same elasticplastic type resistance function was used for the shear walls in each story. The shear stiffness in each story was taken as 30 kips per inch and the ultimate shear strength in each story was taken as 20 kips. The shear wall stiffness was purposely taken to be less than the stiffness of the elastic frame; and the shear wall strength was taken to be quite low. This was done so that one might be able to get an idea of how the response is affected by adding a small amount of shear wall resistance to the stricture. If the shear walls had been assumed to be both very stiff and very strong, the shear wall component of resistance would have dominated the response and invalidated comparison of results. Two earthquake accelerograms were used in the analyses; namely, the E-W component of the earthquake recorded in Helena, Montana, on October 31, 1935, and the N-S component of the earthquake recorded at El Centro, California, on May 18, 1940. The accelerograms were assumed to be true records of the actual ground accelerations. It was found that the recorded accelerograms could be approximated quite closely by piecewise linear functions, which are also convenient for use in the computer.

3 kip sec in I 3 k i p sec in i 3 kip sec in' Figure 9. Interfclor Viscous Dping

TABLE III. MATRIXES, FOUR-STORY STRUCTURE Elastic Stiffness Matrix Moment-Deflection Matrix 136.71 -84.63 9.38 -0.63 -.64053.09448 -.00631.00042 cii -84.63 109.97 -43.75 4.44 kip in1 -.65613.18897 -.01261.00085' 9.38 -43.75 68.94 -33.63.72401 -.59992. o04003 -.00270 0.63 4.44 -33.63 29.76 -.05075.30728 -.02051.00138 Uk j] = |.6o432 -.64oo9.10277 -.00692 in1 -.23923 ~87376 -.70583.04751 Viscous Damping Matrix -.37514.19840.14951 -. 01006 -.10560 ~ 74233 -.75227.10503 6 -3 0 0.01764 -.12401.72314 -.61502 [biJ] l 1-3 6 -3 0 kip sec in-1 L.00725 -.05095.62902 -.58460 0 -3 6 -3 0 0 -3 3Force-Rotation Matrix -80.50 -101.31 104.47 -10.04 91.28 -13.29 -87.33 - 6.01 1.03 0.45 1 11.87 29.18 - 86.56 60.79 -96.68 48.54 46.19 42.27 -7.22 -3.19 kip - 0.79 -1.95 5.78 - 4.o6 15.52 -39.21 34.80 -42.83 42.09 39.36 0.05 0.13 -0.39 0.27 - 1.04 2.64 - 2.34 5.98 -35.80 -36.58

TABLE III. MATRIXES, FOUR-STORY STRUCTURE (Cont'd) Moment-Rotation Matrix i. 0000.5000 -.1575 -.3176 -.0614.oo0089.0587.00oo40 -.0007 -.0003.4070 1.0000 -.3150 -.6351 -.1227.0179.1174.0081 -.0014 -.0006 -.1372 -.3371 1.0000. -.7023.3896 -.0567 -.3727 -.0257.00oo44.0019 -.2017 -.4957 -.5122 1.0000 -.1995.0291.1909.0131 -.0022 -.0010 [vki] = -.0510 -.1254.3721 -.2613 1.0000 -.1456 -.9567 -.0659.0112.0050.0202.0497 -.1473.1035 -.3959 1.0000 -.8876.4524 -.0773 -.0341.0317.0779 -.2310.1622 -.6208 -.2118 1.0000 -.0958.0164.0072.oo0089.0219 -.0o650.0457 -.1747.4414 -.3918 1.0000 -.1708 -.0754 -.0015 -.0037.0109 -.0076.0292 -.0737.o655 -.1671 1.0000.4416 -.000 ooo6 -.0015.00oo45 -.0031.0120 -.0303.0269 -. 0686.4108 1.0000

TABLE IV. MATRIXES, EIGHT-STORY STRUCTURE Elastic Stiffness Matrix 444.43 -295.77 71.77 -11.84 1.60 -0.22 0.02 -0.00 -295.77 421.08 -244.30 61.38 -8.32 1.14 -o0.0o8 0.01 71.77 -244.30 346.55 -203.42 44.74 -6.12 0.43 -0.03 [Cil= -11.84 61.38 -203.42 270.54 -146.44 30.15 -2.12 0.15 kip in-1 1.60 -8.32 44.74 -146.44 209.02 -112.48 13.04 -0.95 -0.22 1.14 -6.12 30.15 -112.48 133.52 -51.62 5.60 0.02 -0.08 0.43 -2.12 13.04 -51.62 78.06 -37.72 -0.00 0.01 -0.03 0.15 -0.95 5.60 -37.72 32.94 Viscous Damping Matrix 6 -3 0 0 0 0 0 0 -3 6 -3 0 0 0 0 0 o 0 -3 6 -3 0 0 0 0 [b1i] 0 0 -3 6 -3 0 0 0 kip sec in' 0 0 0 -3 6 -3 0 0 0 0 0 0 -3 6 -3 0 0 0 0 0 0 -3 6 -3 o0 0 0 0 0 -3 3

-8010) 0 i 2) ~ 2 ii- i 0 k, Figure 10. Index Identification for Matrixes

-81The time-acceleration coordinates of the intersection points of successive line segments were put in punched card form for machine input data. Each accelerogram was "balanced" by shifting the base line slightly to make the computed ground velocity zero at the end of the earthquake. Each card contained the coordinates of seven points of the accelerogram plus an index number. The index number permitted the machine to check each card for proper sequence before accepting its data for use in the calculation. Only the intersection points of the line segments were given to the machine; and the machine was instructed to interpolate linearly between successive points whenever the ground acceleration for an intermediate time was required. The accelerogram for the Helena earthquake, plotted from a tabulation of the punched cards that were used for machine input, is shown in Figure 11. The figure also shows the ground velocity and ground displacement which were computed by integrating the accelerogram numerically on the computer. The response of the four-story undamped elastic system to the Helena ground motion is shown in Figure 12. It is seen in this figure that after the first 2.5 seconds the response tends to be nearly periodic, with a period of 2.0 seconds. The natural period for the fundamental mode of this system is in fact 2.01 seconds. The irregularities in the deflection curves are evidence of vibration in the higher modes. The response of the four-story viscous damped system to the Helena ground motion is shown in Figure 13. As might be expected, there is a marked difference between this and the undamped system. Here the deflections build up to a peak at about 2.9 seconds, which is about the end of the most severe portion of the ground shock. Thereafter, the structure oscillates back and forth with a period of about 2 seconds, and

-82with steadily decreasing amplitudes. The irregularities in the deflection curves are not as pronounced here, indicating that the higher modes of vibration are damped out rapidly. Since the viscous damped elastic system is in fact a linear system and can be solved by the method of normal modes, it is of interest to evaluate the first mode of response for this system and compare it with the total response obtained by direct numerical integration. The eigenvalue and modal column for the first damped mode are: x = (-0.4587 + 3.0987i) secl (6.1).3060 +.0217i (}..5640 +.0354i (6.2).8637 +.0157i 1.0000 The equations of motion for the first damped mode response are, from Equations (2.73), (2.94), (2.99), (2.100) -nd (2.109), - (-.4587 + 3.0987i)r = (.003259 +.2144i)w (6.3) and {x}= 2 Re } (6.4) Figure 14 shows the first damped mode of response of the viscous damped system to the Helena ground motion. A comparison of this and Figure 13 shows that, at least in this case, the response in the first damped mode compares very favorably with the total response. one can also use the more conventional method of modal analysis and find the modal column and frequencies for the undamped system and then add the damping term in the corresponding equation of motion. The

-83undamped frequency and modal columns are: a = 3.1256 sec-1 (6.5).3027 ({}'-.5592 (6.6).8608 1.0000 Using Equations (2.20), (2.35), (2.36), and adding a damping term, one gets the following equations of motion: t + 0.9253t + 9.7691 = -1.3255;w (6.7) {x}=() ()(6.8) In Equation (6.7) the damping term was chosen to get the same fraction of critical damping in this case as prevailed in the first damped mode. Note that there is no way of relating the fraction of critical damping to the properties of the system except by finding the eigenvalue for the damped system. If the same fraction of critical damping is assigned, 0.148 in this case, the response in the undamped mode is nearly the same as the damped mode response. The results are not shown graphically. Figure 15 shows a plot of the response to the Helena quake for the four-story undamped structure with shear walls added. The shear walls act as conservative elements when they are in the elastic range, and become dissipative elements when the elastic range is exceeded. The effect of the conservative forces from the shear walls appears as a shortening of the period of vibration of the structure. After the structure reaches its extreme amplitude, it oscillates with a period of 1.45 seconds, or about three quarters of what the period was when the shear walls were not present. The fundamental period of the system is 1.45 seconds

when the shear walls are in the elastic range. The effect of the dissipative forces due to shear walls appears as a decay in the amplitude of oscillation with time. However, the amplitude here does not approach zero as it did in the case of the viscous damped system, but rather diminishes to the point where the oscillations have an amplitude sufficiently small to permit the shear walls to remain in the elastic range and, thereafter, the system behaves as a conservative undamped system. The dissipative forces here do not tend to damp out the higher modes of vibration. No response is given for the elastic-plastic frame to the Helena ground motion, because the ground motion was not severe enough to induce plasticity in the structure. Figure 16 shows the deflectiontime curves for the four-story undamped elastic structure responding to the El Centro ground motion. The El Centro earthquake was considerably more intense and of longer duration than the Helena earthquake. Only the first ten seconds of response is shown in this curve. The characteristics of the response are much the same as appeared in Figure 12 for the same system and -the Helena ground motion. The amplitude of oscillation increases steadily during the time shown. The calculations were actually carried further; and it was found that the maximum amplitude was about ten per cent greater than that shown on the figure. The response shown here was computed on the assumption that the system would remain elastic regardless of the amplitude of the deflections. This would be the case if the frame members were built of steel having an extremely high elastic limit. ASTM-A7 steel has an average yield strength of approximately 42,000 psi under dynamic load conditions, and would be strained well into the inelastic range by the oscillations shown in Figure 16.

-85For the same framework, with elastic-plastic frame members of 42,000 psi yield strength, the response to the El Centro ground motion is plotted as deflection vs. time in Figure 17. At first glance, these curves appear almost unrelated to the curves for the elastic system, except for the first 2.5 seconds of response. However, a close inspection reveals a great deal of similarity. The periodicity of the oscillation is very nearly the same as for the elastic system; but in the elasticplastic system, the oscillation occurs about a position of equilibrium which shifts away from the zero position and keeps changing as inelastic deformation progresses. Moreover, the shift in equilibrium position is not the same for all floors. Also, the amplitude of oscillation about this shifting equilibrium position remains nearly constant; whereas in the elastic system, the amplitude increased steadily during the first ten seconds of the response. Although the total deflections, measured from thd zero position, are greater for the elastic-plastic system than for the elastic system, the amplitudes of the oscillation are only about half the elastic amplitudes. Here, as in the earlier case where shear walls were considered, the system behaves as a conservative system during the part of its response when the system is elastic, and behaves as a dissipative system when plasticity is present. The occurrence of plasticity induces a decrease in stiffness of the structure, which accounts for the increase in maximum deflection measured from the zero position. Plastic deformation also dissipates energy that would otherwise be stored in the structure as recoverable strain energy, and thus acts to reduce the amplitude of oscillation. The locations at which plastic hinges formed during the response are shown in Figure 17.

-86 - The maximum values obtained for various parameters in the analyses are summarized in the tables below. The important parameters are the maximum deflections at each floor; the shear coefficient for each story, which is defined as the maximum shear in the story divided by the weight of the structure; the shear distribution, which is the shear coefficient for any story divided by the shear coefficient for the bottom story; and the seismic coefficient for each story, which is defined as the maximum shear in the story divided by the weight of that part of the building situated above the story. These parameters are shown in Table V for the Helena ground motion'for the four-story building. Table VI shows the same parameters for the four-story building and the El Centro ground motion, and Table VII shows the parameters for the eight-story building and the Helena ground motion. For comparison, the design parameters computed in accordance with the Los Angeles and San Francisco building codes are also tabulated.

GROUND ACCELERATION, fraction of gravity GROUND DISPLACEMENT, inches GROUND VELOCITY inches per second II 6 b 1 c~~~~~~~~~~~~~~~~~~~~~~~~~ N~~~~~~~~~~~~~~~~~~~~~~~~ - - N, C M A A N ar N lU__ _ _ _ _ _ _ _ I I I ru I I I~ LO N CM CM~~~~~~~~~~~~~~ -I m~~~~~~~~~ 0 CD - S Figure 11. East-West Component of Earthquake at Helena,, Montana, October 51, 1955..7

5 4 3 Q) LuV TIME, seconds. — J ~~~~6-37 8 I 0 c -4 Figure 12. Response of Four-Story Undamped Elastic Frame, Helena Earthquake.

4 -c 0 U IJ LU cc — 2 -- o 1 2 3 4 5 6 7 8 9 10 TIME, seconds Figure 13. Response of Four-Story Viscous Damped Frame, He lena Earthquake.

4 LJ~ 0 -c2 z (.)-'C - 0 4 5 678 91 TIM\^E, seconds. Figure 14. First Mode Response of Four-Story Viscous Damped Frame, Helena Earthquake. 0 C~ n —I w 0 IIiIII I I I, 0 I 2 ~~3 4 5 6 7 8 9 10 TIME, seconds. Figure 14 First Mode Response of Four-Story Viscous Damped. Frame, Helena Earthquake.

4 3 z _o u _0 -3 _ -I lO o 1 2 3 4 5 6 7 8 9 10 TIME, seconds Figure 15. Response of Four-Story Undamped Frame with Shear Walls, Helena Earthquake.

15 u O 10 0 -c tC 5 -J O W r TIME, seconds. Figure 16. Response of Four-Story Undamped Elastic Frame, E1 Centro Earthquake.

10 PLASTIC HINGES FORMED AT THESE LOCATIONS DURING 2. 5 RESPONSE I Z o o IL(~~~~~~~~~~~~~~~~~~~~~~~~N.) LU /, -~~~~~~~~~~~.J~.., /.'~ ~ L-1 1 ~~~ ~~ ~~~~~~~~~~...~ x."''''; trQ~~~ ~~~~~~~~~~~~~~~~~~~~~...'.~. LL~ ~ ~ ~ ~ ~ ~ ~ ~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~'. /.T:' 1 -5 -2"\.''~' ~~.,' I/~~~~~~~ -15 ~3 i / I~~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~~~~~~~~~~'" I. " II 0 I 2 4 5 67 8 9 TIME, seconds Figure 17. Response of Four-Story Elastic-Plastic Frame, El Centro Earthquake.

TABLE V. MAXIMUM MAGNITUDES OF RESPONSE PARAMETERS, FOUR-STORY STRUCTURE, HELENA GROUND MOTION Viscous Damped Frame Frame San Undamped and Total First First Los FranParameter Elastic Shear Response Damped Undamped Angeles cisco Frame Walls Mode Mode (4) Code Code x4 4.13 3.76 3.13 3.11 3.16 X3 3.50 3.30 2.70 2.70 2.72 H r X2 2.39 2.32 1.78 1.78 1.77 0) X1 1.28 2.32.97.97.96 H s4.032.019.014.013.014.024.018'"a S3.062.038:.032.031.033.049.038 hd 9 S2.074.059.046.045.045.o66.052 ~H a)S.073.074.053.053.052.080.060 En D4.44.25.27.25.27.29.29 c N~ D3.85.51.61.60.63.61.63 fl D2 1.01.80.88.86.87.83.87 0 )wO D1 1.00 1.00 1.00 1.00 1.00 1.00 1.00 ~ C4.181.106.080.074.080.133.099 u.v' C~.139.o84.072.070.073.109.085 C_2.103.082.065.063.063.092.073 ( Cl 1 01.073.074.053.053.052.o80. o60 Max. Shear in ith Story (1) Si Weight of Structure th Max. Shear in i Story (2) Di = Max. Shear in First Story Max. Shear in ith Story (3) Ci = Weight of Structure above ith Story (4) Damping Term added to equation to give same fraction of critical damping (14.8%) as for first damped mode.

TABLE VI. MAXIMUM MAGNITUDES OF RESPONSE PARAMETERS FOUR-STORY STRUCTURE, FIRST TEN SECONDS OF EL CENTRO GROUND MQTION Viscous Damped Frame Undamped Elastic Undamped Elastic- Frame and Total First First Los San Parameter Elastic Plastic Shear Response Damped Undamped Angeles Francisco Frame Frame Walls Mode Mode (4) Code Code, x4 15.01 19.06 6.35 6.68 6.79 6.92 X3 11.26 17.65 5.57 5.88 5.88 5.96 CH C. X2 7.96 7.25 3.49 4.o6 3.85 3.87 ~ w H x1 4.75 5.20 1.86 2.290 2.09 " S4.133.071.037.029.029.031.024.018 C H S3.219.082.097.068.069.072.o49.038 Md m S2.180.118.112.102.098.098 o66.052 oo. S1.268.127.109.127.113.113.080.o60 D \D4.50.56.34.23.26.27.29.29 3.82.65.89.54.61.63.61.63 f -A D32.67.93 1.03.80.86.87.83.87'' - D1 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 c C4.754.403.207.161.166.175.133.099. -- C3.490.184.217.153.155.161.0log.085 co C2.251.164.156.142.136 137.092.073 o.5 C1.268.127.109.127.113.113.080.o60 (1) S. Max. Shear in ith Story ) C Max. Shear in ith Story i Weight of Structure i Weight of Structure above ith Story Max. Shear in ith Story (4) Damping term added to equation to give same (2) Di Max. Shear in First Story fraction of critical damping (14.8%) as for first damped mode.

-96 - TABLE VII. MAXIMUM MAGNITUDES OF RESPONSE PARAMETERS, EIGHT-STORY STRUCTURE, HELENA GROUND MOTION Viscous Damped Frame Undamped Total First First Los San Parameter Elastic Response Damped Undamped Angeles Francisco Frame Mode Mode (4) Code Code x8 7.35 4.78 4.28 4.36 x7 6.og 4.36 4.03 4.10 x6 5.00 3.47 3.46 3.50 0.- xI 4.81 2.74~ 2.92 2.95 X4 4.25 2.06 2.23 2.25 -c X3 3.35 1.44 1.57 1.59 m x2 2.11. 88. 94 0.95 S H X1 1.04 o. 40 0.42. 43 S8.026.007.oo4. oo4. oll.005,l S.037.017.010.010.023.011 S.055.023.015.016.032.017 oi S5.o46.026 02.0.020.038.021 U S'. S4.039.027.023.024.043.025 * s.057.026.026.026.0o46.028 S3 o64.027.027.027.049.030 Sl.1073.030.028.028.052.031 D8.36.24.15.15 22.15 D.51.56.36.37.45.36 ~ D7. 75.77.55.56.61.54 D5.63 88.71 72.73.69 Q D4.54.90.84.84.82.81 ~.~ D3.78.87.92.92.89.90 + D2.88.90.98.98.94.96 Q D1 1.00 1.00 1.00 1.00 1.00 1.00 c8.312.0o86.o48.050.133.056 O. C7.174.077. o47.049.109.052 ~ C6.158.067.045.046.092.048 ~U U C5 C.097.055. 042.043.080.045 C5 4 - C4.0o65.044.039.039.071.041. C3.078.035.035.063.038.H k C2.074.031.031.032.057.034 C1.073.030.028.028.052.031 (1) Si Max. Shear in ith Story - Weight of Structure (2) Di = Max. Shear in ith Story i x. Shear in First Story (3) C Max. Shear in ith Story i Weight of Structure above ith Story (4) Damping term added to equation to give same fraction of critical damping (6.25%) as for first damped mode.

CHAPTIER VII ENERGY RELATIONS An explanation of the effect of damping on the motion of a dynamic system lies in energy considerations. The earthquake itself has its origin in the release of energy by fault movement somewhere in the crust of the earth. The ground motion in turn feeds energy into the structure. If the structure were perfectly rigid, its kinetic energy at any time would be equal to one half its mass times the square of its velocity, which would be the same as one half its mass times the square of the ground velocity. If the structure were perfectly flexible laterally, the ground motion would impart to it no energy at all, since it would acquire neither kinetic nor strain energy. Fcr any real structure lying somewhere between the ideals of perfect rigidity and perfect flexibility, the properties of the structure affect the amount of energy it will acquire from the earthquake. The energy per unit mass is not limited to one half the square of the ground velocity, but may exceed this value. The energy imparted to the structure takes three forms; namely, kinetic energy, strain energy, and dissipated energy. Energy dissipation arises from damping forces, which include both viscous damping and inelastic deformation. The three forms of energy can be computed by processing the output data obtained in the numerical integration of the equations of motion. The kinetic energy is given by the equation' 2 T = 1/2 ai(i + (71) -97

-98where T is the kinetic energy of the system, ai is the i-th mass, Xi is the displacement of the i-th mass, w is the ground displacement, and the repeated subscript implies summation. If the frame is elastic, the strain energy in the frame, R, is given by R = 1/2 ij i =1/2 Qi i (7.2) where cij is the stiffness matrix, and Q is the resistance force at the i-th floor. If the frame is elastic-plastic, the same equation gives the strain energy until such time as plasticity is first encountered. After plasticity has occurred anywhere in the frame, Equation(7.2)becomes invalid; but the strain energy can then be obtained from the equation 1 i j in which yij is the inverse of ci..; i.e., c o =c 8 (7.4) ij jk ik (.4) where 5 ik is the Kronecker delta. Equation(7.3)is valid whether plasticity has occurred or not. If plasticity has occurred, the strain energy is fully "recoverable" only if the structure can return to a stress-free configuration (not the zero configuration) without additional plastic deformation. Generally, this is not so.

-99If the shear walls have not been deformed beyond the elastic range, the strain energy contained in the shear walls, S, is given by the equation = 1/2 Vi i (7 5) where Vi is the shear resistance force in the i-th story, and = x - x is the deflection in the i-th story. i i i-l If inelastic deformation has occurred, the strain energy is given by the equation 2 V. 1 1 - (7.6) 2 G Gi being the shear rigidity of the walls in the i-th story. Equation(7.6 ) is also valid before inelasticity is reached. From the above equations, the strain energy contained in the structure at any time for any configuration of the structure, within the limitations of this report, can be calculated. Two points should be noted. First, the strain energy is computed from a knowledge of the constant properties of the structure, and the configuration at any instant of time. It does not require a knowledge of the past history of the configuration. Second, the strain energy R for the elastic-plastic frame is the strain energy which would be recovered if the frame could return without further plastic deformation to its equilibrium configuration of that particular instant. Likewise, the strain energy S for the shear walls is the strain energy that would be recovered if the walls returned to their equilibrium configuration of that instant. However, in general, if inelastic deformation has occurred, the equilibrium configuration for the frame and the

-100equilibrium configuration for the shear walls at any instant of time will not be the same. Therefore, the strain energy is not fully "recoverable." The dissipative forces considered in this report are viscous damping, inelastic deformation of the shear walls, and plastic hinge rotation in the structural frame. The rate of dissipation of energy in viscous damping, L, is given by the equation L = Di x (777) where Di is the viscous damping force at the i-th floor. Since the parameters appearing in the right side of this equation appear in the computer output for equally spaced values of time, it is easy to compute the rate of dissipation for each value of time and integrate the results numerically to find the total amount of energy dissipated in viscous damping at any time. The rate of energy dissipation due to inelastic deformation of the shear walls is given by the equations K i i= (7.8) and Vi if lvii = Ui W (7v9) 0i if [vil <u. In the above equations, K is the rate of energy dissipation in inelastic deformation of shear walls, W is the shear dissipation force in the i-th story, Ui is the strength of the shear walls in the i-th story, and other terms are as defined earlier.

-101Because the terms Wi are either zero or constant non-zero values dictated by the properties of the structure, it is more convenient to consider the work done by the dissipative forces in a finite time interval, rather than the rate of dissipation at any instant of time. Using this approach, the energy dissipated by inelastic deformation of the shear walls in a small finite time interval is AK = ZA K (7.10) i i where, dropping the summation convention temporarily, {U(L)(I i 4s l-tJsi- I\sK ~if IV~= U \AK = - (7.11) if i,(7.10) and (.11), In Equations (7.10) and (7.11), 6K is the energy dissipation in the time interval ts < t < ts+1 AKi is the energy dissipation in the i-th story shear walls during the same interval, and the subscripts, s and s+l, denote the values of the parameters at times ts and ts+l. Other terms are as defined above. Here again, since the required data is produced in the computer output at equal intervals of time, it is simple to process the data and obtain the energy dissipated in inelastic shear wall deformation.

-102The energy dissipation due to plastic hinge rotation can be computed in a similar manner. At each location, if energy is being dissipated at all, the dissipative force, in this case having the dimension of a moment, is constant. The rate of energy dissipation is given by the equations J = Nk (7.12) and Nk k if {Mk[ = Pk (7.13) 0 if IMkI < Pk where J is the rate of energy dissipation through plastic hinge rotation,.k is the plastic hinge rotation at location k, M is the moment at location k, and Pk is the plastic hinge moment at location k. Since the dissipative force, Nk, is either zero or constant, it would be more convenient to use an equation which gives the work done in plastic hinge rotation during a small finite time interval. This equation is 6J - Z'Jk (7.14) k where Jj= { P(k)Ik,s+l - Pk,sI if IMk,s+l Pk (7.15) 0 if lMk,s+ll < Pk

-103In these equations, z\J is the energy dissipation in the time interval ts<t<ts+l lJk is the energy dissipated in plastic hinge rotation at location k in the same interval, and other terms are as defined above. Equation(7.15)differs in form from Equation(7.11)because p k is a measure of plastic deformation while 5i is a measure of total (elastic plus inelastic) deformation. The data necessary for computing the energy dissipated in plastic hinge rotation is generated at each time step in solving the equations of motion. Since the use of this data was not anticipated at the time the computer program was written, it was not included in the computer output. Therefore, energy calculations involving this source of dissipation have not been made. Figure 18 shows the energy curve for the four-story undamped elastic structure subjected to the Helena ground motion. The curve is plotted as energy per unit mass vs. time. Since the structure is elastic and undamped, the total energy input is equal to the total energy stored in the structure in the form of kinetic plus strain energy, and only one curve appears. Figure 19 shows the energy curve for the four-story viscous damped structure in the same earthquake. In this figure two curves appear, the upper one being the total energy fed into the structure, and the lower one being the total stored energy (strain plus kinetic). The difference between the two curves represents the energy dissipated in viscous damping. The energy curves for the four-story structure with elastic-plastic shear walls and the Helena ground motion is shown in Figure 20. Here again two

branches appear, the upper branch representing the total energy input, and the lower branch representing the stored energy, strain plus kinetic. Figures 21 and 22 are the energy curves for the eight-story undamped and viscous damped structures. For the four-story structure, the total energy input shows a tendency to increase as damping forces are added to the structural system. This is not evident in the eight-story structure. In any case, the dissipation of energy by the same damping forces far offsets any increase in total input. The fact that the peak value of the energy stored in the structure with shear walls is markedly higher than even for the undamped structure is probably due to the fact that the additional stiffness, due to the shear walls, reduced the fundamental period of vibration of the structure. It is known that small changes in the natural period of the system may result in substantial differences in such elements of response as maximum deflection, maximum energy, etc. It should not be inferred that a decrease in the natural period is necessarily accompanied by an increase in maximum stored energy. If the fundamental period were reduced to zero, which would be the case of a perfectly rigid structure, the maximum energy per unit mass would be half the square of the maximum ground velocity, which is in this case 23.4 in 2 sec. This is far below the maximum energy input in any of these five cases, and in fact is below the peak value of stored energy in each case. Although the complete energy functions for the elastic-plastic frame were not evaluated, the strain energy and kinetic energy were computed for isolated points of time during the last few seconds of the computed response to the El Centro ground motion. Energy calculations for these same times were made for the undamped elastic system, also subjected to the El

-105Centro earthquake. The energy stored in the elastic structure turned out to be about three or four times as great as the energy stored in the elasticplastic structure. Consistent with this, the amplitudes of vibration for the elastic structure were found to be considerably greater than for the elastic-plastic structure, even though the latter suffered greater deflections measured from the zero position.

60 40 2n 0 0 1 2 3 4 5 6 7 8 9 10 TIME, seconds. Figure 18. Energy Per Unit Mass, Four-Story Undamped Elastic Frame, Helena Earthquake.

6 VA | |TOTAL INPUT ~S40- _ z w 0 STORED ENERGY (STRAIN PLUS KINETIC) i00 0 1 2 3 4 5 6 7 8 9 10 TIME, seconds. Figure 19. Energy Per Unit Mass, Four-Story Viscous Damped Frame, Helena Earthquake.

90 TOTAL INPUT 80 _ 60 0 u) o 40 TIME seconds Figure 20. Energy Per Unit Mass, Four-Story Undamped Frame with Shear Walls, Helena Earthquake. with Shear Walls, Helena Earthquake.

60 50 r30 Cr w az 20 0 L 10 0I 0 1 2 3 4 5 6 7 8 9 10 TIME, seconds. Figure 21. Energy Per Unit Mass, Eight-Story Undamped Elastic Frame, Helena Earthquake.

60 50 *I~ ~~~~~ / r —-~~~~~~~~~~~~TOTAL INPUT.E40 (n z 30 wr 10, 20 0 0 I. 0 I 2 3 4 5 6 7 8 9 10 TIME, seconds. Figure 22. Energy Per Unit Mass, Eight-Story Viscous Damped Frame, Helena Earthquake.

CHAPTER VIII SUMMARY AND CONCLUSIONS In the foregoing chapters, several methods of solving the equations of motion for dynamic systems have been presented. In Chapter II, there was shown the classical method of normal modes for finding the motion of undamped linear systems. It was shown that the undamped frequencies and modal columns can be used to generate the exact solution for a damped linear system if the damping matrix is a linear combination of the inertia and stiffness matrixes. For this special case, the equations relating the undamped modes and frequencies to the damped system were given. If the damping matrix is linearly independent of the inertia and stiffness matrixes, one can still find a set of eigenvalues and modal columns for the damped system, which then leads to uncoupled equations of motion. In this case, the modal columns and eigenvalues occur as complex conjugates. The real part of each eigenvalue indicates the damping term and the imaginary part indicates the damped natural frequency for that particular mode. In Chapter III, methods were presented for evaluating the eigenvalues and modal columns both for undamped and damped systems; and a method of integrating the uncoupled equations of motion was given. In Chapter IV, a general numerical method was established for solving the differential equations of motion for non-linear systems. The method is well adapted to high-speed digital computation and imposes no restrictions on the non-linearities other than piecewise continuity of the functions. In Chapter V, the equations of motion were set up for a multistory building, treated as a lumped mass system, in which elastic and -111

-112inelastic flexural deformation of the frame members, elastic and inelastic shear deformation of the walls, and viscous damping, were all taken into account. The properties of the structure were, of course, idealized to produce a solvable dynamic system. It was proved that for a rigid-jointed rectangular elastic-plastic bent in lateral motion, the locations of the plastic hinges and the lateral velocities of the floors determine uniquely the plastic strain rates throughout the bent. The response of several variations of two multi-story buildings to two different earthquake ground motions were computed by the methods herein developed, and were summarized in Chapter VI. Some of the results were shown in graphical form to illustrate the character of the response. Other results were shown as a tabulation of the maximum magnitudes of the critical parameters. The energy relations, both for stored energy and for dissipated energy, have been developed in Chapter VII. From the results of the studies presented herein, the following conclusions can be drawn: The methods of analysis, both for linear and non-linear systems, are practical and can be carried out on computing equipment available today. The solutions summarized in Chapter VI and the energy calculations for Chapter VII were carried out with the use of an IBM- type 650 computer, which is a magnetic drum computer with a 2000-word memory. The first mode integrations each required about ten to twenty minutes of machine time. The multi-degree of freedom integrations each required from one to ten hours of machine time, depending on the complexity of the system. Machines a hundred times faster and a hundred times larger than the IBM type 650 exist today.

-113 - The displacement response of a damped linear system is approximated reasonably well by the response of the first mode alone. For each of the analyses reported herein, the maximum displacement in the first mode differed from the total maximum displacement by about ten per cent at most. For the same system, the deviation was greater for the more intense earthquake. However, there is no reason to believe that this would be true in general. Fcr the same earthquake, the deviation was greater for the taller structure. One would expect the modes beyond the first to be more important for a taller structure, simply because there are more of them. In terms of seismic coefficients, which are the parameters of first significance, the first mode approximation is less attractive, espe. cially in the upper stories. For the four-story structure, deviation of first mode from total response in terms of seismic coefficients was eleven per cent at most; but for the eight-story structure, the deviation ranged from practically nil in the bottom three stories to 44 per cent in the top story. Any form of energy dissipation in the system tends to reduce the amplitude of oscillation below that which would prevail if the dissipative force were not there. This is true even for the case in which the energy dissipation occurs through the formation of plastic hinges in the structural members of a frame. The frame becomes less stiff when the plastic hinges are present, but the plastic deformation still results in a decreased amplitude of oscillation. The decrease in amplitude is accompanied by a decrease in the shears and seismic coefficients, as can be seen clearly from Table VI. The results here partially substantiate a thesis proposed by Housner for limit design of structures(20). In brief, Housner's thesis is that if a structure can be built in such a way that inelastic deformation can progress

-114far enough to dissipate a substantial amount of energy, one can then design the structure on the basis of lower seismic coefficients than would have to be considered for a comparable elastic structure, simply because the energy dissipation will act to reduce the magnitude of the maximum shears in the structure. However, it also appears that if the ability of a structure to absorb energy is increased, its appetite for energy may also increase. Energy absorption then may be only partially effective in reducing stored energy, for it may also increase the total energy input. Information presently available is not sufficient to enable one to determine the value of "equivalent viscous damping" to be assigned to any given structure. More information is needed in this regard. Certainly no one can question the desirability of treating a structural system as a viscous damped system, if that can be accomplished. Damping forces in a building are present in both the structural and non-structural elements, and in many buildings the non-structural elements are apt to be the greater contributors. It seems reasonable that the damping forces in the structure, if they are to be considered as viscous damping forces, are probably reasonably well represented by inter-floor damping forces, or symbolically by interfloor dashpots. If a numerical value of force per unit of velocity can be assigned to each of these interfloor dashpots, one can then evaluate the viscous damping matrix. Having this, the eigenvalues and modal columns for the system can be found, giving the frequency and fraction of critical damping for each mode. The author feels that this link connecting the system of damping forces to the fraction of critical damping in each mode is an important one. It seems inherently simpler to assign a value or range of values of force per unit velocity to a structural or non-structural element

-115 - than to assign a fraction of critical damping to a structure as a whole. It is interesting to note that the same system of interfloor damping forces produced 14.8 per cent of critical damping in the first mode for the four-story structure, but only 6.25 per cent of critical damping for the eight-story structure. Since the damping forces are produced in large part by non-structural elements, they might not be significantly different for similar buildings of different height. If this is so, it can be expected that tall buildings are less effectively damped than short buildings of similar type. The actual damping forces in a structure are of course not viscous. For the most part they are attributed to internal friction within the components of the structure and thus might behave more nearly as Coulomb damping forces or some modification thereof. Jacobsen(5) has shown that a single degree system with non-viscous damping, subjected to a sinusoidal input force, can be treated as a viscous damped system without any appreciable loss in accuracy, provided the damping coefficient is chosen so that the energy dissipation in each cycle of oscillation is the same for the viscous damping forces as for the actual damping forces. It seems logical that this idea might apply also to a multi-degree structural system. An additional difficulty then presents itself, for the fraction of critical damping for equal energy dissipation must vary as the amplitude of the motion varies. This is easily demonstrated for the single degree oscillatory system by considering the energy dissipated per cycle of oscillation. The same tendency exhibits itself in the effect of the inelastic deformation of the shear walls for the examples reported in Chapter VI. When elastic-plastic shear walls were added to the undamped system the reduction in maximum response was much more pronounced in the

more intense earthquake. If the energy dissipation by the acticn cf the shear walls were to be represented as equivalent viscous damping, a greater fraction of critical damping would need to be assigned for more intense ground motion. This agrees with the results cf large amplitude vibration tests reported by Hisada and Nakagawa(24), which indicate that the fraction of critical damping as determined experimentally tends to increase with an increase in amplitude. One strong reason for the desirability of reducing a structure to an equivalent viscous damped system is that this puts the structure in a condition to be analyzed by response spectrum techniques, which are applicable only to linear systems.(2). The methods and procedures presented in this report deal with analysis c f structures. Given a specific structure and a specific dynamic force, the problem of analysis is to find the response in terms of deflections, stresses, moments or whatever parameters may be cf interest. Design is an entirely different problem. The designer has no way of foretelling what dynamic loads will actually occur during the life of the structure. Moreover, he does not know the dynamic properties of the structure until his design is complete. The Civil Engineering approach to design for dynamic loads is to emplcy equivalent static loads as far as possible. In highway bridge design, for example, impact is expressed as a static load equal to a fraction of the weight of the vehicle, the magnitude of the fraction being a function of the span length. Similarly, in seismic structural design, earthquake forces are represented as static lateral forces equal to some fraction of the weight of the structure, the fraction itself being some function of the structural characteristics. This approach is both reason

-117able and desirable. Static analysis is much simpler and less costly than dynamic analysis, and it is understood by the profession as a whole. In designing by statics, the designer's intuition serves as a guide to best procedures, and as a rough check on the results. Certainly these attributes are of great value. Dynamic analysis as a Civil Engineering design tool must ordinarily be limited to use as a final check on the design. The primary function of dynamic analysis is to lead to a better understanding of response phenomena, to an improved evaluation of the structural parameters, and to a more complete knowledge of the effects of load parameters and structural parameters on the response. These in turn lead to better design procedures. In this dissertation, the application of the complex eigenvalues and modal columns of a damped linear system to building vibrations, a workable machine method for analyzing the response of a non-linear multidegree-of-freedom system to earthquake, and the uniqueness theorem for plastic strain rates, have been presented for the first time. The author feels that by these features, together with the response and energy calculations and the qualitative observations based on them, this dissertation contributes to a better understanding of structural response to earthquake, and thereby fulfills its original purpose.

APPENDIX A PARTITION THEOREM FOR EUCLIDEAN N-SPACE By Hans Samelson, R.M. Thrall, Oscar Wesler Let En denote the usual real n-dimensional Euclidean vector space of ordered n-tuples. Let 1'''"' in' v1''''n be 2n vectors in n E such that every sequence of vectors a a1..., c, where a is either n i t or q, is a linearly independent set. Let <a,...,c > denote the cone spanned by these a's with non-negative coefficients. The 2n cones in En spanned by the 2n such sequences of a's will be said to partition E if their union is all of En and if the intersection of every pair of distinct cones <al a... cCn>a 1... n> is precisely the lowerdimensional cone ("face", "edge", etc.) spanned by the common vectors, i.e. by those ali for which ai = hi, or {O} if there are no such al1 The most familiar example of a partition is the one with the unit vectors (1,0,0,...,0), (0,1,0,...,O),..., (0,...,O,0,1) in the role of the,'s, and their:negatives in the role of the Ir's. The partitioning cones are then the 2n orthants into which the rectangular coordinate axes divide the space. Although a simple necessary and sufficient condition that the i's and Ir's provide a partition would seem a very natural and desirable thing to have, it appears never to have been written down. An extremely simple characterization does indeed exist, and is given in the following theorem. While several proofs have occurred to us, including one by induction on the dimension n, and none of them excessively difficult (others may occur to the reader) it is interesting that they tend to go to a surprising lot of trouble to settle a matter that seems so "obvious" in the first place. We present the most elementary of our proofs.

-120Let the orientation of an ordered set of vectors a1, c2,...,a be given as usual by the sign of the determinant, sgn det (a,...,). Then any n-l of the a's determine a hyperplane separating En into two halfspaces oppositely oriented with respect to these a's, that is to say, into the two half-spaces lying on opposite sides of the hyperplane. For the sake of simplicity in the statement of the theorem, we assume without loss of generality that the ordered set 1' Y2'''. t is positively n oriented. Theorem: The 2n cones <a,...,a > partition En if and only if sgn det 1 * n (a,'...'a) = (-1)s where s = the number of Tr's among the a's; equivalently, if and only if, with respect to every choice of n-1 a's, the t and j for the omitted index lie separately in the opposite half-spaces. Proof: The equivalence of the determinant and half-space conditions is obvious; we shall give the proof of the theorem using the geometric language of half-spaces wherever possible. Necessity is immediate, for if the t and q of any index, say thej, were to lie in the same half-space relative to any choice of vectors 11... ajl l... a for the remaining indices, then the two cones 1'' a,...,'a> and <,,a*,1,a,...,a > J- j j j+l' n 1 j-l j j+l n clearly have an n-dimensional intersection, contradicting partition. Tc prove sufficiency, consider first the union of the 2n cones, which is a closed set in En. Its complement, if not empty, is an open set whose boundary must be made up of n-l dimensional pieces of certain n-l dimensional cones or faces spanned by certain sets of n-l a's. But no n-l dimensional face can contribute such pieces to the boundary, for the fact that the omitted 5 and r lie in opposite half-spaces clearly

-121plades the relative interior of a face in the interior of the union of the two cones spanned by the face and the omitted 5 and q respectively. Hence the complement is empty, and the union of the 2n cones must be all of En. To show that the intersection property holds, assume the contrary and let r be the largest integer such that two cones exist having exactly r spanning vectors in common but whose intersection contains a point y not in the r-dimensional face common to both. We may, without loss of generality, take the first r vectors of both cones to be the common ones and write them as Q,..., rZ,,.. and <X'' r+'' q > Note that r may be any one of the integers 0,1,2,...,n-2 but may not be n-l, for in that case the cones would necessarily lie in separate half-spaces and then not possibly have in common such a point y as described above. Taking such a y, it may be written in two ways relative to the two cones as follows: (1) Z - al C +... + a ar + ar a +... +a a 1 r r r+l r+l n n bl c1 +... + b i b+ x+... + b n with all coefficients ai,bi non-negative, but at least one of the a r+l'...,a and at least one of the b.*.,b positive. As a matter of n r+l n fact, a,+l...a and b +l...,b must all be positive, for if any one of them, say a, were zero, we could replace ~ by r in (1) thereby contradictJ j j ing the maximality of r. Futher, if any of the a l,..,a, bl,...,b ares zero, then by simply adding to y a point in the common face spanned by a,...~ we can get a new point 7' all of whose coefficients in both 1 r representations are positive. We may therefore assume that y in (1) has all its coefficients positive, noting only in passing that it is then an

-122interior point of both cones, whose intersection is thus n-dimensional, containing a sphere about y. Now consider the unique representation of ~r+1 in terms of the independent set a,...p,,.. r+l 1 r r+l n =x + +x I +... + +... +x r+l 1 1 r r r+l r+l n n We say nothing about the signs of the x in general, but the coefficient i x must certainly be negative, as may be seen at once from the halfr+l space condition, or equivalently from the simple relation det (a,..., Ca, 5,,..., q ) x det (a,..., a,r 1 r r+l r+2 n r+l 1 r r+l r+2 n and the sgn det condition. The point y + t,r+l p which belongs to the first cone for all positive values of t, may be written using its second representation in (1) as (2) 7 + t~ = (b + tx ) a +... + (b + tx ) a + r+l 1 1 r r r (b + tx ) +... + (b + tx ) r+l r+l r+l n n n Taking t to be the smallest positive number which makes one of the coefficients of t'n... in (2) equal to zero, and adding to + t r+ a suitable point from the common face if any of the coefficients of c1'.'.. r in (2) should then be negative, we would once again have a common point 7' with two representations as in (1) above, but this time with one of the coefficients of Tn,..., equal to zero. As this' is in r+l n contradiction to the maximality of r, we conclude tha the intersection property must hold and the theorem is proved.

-123 - It is instructive as well as useful in facilitating applications and for computational purposes to view the foregoing matter in the light of matrices. Accordingly, let A and B be the square matrices whose columns are, respectively, the components of the ~i and ri. We will say that the matrices A and B give a partition whenever their columns do. The partitioning condition on the S's and r's is clearly a centro-affine invariant, that is to say invariant under the full linear group of all non-singular linear transformations on En; under a non-singular linear transformation with matrix P, the matrices A and B are replaced by PA and PB, and the matrices PA and PB give a partition if and only if the matrices A and B do. Further, if sgn det P = 1, then the positive orientation of the first set of vectors is preserved. Under the assumption that matrices A and B give a partition and sgn det A = 1, it follows from our main theorem that sgll det _ = (-l)n, or that sgn det (-B) = 1. Hence sgn det (-B-1) = 1, so that the pair C, -I with C = -B A again gives a partition, with the first set of vectors in a positive orientation. (As a matter of fact, C will be positively oriented even when A is not, as may be verified immediately.) Now let D be obtained from C by replacing some h columns of C by the corresponding columns of -I. Then det D = (-1)h det C' where C' is the principal submatrix of C obtained by deleting the given h columns and corresponding rows. By the main theorem, however, sgn det D = (-1)h, so that sgn det C' = 1. Putting the argument together, we conclude with the simple condition that two matrices A and B, equivalently the matrices C = -B A and -I, give a partition if and only if every principal minor of C is positive.

REFERENCES 1. Frazer, R. A., Duncan, W. J., and Collar, A. R., Elementary Matrices, Cambridge University Press, Cambridge (1955). 2. Foss, Kenneth A., "Ccordinates Which Uncouple the Equations of Motion of Damped Linear Dynamic Systems," Technical Report 25-20, Massachusetts Institute of Technology, March, 1956. 3. Den Hartog, J. P., Mechanical Vibrations, McGraw-Hill Book Co., Inc., New York (1940). 4. Stoker, J. J., Nonlinear Vibrations, Interscience Publishers, Inc., New York (1950). 5. Jacobsen, L. S., "Steady Forced Vibration as Influenced by Damping," Transactions, ASME, 1930, Paper APM-54-15. 6. Schenker, Leo, "The Dynamic Analysis of Tall Structures in the Elastic and Inelastic Ranges," Doctoral Dissertation, University of Michigan, 1954. 7. Newmark, N. M., "Methods of Analysis for Structures Subjected to Dynamic Loading," Report to Directorate of Intelligence, U. S. Air Force, December, 1950. 8. Whitney, C. S., Anderson, B. G., and Salvadori, M. G., "Comprehensive Numerical Method for the Analysis of Earthquake Resistant Structures," Journal of ACI, September, 1951. 9. Criscione, E. S., Hobbs, N. P., and Witmer, E. A., "Outline of a Method of Calculating the Elastic and the Post-Buckling Responses of an F-80 Horizontal Stabilizer to Blast," Technical Note 56-137, Massachusetts Institute of Technology, May, 1956. -125

-126 - 10. Isada, Nelson M., "Design and Analyses of Tall Tapered Reinforced Concrete Chimneys Subjected to Earthquake," Doctoral Dissertation, University of Michigan, 1955. 11.. Cohen, E., Levy, L. S., and Smollen, L. E., "Impulsive Motion of Elasto-Plastic Shear Buildings," Transactions, ASCE, Vol. 122, 1957, Paper No. 2861. 12. Forsyth, A. R., A Treatise on Differential Equations, Macmillan and Co., Ltd., London (1933). 13. Thrall, Robert M., and Tornheim, Leonard, Vector Spaces and Matrices, John Wiley and Sons, Inc., New York (1957). 14. Guest, J., "The Solution of High Order Polynomial Equations by Matrix Iteration," Report SM.236, Aeronautical Research Laboratories, Melbourne, December, 1955. 15. Conte, S., and Reeves, R., "A Kutta Third-Order Procedure for Solving Differential Equations Requiring Minimum Storage," Technical Note No. 1037, Ballistics Research Laboratories, Aberdeen,.Maryland, January, 1956. 16. Ayre, Robert S., "Methods for Calculating the Earthquake Response of'Shear' Buildings," Proceedings of World Conference on Earthquake Engineering, Berkeley, California, June, 1956. 17. Masur, E. F., "On the Fundamental Frequencies of Vibration of Rigid Frames," Proceedings of the First Midwestern Conference on Solid Mechanics, Urbana, Illinois, April, 1953. 18. Hildebrand, F. B., Introduction to Numerical Analysis, McGraw-Hill Book Co., Inc., New York (1956).

-12719. Blume, John A., "Period Determinations and Other Earthquake Studies of a Fifteen-Story Building," Proceedings of World Conference on Earthquake Engineering, Berkely, California, June, 1956. 20. Housner, Gecrge W., "Limit Design of Structures to Resist Earthquakes, " Proceedings of World Conference on Earthquake Engineering, Berkeley, California, June, 1956. 21. Walter, J. K., Benjamin, J. R., and Williams, H. A., "Continued Experimental and Mathematical Studies of the Behavior of Brick Walled Bents Under Static Shear Loading," Technical Report No. 5, Stanford University, August, 1954. 22. Stivers, R. H., Benjamin, J. R., and Williams, H. A., "Stresses and Deflections in Reinforced Concrete Shear Walls Containing Rectangular Openings," Technical Report No. 6, Stanford University, August, 1954. 23. White, Merit P., "Friction in Buildings: Its Magnitude and Its Importance in Limiting Earthquake Stresses," Bulletin, Seismological Society of America, Vol. 31, April, 1941. 24. Hisada, T., and Nakagawa, K., "Vibration Tests on Various Types of Building Structures up to Failure," Proceedings of World Conference on Earthquake Engineering, Berkeley, California, June, 1956. 25. Hudson, D. E., "Response Spectrum Techniques in Engineering Seismology," Proceedings of World Conference on Earthquake Engineering, Berkeley, California, June, 1956.

UNIVERSITY OF MICHIGAN 0Il I11111115 11 llli II 1 67 I7I IILI 3 9015 02223 1677