THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING DYNAMIC ANALYSIS OF ELASTO-PLASTIC STRUCTURES Glen V, Berg Donald A, DaDeppo April, 1959 IP-364

TABLE OF CONTENTS Page LIST OF FIGURES............ i....... i INTRODUCTION o o o o o o o o o o...............................0 THE DIFFERENTIAL EQUATIONS OF MOTION FOR DYNAMIC SYSTEMS......... 3 NUMERICAL SOLUTION OF THE EQUATIONS OF MOTION.....oo........o.... 6 MILNE PREDICTOR-CORRECTOR METHOD o.............o...........oo...o 6 RUNGE-KUTTA METHOD......o..............o......................... 8 THE RESISTANCE FORCES FOR AN ELASTO-PLASTIC BENT..,,,o..o..oo..oo 10 ADAPTATION FOR THE COMPUTER, ooo....o..........oo...oo....oo...... 26 NUMERICAL EXAMPLE..o.................................O...oo...... 28 APPLICATION...................................................... 39 CONCLUSIONS ooooo aoooooaoo............................. a o o o o 4o ACKNOWLEDGMENT............................................... O.. 41 NOTATION................aa....................................... 42 REFERENCES.a..............a..aa............................ o....... 44 ii

LIST OF FIGURES Figure Page 1 Damped Linear Dynamic Systemo.4.oooooooooo.ooo.o.oooo 4 2 Idealized Moment-Rotation Diagram.ooooooo....o..o. 13 3 ElastoPlastic Solution by Superposition.........000 O..00 15 4 Partitioning Vectors o oooO........oooooooooooooooooooooooo 20 5 Bent for Numerical Example o..oo.....ooo....ooo o.oo o oo...6 29 6 Matrix Elements k..and j. o 0o o.............................. 30 7 Matrix Elements bigaand Vk o.oo ooo0 o...................... 31 iii

INTRODUCTION The basic tool of the theory of linear vibrations was forged by Daniel Bernoulli in 17535 when he enunciated the principle of resolution of vibrations into independent modeso 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 infintely many degrees of freedomo The theory of vibrations of nonlinear systems has not fared as well as its linear counterparto Bernoulli's principle of resolution of vibrations into independent modes, which plays such an important role in linear theory, is inapplicable to nonlinear systems. Exact solutions have been obtained only for a few special cases of single degree of freedom systems. For steady-state forced vibrations of single degree of freedom systems with nonlinear damping, methods have been devised for finding approximately equivalent linear systems (1) For vibrations due to transient forces it may sometimes be possible to approximate the system by a linear system over a limited range, or by one of the nonlinear systems for which the solution is knowro If strong nonlinearities are present, this approach is apt to be unsatisfactory; in such cases, for single degree of freedom systems,, the solution can often be obtained by the

use of an electronic analog computero Alternatively, numerical methods can be employed, For a nonlinear system having many degrees of freedom, the complex circuitry required for tuhe analog computer becomes prohibitive, and numerical methods provide the only suitable approacho A high-speed digital computer is essential. This paper presents two numerical methods for solving the differential equations of motion for nonlinear systemso The differential equations of motion are formulated for a lumped mass multi-story elasto-plastic framework subjected to dynamic lateral forces, taking into account elastoplastic deformation and viscous damping. Uniqueness of the solution for an elasto-plastic frame is provedo

-3 - THE DIFFERENTIAL EQUATIONS OF MOTION FOR 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 1o 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, One then has a linear dynamic system of n degrees of freedom. Let mi = mass of the ith body, xi = displacement of the ith body, relative to a coordinate system fixed in the base, y(t) = displacement of the base, relative to a "fixed" frame of reference, Fi(t) = external force applied to the ith body, and let dots denote differentiation with respect to time, Define the stiffness coefficient kij as the force exerted on the th2a, i body by the springs when the configuration of the system is xj =1 unit length; xk = O, k f j; and xk = 0 for all ko Also, define the damping coefficient c.. as the force exerted on the ith body by the dampers when the ij

-4 - t- x1 ^X2 — y (t) xl and x2 are displacements relative to base y is displacement of base Figure 1. Damped Linear Dynamic System.

configuration of the system is x.j = unit velocity; xk = 0 k ~ j; and xk = 0 for all ko 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: n n ijr ) ci + k kjxj Lji -j =l The Newtonian equation of motion for the ith body is then n n mi (~i + jY) + + j kijxj = Fi(t) (1) j1 j=1 Thus, for the entire system, m i ciji + kijxj Fi(t) - miY(t) (i = 1,2Pn) (2) j=I j =l Let fi(t) Fi(t) - miY(t) (i 1 i 2,o.n) (3) Equations (2) then become n n Mix + c_ c Xj + j k jX = fi(t) (i -=,2,.n)) (4) j=l j=l which are the general differential equations of motion for a damped linear dynamic system with n degrees of freedom. In Equation (4) the term cijx. will be called the damping n L 1J J force, and the term kijxJ will be'called the restoring force. The sum of the damping and restoring forces will be called the resistance force Rio

In a general nonlinear system the resistance force Ri for each mass may be a nonlinear function of time and of the velocities and displacements of all masses. The equations of motion for a nonlinear system are then mixi + Ri (xl1 x2 ' ~ Xn, Xl, x2,~~ Xn' t) = fi(t) (i = 1,2on) (5) NUMERICAL SOLUTION OF THE EQUATIONS OF MOTION There are many methods of approximating the solution to Equations (5) numerically, While these differ in details they all employ the same basic idea. Let h be a small finite increment of time. One takes the known discrete values of xi, xi, and xi (i = l,2,,on) at times t, t-h, t-2h, etc., as required by the particular numerical process, and the values of f.i which are explicit functions of time, and projects these to evaluate xi, xi, and xi at time t+ho In this manner, one advances step by step through the solution. No one method of numerical integration can be claimed to be the best for all problems. The following two methods have exhibited excellent behavior in machine computation of structural response to dynamic loads. MILNE PREDICTOR-CORRECTOR METHOD The Milne Predictor-Corrector method(2) is well suited to problems in which (a) damping forces are absent or negligible, and therefore the resistance functions Ri do not involve the velocities x; and (b) the driving forces fi can be approximated by smooth curves through sets of discrete points equally spaced in time. Blast problems, for example, often meet these requirements~ To employ the Milne Predictor-Corrector method one

must know the discrete values of xi and xi at the ends of four consecutive time intervals of equal length h. The steps of the method are: Predictor xi (t+h) x(t) + xi(t-2h) - xi(t-3h) + 5i(t) + 2xi(t-h) + 5xi(t-2h) (6) Corrector 2 xi(t+h) 2xi(t ) xi(t-h) + i h (t+h) + l2xi) (t) + xi (t-h)] (7) By recursive use of these two formulas one advances step by step through the solution, The Milne Predictor-Corrector process does not require calculation of velocitieso Also, being noniterative, it is not accompanied by a conver17 6 vi gence problem. The dominant error term in the predictor is + 2T h x, 1 6 vi while in the corrector the dominant error term is - 2 h xi The truncation error committed in each step is therefore approximately 1I |xi - xi|. One can easily program the computer to accumulate a set of error functions -il 1 Xi l (8) and halt when any ci exceeds a present maximum, The error functions are not absolute error bounds, for they consider only the dominant term in the truncation error in each step. Nevertheless, they are worthwhile checks to incorporate in the program. Because each step of the Milne Predictor-Corrector process uses information from three preceding time steps, some "starter" method must be employed to calculate the initial stepso Also, the length of the time interval h must remain constant. Any time after the sixth cycle the time interval can be doubled in an obvious manner. However, reducing the time interval would require interpolating to find the needed values of displacement at intermediate points of the preceding intervals.

RUNGE-KUTTA METHOD For problems which involve very irregular driving forces, such as the earthquake response problem, it is advantageous to use a singlestep 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 procedure(3) is a single-step method well suited to high-speed computers. For the equations Yi Y yi(:Y1 Y2o- yn (i = l,2,.,n) (9) the formulas for the Runge-Kutta third order procedure are Km ~ hYi(YL'Y1,y2 y1,9t) Kil = hjY(y1 + hK + P10 Y2 PK20. y Yn + PKnO, t + ph) (10) i2 =1 hYi(y] + -r] 10 + rKlly n + [q-r] KnO + rnl t + qh) and Yi(t+h) = Yi(t) + Kilo + mKil + ni2 (11) The relations governing the constants ~, m, n, p, q, and r are derived by expanding both sides of Equation (11) in power series in h, and equating coefficients of the powers up through h35 The relations are Q +m+n =1 mp + nq = 1/2 2 2 (12) mp + nq =1/3 npr 1/6

to which there are infinitely many solutions. One of the more useful solutions is = 1/4 m = 0 n = 3/4 p =1/3 q = 2/3 (13) r = 2/3 This makes m = (q-r) = 0 and thus reduces the number of arithmetic opera(4) tions required. A second solution, due to Conte and Reeves,( is = o 62653829327 m 85614352807 n = -48268182134 (14) p 62635829327 q =.07542588774 r = -o55111240553 In this solution, = p = (q-r), which leads to a reduced storage require ment for the computero To adapt Equations (5) them in the form to the Runge-Kutta process, one can write i - z1 z i f i(t ) - Ri (xl, X2,oXn, l, z2.o* Zn) t) (1 5) The solution is then straightforward.

The numerical methods above are applicable to any system, subject only to the restrictions that the number of degrees of freedom must be finite and the resistance functions R. must be single-valued and pieces wise continuous. THE RESISTANCE FORCES FOR AN ELASTO-PLASTIC BENT The specific problem considered in the remainder of this paper is the response of a multi-story bent to lateral dynamic forces. The bent is a rectangular plane framework of elasto-plastic members, loaded by lateral dynamic forces applied at the joints and in the plane of the bent, Connections may be either pinned or fully restrained. The mass of the structure is assumed to be concentrated at the joints. Damping is assumed to be vise cous, Deformations due to axial forces and shear forces in the members are neglected, and the effect of axial forces upon the stiffnesses and plastic hinge moments of the members is neglected. These assumptions yield a dynamic system having degrees of freedom equal to the number of stories in the bent. The principles employed can be extended to more general problems, but such extension is not considered herein. Consistent with the assumption of viscous damping, the resistance functions R. are separable into damping forces Di, which are linear functions of the velocities, and restoring forces Qio As in Equations (4), the damping forces are n Di -= c (16) L e j=l The coefficient cij is the damping force that would exist at the ith floor

-11 - if the structure were passing through the equilibrium configuration with unit velocity at the j-th floor and zero velocity at all other floors. The restoring forces for an elastic frame can be expressed in terms of a matrix of stiffness coefficients k such that n Qi = kjx, (17) j=l just as in Equation (4). The forces Q (i = 1,2,.on) are the forces necessary to hold the frame in static equilibrium in the configuration defined by the lateral deflections xj (j = 1,2,n)o The j-th column of the matrix k 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 eccentricity of dead load (the "overturning" effect) can be taken into account in evaluating the matrix if desiredo Let the subscripts k and & denote locations in the frame, such that each value of the subscript denotes the end of a specific structural member at a specific joint; and let Mk be the moment at location k, considered positive when it tends to rotate the end of the member clockwise. Then the end moments in the members of an elastic frame can be expressed in terms of a matrix of influence coefficients p, such that n M'k =j I Jkxj a (18) j = The j-th column of the matrix L is the set of end moments that would exist in the members if the frame had unit lateral deflection at the j-th floor

-12" and zero lateral deflection at all other floorsO The matrix | can be evaluated by conventional static analysiso If 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 restoring forces and moments. For this paper it is assumed that the connections are either pinned or fully restrained, and that all members in the frame, columns and girders alike, have ideal elasto-plastic moment-rotation characteristics. The typical moment-rotation curve follows a linear elastic branch until the moment reaches the plastic hinge moment, Upon further deformation, the moment remains constant and a plastic hinge forms, making a "kink" in the member. Upon reversal of strain after reaching plasticity, the momentrotation curve follows a path parallel to the original elastic brancho 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 loopo This idealized behavior is illustrated in Figure 2o 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 disappearO Cohen, Levy and Smollen developed a procedure for adapting the method of normal modes to an elasto-plastic frame with infinitely rigid (5) girders, and Schenker indicated a method of extending this to include

-13 - I 0\ Figure 2. Idealized Moment - Rotation Diagram

(6) the frame with flexible girderso ) 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 accordinglyo In this way a continuous solution can be foundo For a step by step numerical solution, an analogous and somewhat simpler 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 stepo 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 plastic hinge rotations did not occur. The resulting moments may exceed the plastic hinge moments at one or more locations in the frame. To remedy this, "corrector~ solutions are superimposed in which the frame has hinges at these locations, and these hinges are rotated in such a way that the total moments, obtained from superposition of the elastic and corrector solutions, nowhere exceed the plastic hinge momentso The superposition of elastic and corrector solutions is illustrated in Figure 35 The effects of plastic hinge rotations upon the forces at the floors and the bending moments in the frame members can be expressed in terms of two matrices of influence coefficients b and v, evaluated as described belowo

-15 - Elastic I 1 7 Corrector Elasto-Plastic Figure 3. Elasto-Plastic Solution by Superposition.

-16 - Let cp be the plastic hinge rotation at location J, taken as positive when it corresponds to clockwise rotation of the end of the membero (The hinge rotation in Figure 3 is positive,) Take the bent in its initial (unstrained) position, insert a hinge at location a, and rotate the hinge through an angle cpq = 1, with zero lateral deflection at all joints. Evaluate the lateral forces at all floors and the end moments for the members in the frame, The lateral force at the i-th floor is the influence coefficient bi, and the moment at location k is the influence coefficient vkO Both matrices b and v can be evaluated by conventional static analysis. It will be seen that only those locations at which plastic hinges form are significant to the dynamic responseo If the lateral deflections x and the plastic hinge rotations cp were known. for some instant of time, one could evaluate the restoring forces and moments as Qi -- kijxj + bi (19) J S and 'Mk Lkjxj + Vk ) (20) j g The deflections x can be found by integrating the equations of motion., but one must also find the plastic hinge rotations which satisfy all the constraints of the elasto-plastic system, An iterative method of solution is developed below, Consider the time rates of change of the forces Qi, the moments Mk, and the plastic hinge rotations qko Let the plastic hinge moment at location k be Pko The elasto-plastic constraints are:

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 decreasingo If the moment is constant, 'the hinge may rotate in a direction. consistent with the momento If the moment is decreasing, the hinge cannot rotate, 3o Equivalent conditions exist for a negative plastic hinge o 4o If the moment is less than the plastic hinge moment in magnitude, the hinge cannot rotateo These conditions may be expressed mathematically as the following set of constraints: l ~PkS Mk < P 20 If Mk= Pk the n M < 0, O. - 0, and Mq = 0 (21) 5. If M~ "Pk then Mk > 0, cpk > 0 and Mkc =0 i 40 If ~Pk < Mk < Pk then = 0 Differentiating Equation (20), one gets Mk = 'kjXj + vkP (22) The following uniqueness theorem can be established: Given any rigld-jointed rectangular plane framework of ideal elasto-plastic flexural members, and given any xIs whatever, and any MIs not exceeding the plastic hinge moments* there exists one and only one set of MNs and cpls which satisfies Equation (22) and Constraints (21)o The proof followso * The case in which all members at a joint become hinges simultaneously is excluded,

-18 - First, observe that if, for one or more values of k, -Pk < Mk < Pk then k = O and there is no constraint on Mk. Thus the solution for cp 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 ak = (Sgn Mk) k jXj (23) ki= (Sgn Mk)(Sgn M) vk 2) kR Zg = (Sgn M )(P" Equation (22) then can be written Z 6-k5k y = ak, (214) a kz (24) where 6k is the Kronecker delta, 5k = 1, H k = R 5 i o0, k, Constraints (21) become y ~k Zk > ~0 (25) Ykzk = 0O The vector a is unrestrictedo To establish uniqueness, it must be shown that the matrix ~ has properties which guarantee that, corresponding to any vector a, there exists one and only one pair of vectors y and z that satisfies Equation (24) and Constraints (25)o This amounts to a partition problem in n-dimensional Euclidean space.

Consider the columns of the matrices -6 and t as vectors l, n 2' n; 5l' 2','oo'oan That is, let {Ti} = {-ji} and {i } {=ji Then let al, a2,~~ on be a set of vectors such that every ai is either Ti or if, There are 2n such sets, If the t's and Ids partition the space, every vector a in the space can be formed by linear combination of one and only one set of ai's with nonnegative coefficients,* In other words, for any given vector a there is one and only one set of Qi's for which {a} = j qi {ai}, (26) i where every qi is nonnegative. If this is so, then for those i for which i T = ii, one has yi = qi and zi = 0; and for those i for which a. = i, one has yi = 0 and zi = qio Thus the vectors y and z corresponding to any given vector a are unique and satisfy the constraintso To illustrate this concept, consider the two-dimensional case. Suppose positive plastic hinges exist at locations 1 and 2, and all other locations are elastic. Further, suppose that the matrix v is (see example, page 32), r 69,700 -16,580. k = 169580 -10,5860 in-kips/radian k L -16,580 105560] In this case, Vkp = bkn and Equation (24) becomes r 69,700 L16,580- {} + ~{y} = {a} L -16,580 105,56o0 L0 1 The columns of the matrices in this equation are plotted as the vectors tl, A2, ql and '2 in Figure 4. The four vectors partition the plane so that * If any coefficient is zero, say m = 0, obviously one can have either am = qm or = ~mo However, in this case ym = Zm = 0 and uniqueness is preserved.

-20 - 2 Figure 4. Partitioning Vectors.

-21 - any vector {a} in the entire plane can be formed by linear combination of one and only one set of ai's with nonnegative coefficients. (7) By virtue of a recent theorem due to Samelson, Thrall and Wesler, necessary and sufficient conditions for the columns of | and -5 to partition the space are that every principal minor of the matrix t be positiveo 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, aboveo Let the hinged joints in this frame undergo hinge rotations cp, with lateral deflection of the floors preventedo The resulting strain energy in the frame is E - P k o (27) k The moments are k- \ v k o (28) Combining Equations (27) and (28), one gets E = I Z I ke' (29) k Because Equation (29) is a strain energy equation, it is a positive definite quadratic form. The principal minors of v are therefore positiveo Finally, since t was obtained from. v simply by changing the signs of the off-diagonal elements of certain rows and the corresponding columns (see Equations (23)), the principal minors of ~ are positive0 It follows from the Samelson-ThrallWesler theorem, that a unique solution exists, The solution can be found by the following iterative scheme0 Let zP be the p-th approximationv to the vector z, with all elements of zP nonnegativeO (Zero is a convenient and satisfactory first approximation,)

-22 - For convenience in notation, define kO = kkn+l = z0 = Zn+l = 0 Starting with k = 1, let k-l n+l P+l- V a P+l V 30) Uk =ak- L ~kg z (3~0) 2=0 =,=k+1 Then, p+l if p+l p+l k p+1 p+l and if uk < 0 let zk =0 Repeat this for k = 2, k = 3,o.k = n. Then let n k+1 - C ikZi ak (k 1,2,. n). (32) 1=l The (p+l)st approximation to the vectors y and z is then complete, It satisfies Equation (24) and the constraint zk > 0, but may fail to satisfy the other constraints in (25). The process is repeated until the amount by which the approximation fails to satisfy all the constraints is insignificanto Once the vector z is determined (and therefore the p's) the rates of change of the restoring forces are found by the equation Q! ' ^ ^ j k ^ b~ (33) j 2 The foregoing procedure is, in essence, an adaptation of the GaussSeidel iterative method of solving linear equations Equation. (24) and Constraints (25) reduce finally to a system of linear equations, since at least

-23 - one of zk, Yk must be zero for each value of ko The nonzero 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 rapido Using this procedure, at each step one must perform three matrix multiplications and solve a system of equations equal in number to the number of plastic hinges, An obvious alternate is moment distribution, whereby at each step one must solve a set of equations equal in number to twice the number of members in the structureo In developing the above procedure, rates of change of deflection, hinge rotation, and moment have been used. As soon as one introduces finite time intervals in place of differentials, one loses the claim to uniqueness, since it cannot be claimed that reversals in strain do not occur within the time interval. However, by taking the intervals sufficiently small and considering the strain rates as 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 system0 To set up the finite difference system, let the subscripts s and s+l denote a condition at the end of the sth and s+lt time step; eog,, Xi s+l is the value of the variable oi at the end of the s+lst time stepo Also, let Axi = X s+l " xis s AQ = Qi,s+l Qi,s ' AM = Mk s+l Mk, s an I~ = ~, ~

Then the finite difference equations are AQi - j k.Ax. + b, (55): =j and AMk L= 1k iX I VkApg, (36) j Q which correspond to Equations (33) and (22)o The elasto-plastic constraints are: 1. If 0 < Mks -< Pk then either AMk =Pk Mk s and Acp < 0 (37) or -Pk Mk s < Mk < k - Mks and Ack 0 and 2o If P < Mk < then either AMk = Pk Mk and APk 0 (38) or Pk Mk, < AMk < Pk Mk s and Ack = In the finite difference system, by taking Axj sufficiently large, one can construct mathematically a situation for which no solution existso For example, if, for some value of k, Mk s were negative, one could choose the Axj so that j tkjAj would be a very large positive quantity, so large that a solution of Equation (36) would require either that AMk exceed Pk - Mks or that Acp be negative, both of which are in violation of Constraints (38) for negative Mk s The occurrence of such a situation in the course u 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 significance of the numerical resultsO

-25 -Should such a situation arise, it would simply mean that the time interval was taken too large in the first placeo

-26 - ADAPTATION FOR THE COMPUTER In using the computer, both storage requirements and computation time can be reduced by expressing the moments in nondimensional form, To do this, one can define Mk k = Pk VkkK k ^kkj k2 V 0k(39) Xk j Pk notation, Equations (35) and (36) become AQi = J ki jAxj+ ie (4o) j ~ and the elasto-plastic constraints (37) and (38) become o1 If 0 < p < 1, then - ks k either Ak = (1 - Pk s) and A_ <K 0 (42) or (-1 - Pk s) < APk < (1- Pk s) and A = 0 2. If -1 < pk < 0 then either Apk = (-1 - Pk and Ak > 0 (43) or (-1 - Pk ) < p < (1 - Ps) and Ak = 0. kys- k ~ k~s Fk

-27 - The iterative procedure of Equations (30), to this formulation of the problem. Let r = number of moment locations (31) and (32) is readily adapted (the order of the matrix a) and A = p-th approximation to LAk k Define 0 k kO o = 0 for all = a 0 kr+l = mAP = 0 r+l Then, starting with k = 1 uP = 3 X Ax Pk iZ Ik J j If (Pks + UkP) > 1 if -1 (Pks + Uk) 1 if (Pk s + ukP) < -1, Repeat this for k = 2, k = tiono Then let k for all k, for all p and p = 1, compute k-1 r+l + I ~a Ap + ak2 A (44) kg O kg Y 2=0 =k+l let APp = (1- ks ukP); let Ar = 0; and (45) k let Ap = (-1- Pk,s - ukP) = 3,..ok = ro This completes the p-th approxima r e - 1\fk -tk I k=l (46) which is an indicator of the amount by which the solution fails to satisfy the constraints (42) and (43). If eP is larger than some predetermined error test 6, one returns to k=l and computes the p+l st approximation. When the convergence test is satisfied, one accepts the last approximation as the correct solution for the Ark and computes AQi and Apk from Equations (40) and (41

-28 - NUMERICAL EXAMPLE The following numerical example illustrates the foregoing processes. Consider the undamped bent shown in Figure 5. Taking fy = 36,000 psi, one gets for plastic hinge moments P =P2 P3 = 1398 in-kips P4 = 2072 in-kips P5 a 1583 in-kips Because the bent is symmetric and the effects of axial forces on the member stiffnesses are ignored, the stress distribution will also be symmetric, Hence only half the moments need be recorded. Moreover, the moment at location 5 cannot reach plasticity and therefore does not need to be recorded. The matrices of influence coefficients k.i, bi;/ -kj, and vk which describe the elasto-plastic properties of the system are the forces and moments shown in Figures 6 and 7. The effect of eccentricity of the dead load on the moments may be taken into account if desired. It is neglected in this example. The complete matrices are r 55 73 -24.11] l -1 ki = kip in-1 J L-2411 20.20 r 1425 2019 -1674 -345 -b. kips/radian = L-1325 -1584 563 1021o 726.1 -662.6 1009.6 -791.8 tkj = -83569 281.3 in-kips/inch -172,7 510.5

-29 - W = 20 kips 2 m2 =.05181 kip sec2 in-1 W1 = 30 kips m1 =.07772 kip sec2 in'l 24 --- — 424' - Figure 5. Bent for Numerical Example.

FrX( vm"X gulsn*oE3 Xy;-sqW *9 amnT2?Z, K-ZX 1-i1 tIrl T= x

-31 - b21 b2 2~^1'1 1 14v1 b, /v21 b1 7 7 o of b 22 ofII Figure 7. Matrix Elements bil and vkl

69,700 25,720 Vk-, 514 16, 580 25,720 88, 300 -31,370 -56,930 -9,140 -31,370 80,010 -48,640 -16,580 -56,950 -48,640 105,560 in-kips/radian Converting the moments to nondimensional form according to Equations (39), one obtains 29. 12 26, 58.5194.7222 Xkj -.5987 -.0833 1.0000,3690 k o1311 -.16o5 31.97 -25.08 -.4740 -.5664.2012.2464 -29.25 9.83 -6,78 kip 20.04.-1 in -o2328 -.7993 -.6829 1,0000.2913 1.0000 -3553 -.4350 - 1142 - 3921 1.0000 -.4102 Suppose that at time t0 the system is in the configuration xi = 2.600 in x = 10,00 in/sec. x2 = 4.900 in x2 25.00 in/sec. and that plasticity has not yet been encountered at any location, Suppose also that the driving forces at time t are fl = f2 = 10,00 kips,

and that these forces decay linearly to zero at time tO + 0ol second. The complete process for determing inelastic behavior and executing one time step of Runge-Kutta integration is given belowo The Runge-Kutta constants for this example are those given in Equation (13), namely, =, 1/4 p = 1/3 m = 0 q = 2/3 n = 3/4 r = 2/3 and a time interval of h = 0.02 seco is usedo These initial conditions and time interval have been chosen for illustrative purposes, and are not necessarily typical or optimum. Because plasticity has not yet been reached, the restoring forces at time t0 are given by Equation (17), Q kx r 35 73 -24.11 2.600] Qi = kijxj - L -24,ll 20.20 L4.90o0 or Q1 = a25y24 kips Q2 = 36.29 kips and the nondimensional moments at time t are.5194 -.4740 - =\ ix.7222 -.5664 2.60oo Pk = / -jx. -= 4 L Jk J -.5987 2012 L4.900J -.0833.2464

or p1 P2 p3 P4 From Equations (15) an( xl (t0) x2 (t 0) = Z' (to) '2 (to) = The first of Equations x O K = h x 10 1 x z 10 h z1 Z h 20 2 =:. 9722 = - 08976 - -5707 9908 d the initial conditions, 10O00 in/sec 25.00 in/sec 2 (10 + 25.24) = +45 ~07772 1 (10 - 36.29) = -50~.o5181 (10) yields (t) = (.02)(10.00) = 0 354 in/sec2 7,4 in/sec2 0 200 in 0.500.in 9.068 in/sec -10o148 in/sec, (to) = (,02)(25.00) (to) = (.o2)(453.4) (t ) = (.02)(-507.4 0 The superscripts x and z identify the variable associated with Ko Before proceeding to the second of Equations (10), one must evaluate the functions Q(Xi + - Kx ), which in turn require the evaluation of the AO's. The first terms in Equations (44) are z j 1 3 Xkj=j j Xkj (3 K jo) j 05194.7222 -.5987 -.0833 -,4740 -. 5664.2012.2464.200 s5oo -,0444 -.462 - 0064 o0355

-35 - The iteration of Equations (44) and (45) yields, for p = 1, u1 =-.0444 p + u = -.9722 0444= -.o0166 Ar1 +.o0166 U2 =- o462 + (3690)(.o0166) = -.o4oi 1 p2 +u = -.8976 - o401 = -.9377 A1 o I42 0 1 3 P3 + A1 Al! = 1 U4 = P4 + Al4 = From Equation (46), 1 e = -00o64 + (-.1311)(o0166) = -.oo86 u3 = -.5707 oo0086 = -.5793 0 +.0355 + (-.1605)(o0166) = +,0328 u4 =.9908 +.0328 = 1.0236 -.0236 and the process must 2 2. 2 bE j |At - AL I = o0402, i e repeated. The second cycle t =.0111 0- gives AL3 = 0 Al!4 = -.0245 and e = o0064 As the process is repeated, e approaches zero and the iteration converges to

-36 - mA2 At3 A4 (40) AQs = o.1009 - 0 0 -= 0246 then yields - k 3 J i J P1 + Ap1 P2 + 42 P3 + Ap3 P4 + AP4 = -10000 = 9201 = - 5617 1.0000 Equation Kj ) + Y Pi pte or AQ2 LQ2 = -1.636 = o759 + o.484 = -1.15 kips - 0.783 = 0.98 kips and 1 + AQ1 Q2 + AQ2 = -26.39 = 37.27 kips kips One can now return to the second of Equations (10) and compute i h Lzi + i 3 io and Kz h f i (to + ) Qi (xj + )] The results are Kx= (.02)(10.00 + o068) 0.260 in 11 ) x a10148 = (.02)(25,00 - 18) = 0o432 in 21 3 zK.02 (10.002 - x 2.00) + 26.39 = 9.193 in/sec 21= (o0517) (looo x 27 = 103784 in/sec,021 C21,i8o) \jo.a -- x 2.00) - 37.27 = 10784 in/seco 21 o5181 3~~~~

-37 - Before going to the last of Equations (10) one must evaluate Q(xi + 2 ) ). Using Equations (40), (41), (44), (45), and (46) again, 3 il just as illustrated earlier, one gets i1 = oo08o pL_ + AP1 =.000o A2 =-0 P2 e+ ~2 =.8958 A3 - o0 P3 + 03 =-.5861 A44 -.o460o P4 + Ap4 = 1o0000 1Q2 b42 = -021 kips = 0 50 kips Qi + AQ1 Q2 + AQ2 = -25.45 3 56.79 kips kips The last of Equations x 12. x r22 z 12 z 22 (10) is now executed. = (o02)(iooo + 2 x 9.193)= 3 0.323 in (.02)(25.00 - x 10.784) = 0o356 in.07772 )(10.00 x +.02) 2 - (.077) (10o00 L x 2000) 60.79 - 10.856 in/sec. = ( o5, 81 ) [(lOoOO 3 x200 05 ] Equation (11) now yields the displacements and velocities for the end of the time step. x1 (t0+)= XI (t) + Xi + 4 2 = 20892 in 10410 ~ r lo 4 "12 X2 (to+h) = 5o292 in xl (t0+h) = Zl(O + I Kl0 + 4 K12 = 18.85 in/sec x,(t +h) = 14.52 in/see

-38 - The terminal values of Q and puted from Equations (40), (41), (44), are Ag!1 Ag2 mA3 Ag!4 0 0 0 - o0630 p for this time step must now be com(45), and (46), as before. The results pi (t +h) - -.9916 P2 (t0+h) = -.8584 p (t0+h) = -.6236 P4 (t0+h) = 1.0000 kips kips Q1 (to0h) QS (t0-+h) Q 0(Q~. = -23,83 = 35.91 One now advances to the next as the initial values for the next stey time step, taking these terminal values

-39 - APPLICATION The multi-degree elasto-plastic response problem has been programmed for solution on an IBM Type 650 computer at the University of Michigan. Both the Milne Predictor-Corrector Method and the Runge-Kutta Method have been used, The programs were written for the basic IBM 650, without index registers, automatic address modification, or floating point hardware, If the programs were to be rewritten today, they could be simplified a great deal by taking advantage of these improvementso The 29000-word memory was able to accommodate a system with as many as 16 degrees of freedom and 16 plastic hinges using the Milne method; and with the Runge-Kutta method the system could go to 12 degrees of freedom and 16 plastic hinges. The more severe limitation for Runge-Kutta occurred partly because viscous damping and a form of inelastic deformation not considered in this paper were taken into account in that programo The four matrices, kij, ip Xikj and akg, were stored in the memory initially as part of the data for the main programO A separate program was written to evaluate all of these matriceso Input for the matrix evaluation program comprised the number of stories, number of bays, story weights, member stiffnesses, plastic hinge locations, and plastic hinge moments. Acting on this information, the machine wrote the slope deflection equations, solved them by iteration, evaluated the matrix elements, and printed the matrices on load cards ready for use with the main program. For the matrices Pi^ XkJ and ki, one need consider only those locations at which plastic hinges form during the responseo There is no need to calculate or record the moments at locations that remain elastico

Thus the limitation of 16 plastic hinges is not as severe as it might seem, Intuition and experience with smaller systems may enable one to predict which locations are most likely to become plastic during the response. The unlikely locations can be omitted from the matrices. One can process the output data after the run is complete to check whether or not all of the unrecorded locations did, in fact, remain elastic. The programs were applied to the problems of structural response to blast and earthquake. In each case the driving forces were put in punched cards which were read by the computer as the solution progressed. For the blast problem, the equations governing the dynamic pressures were given to the machine in a separate program and the machine calculated the driving forces and converted them to punched cards ready for unput for the main program. For the earthquake problem, recorded accelerograms were approximated by piece-wise linear functions, and the time-acceleration coordinates of the intersection points of successive line segments were put in punched card form, The machine accepted this input data and interpolated where necessary to find the desired driving forceso CONCLUSIONS Experience with both the blast and earthquake problems shows that the methods of analysis presented herein are practical and can be carried out on equipment that exists todayo Tomorrow's equipment, or even a late version of todays equipment, greatly facilitates programming and permits many of the restrictions to be relaxed.

ACKNOWLEDGMENT This paper is based in part upon research conducted for the University of Michigan Research Institute under contract with Air Force Special Weapons Center, Albuquerque, New Mexico, and in part upon research conducted by the senior author for his doctoral dissertation at the University of Michigan under a doctoral committee of which Dr. Bruce Go Johnston was chairmano Computation was performed at the Statistical and Computation Laboratory of the University of Michigano The authors wish especially to thank Mr, Bruce Arden and Mrs. Sarah Brando, of the Computation Laboratory staff, for their unlimited cooperation.

NOTATION Symbols are defined where they first appear in the texto Those used frequently are summarized here for reference. ij subscripts associated with mass or force k,2 subscripts associated with bending moment or plastic hinge h time interval in numerical integration,,m,n,p,q,r constants in Runge-Kutta method of integration fi driving force mi mass xi displacement D damping force Mk bending moment Pk plastic hinge moment Q. restoring force R. resistance force 1 cpk plastic hinge rotation Pk dimensionless bending moment plastic hinge rotation (with dimensionless bending moment) bi rotation-force coefficient cij o damping coefficient k. stiffness coefficient i0 BPi, rotation-force coefficient (with dimensionless moments) 6bk Kronecker delta Xkji displacement - moment coefficient (dimensionless moments) Akj

-43 -u, displacement - moment coefficient k.j vk rotation - moment coefficient akg rotation - moment coefficient (dimensionless moments) -kRA

REFERENCES 1o Jacobsen, L. SO, and Ayre, Ro S, Engineering Vibrations, McGrawHill Book Co. Inc,, New York (1958) 226-232. 2. Milne, W. E, Numerical Calculus, Princeton University Press, Princeton, 1949. 3. Levy, H., and Baggott, E. A. Numerical Solutions of Differential Equations, Dover Publications Inc,, New York, 1950. 4. Conte, So, 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. 5. Cohen, E,, Levy, Le S., and Smollen, L. E. "Impulsive Motion of Elasto-Plastic Shear Buildings," Transactions, ASCE, 122, 1957, Paper No, 2861. 6. Schenker, Leo, The Dynamic Analysis of Tall Structures in the Elastic and Inelastic Ranges, Doctoral Dissertation, University of Michigan, 1954. 7. Samelson, H,, Thrall, R. Mo, and Wesler, 0. "A Partition Theorem for Euclidean n-Space," Proceedings, AMS, 9, No. 5, (October 1958) 805-807.

UNIVERSITY OF MICHIGAN 3 9015 02223 1750