APPLICATION OF A DISCRETE ELEMENT MODEL TO THE STUDY OF THE STATIC AND DYNAMIC STABILITY OF BARS, ARCHES, AND RINGS Thesis for the Degree of Ph. D‘ MICHIGAN STATE UNIVERSITY PHILIP C. RYMERS .1968 ’17 0-169 Date This is to certify that the thesis entitled APPLICATION OF A DISCRETE ELEMENT MODEL TO THE STUDY OF THE STATIC AND DVNANIC STABILITY OF BARS, ARCHES, AND RINGS presented by Philip C. Rymers has been accepted towards fulfillment of the requirements for Ph.D. Mechanic 5 degree in Elba... a 3 Major professor August 8, 1968 \ . ImmmflA‘L. .1 ‘- I ‘ ‘Lvtfi . L 3.533.; R y f K v-1 - 5 . I Mudu‘éan State 5.., Umvez sity .-__. -"!"-I “or ¥' 1 ABSTRACT APPLICATION OF A DISCRETE ELEMENT MODEL TO THE STUDY OF THE STATIC AND DYNAMIC STABILITY OF BARS, ARCHES, AND RINGS By Philip C. Rymers A physical model of a structural component consist- ing of a series of mass points connected by bars is used to investigate stability of the structural component. The bars are incapable of bending but exhibit axial strain capability. Bending is accounted for by considering each mass point to be a nodal point and by attaching a torsion spring across the node. Equations of motion are written for each mass point, employing a right-hand cartesian coordinate system. Provi- sion is made for non-uniform elastic characteristics in bending and axial deformation. Unsymmetrical force applica- tion, physical properties,and geometry capabilities are also possible with this model. Solutions of the equations derived for a specific structural component are presented employing a digital com— puter for certain problems and using conventional direct methods for other problems when feasible. Convergence of Philip C. Rymers solutions obtained toward known solutions is demonstrated where known solutions exist. Unsymmetric capabilities of the model are demonstrated. The problems investigated include; the buckling of a straight, pin-ended bar, including the effect of axial strain; the static buckling pressure on a circular ring, considering both normal and central loading; the buckling of a column, including the follower force type of loading; the determination of the period of oscillation of a beam; the determination of the period of oscillation of an inex— tensible circular arch; and the determination of the pure impulse required to cause snap-buckling of an arch. The investigation' of static problems is sufficiently comprehensive, the method so straightforward, and the results obtained compare so favorably with known solutions presented in the literature, that the model used is established as a reliable practical model. When the model is used in the solution of dynamic problems, particularly the determination of oscillation periods of beams, predictable accuracy is demonstrated. The physical model, when applied to arches, is shown to respond to unsymmetric loading or unsymmetric physical properties and geometry. However, no solutions of such problems have been found in the literature for comparison purposes. The computer program and flow chart used in the arch analysis are included as the Appendix of this thesis. APPLICATION OF A DISCRETE ELEMENT MODEL TO THE STUDY OF THE STATIC AND DYNAMIC STABILITY OF BARS, ARCHES, AND RINGS BY .‘ A. “ Philip C. Rymers A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Metallurgy, Mechanics and Material Science 1968 ACKNOWLEDGMENTS I offer my most grateful thanks to the many people who have helped me in this project. Especial thanks are extended to Dr. W. A. Bradley for his helpful guidance and counsel throughout the study which made this thesis possible, and to the members of my guidance committee, Dr. G. E. Mase, and Dr. L. E. Malvern, and Dr. P. Wong of the mathematics department. Additional mention is also due those at the Univer- sity of Nevada, eSpecially J. Bird and Dr. G. E. Sutton, whose comments and constant encouragement have been a great help to me. And finally, the patience and understanding of my wife and family must be included in any acknowledgment of assistance. Their steady confidence has certainly contri- buted to the successful completion of this project. ii TABLE OF CONTENTS ACKNOWLEDGMENTS O O O O O O O O O O O O O O O O 0 LI ST OF TABLES O O O 0 O O O O O O O O O O O O 0 LIST OF ILLUSTRATIONS O I O O O O O O O O O C O 0 CHAPTER I. II. III. INTRODUCTION 0 O O O O O O O O O O O O I 1.1 1.2 1.3 Motivation 0 O O O O O I O O O O O 0 History of the Problem . . . . . . . Objective . . . . . . . . . . . . . DEVELOPMENT OF THE MODEL . . . . . . . . N NNNN O 0.. 0 U1 DWNH NNN mum Method of Approach . . . . . . . Formation of the Physical Model . . Equations of Motion . . . . . . . Shear Force, Bending Moment, and Normal Force . . . . . . . . . Axial Deformation and Rotation of a Bar . . . . . . . . . . . . . . Reduced Form of the Equations . . . The Column Problem . . . . . . . . . Linearized Form of the Equations . . APPLICATIONS AND RESULTS . . . . . . . . 3.1 Axial Strain Effect on the Buckling of a Straight, Pin-Ended Bar . . 3.1.1 Bifurcation Values . . . . 3 1.2 Deformed Shape of a Pin-Ended Bar 0 O O O O O O O O O 0 Static Analysis of a Circular Ring . 3.2.1 Central Loading of a Circular Ring 0 O O O O O O O O O 0 iii Page ii vii H m whopa 00K) 11 14 18 25 28 29 29 32 38 38 CHAPTER Page 3.2.2 Methods of Pressure Distribution . . . . . . . . 40 3.2.3 Normal Loading of a Circular Ring . . . . . . . . . . . . 41 3.3 Follower Force Problem . . . . . . . 43 3.4 Computer Solutions of Nonlinear PrOblemS O O O O O O O O O O O O O 50 3.4.1 Application to the Period of Vibration of Beams . . . . . 54 3.4.2 Application to a Circular Arch . . . . . . . . . . . . 57 3.4.3 Application to a Parabolic Arch . . . . . . . . . . . . 59 3.4.4 Application to a Column . . . . 81 IV. DISCUSSION, CONCLUSIONS,AND FUTURE STUDIES 0 O O O O O O O O O O O O O O 84 4.1 The Problem Summary . . . . . . . . . 84 4.2 Results of the Analysis . . . . . . . 85 4.3 Conclusions . . . . . . . . . . . . . 86 4.4 Future Studies . . . . . . . . . . . . 88 APPENDICES O O O O O O O O O O O O O O I O O O O O 9 0 A. Flow Chart 0 O O O O O O O O O O O O O 91 B. Computer Program Listing . . . . . . . 94 BIBLIOGRAPHY O O O O O O O O O O O O O O O O O O O 103 iv Table 3.1. LIST OF TABLES Numerical values of the ratio of applied load to critical buckling load for a pin-ended column considering the in- fluence of axial strain . . . . . . . . Numerical values of deformed shape of the two element pin-ended column consider- ing axial strain . . . . . . . . . . . Numerical values of deformed shape of the three element pin-ended column consider- ing axial strain . . . . . . . . . . . Critical pressure coefficients for a stat- ically loaded uniform circular ring under uniform external pressure . . . . Buckling and natural frequency data for follower force problem . . . . . . . . Typical computer output. Five mass point beam oscillating in its fundamental mode in the region of the one quarter period of oscillation time . . . . . . Results of computations of periods of oscillations for transversely loaded beams O O O I O O O O O O O O O O O O I Results of extrapolation of data from Table 3.7. Transversely loaded beam excited in the fundamental mode . . . . Inextensible circular arch of seven mass points. Fundamental period of oscil- lation for each mass point that moves and deviation from linear continuum results . . . . . . . . . . . . . . . . Page 31 34 37 42 48 54 54 57 58 Table 3.10. 3.11. 3.13. 3.14. 3.15. 3.16. 3.17. One degree of freedom arch--pure impulse One degree of freedom arch--sudden1y applied force 0 O O O O I O O O O O 0 One degree of freedom arch--impulsive force 0 O O I O I O O O O O O O O I 0 One degree of freedom arch--load applied over 40 time cycles . . . . . . . . . One degree of freedom arch--load applied over 80 time cycles . . . . . . . . . Two degrees of freedom arch--unsymmetri- cal mass--initia1 velocity of 1000 in/Sec O I O O O O I O O O O O O O 0 Control problem--one degree of freedom-- initial velocity of 1000 in/sec . . . Two degrees of freedom arch--bars of unequal stiffness--initial velocity Of 1000 in/Sec O O O O O O O O O O O 0 Two degrees of freedom arch--horizonta1 force of 200 pounds applied to the mass point--initial velocity of 1000 in/sec . . . . . . . . . . . . . . . Three mass point, six degree of freedom arch-—vertica1 impulse at tOp mass point of 19000 in/sec . . . . . . . . vi Page 63 63 64 64 65 70 70 71 71 80 3.8 3.9 3.10 3.11 LIST OF ILLUSTRATIONS Typical mass point arrangement . . . . . . Force and moment distribution on a typical bar 0 O O O I O O O O O C O O O O I O 0 Displacement of a typical bar . . . . . . The column model . . . . . . . . . . . . . The two-bar model of a deformed, pin-ended bar 0 O O O O O O O O O O O O O O O O 0 Bar with axial deformations considered, post buckling . . . . . . . . . . . . . Bar with axial deformations considered, bifurcation point and post buckling . . Four and six bar approximations to a circular ring . . . . . . . . . . . . . Method of force determination at a joint for a uniformly loaded ring . . . . . . Two mass point approximation of the fol- lower load problem . . . . . . . . . . Influence of mass points added for fol- lower force column problem . . . . . . Acceleration, velocity and displacement vs. time curves for a representative mass point . . . . . . . . . . . . . . Inextensible arch of seven mass points . . Single degree of freedom arch, pure impulse loading . . . . . . . . . . . . Single degree of freedom arch, impulsive force 0 O O O O O O O I O O O O O O O 0 vii Page 10 12 15 22 32 35 36 39 41 44 47 52 59 65 66 Figure Page 3.12 Single degree of freedom arch'suddenly applied force . . . . . . . . . . . . . . 66 3.13 Single degree of freedom arch, force increasing to maximum over 40 time increments . . . . . . . . . . . . . . . 67 3.14 Single degree of freedom arch, force increasing to maximum over 80 time increments . . . . . . . . . . . . . . . 67 3.15 Unsymmetrical behavior of single degree of freedom arch approximation . . . . . . 69 3.16 Three mass point arch, restricted mOtion O O O O O O O O O O O O O O O O O 75 3.17 Five mass point arch, restricted motion . . 76 3.18 Three mass point arch, mass points free . . 77 3.19 Displacement response of three and five mass point arch approximations . . . . . 78 viii CHAPTER I INTRODUCTION 1.1 Motivation The tOpic for investigation presented here origi- nated from the First International Conference on Dynamic Stability of Structures held at Northwestern University in 1965. Of particular interest at that time was a pres- entation by Budiansky [17] in which a dynamic buckling analysis of a single mass point, two bar approximation to an arch was considered. Efforts to extend this analysis to include more mass points and bars as well as to beams and columns led to the use of a digital computer as a tool to solve the equations which arise. This approach was partially motivated by the work of Eppink and Veletsos [24] in which the dynamic analysis of a symmetric, shallow circular arch is presented. Prior to these events, attention was drawn to the general tOpic of elastic stability in a formal course of instruction presented by the author's graduate advisor, Dr. W. A. Bradley. Subsequent study and investigation have served to greatly increase interest in the subject of static and dynamic elastic stability, resulting in the thesis that is presented here. 1.2 History of the Problem, The investigation of the stability of structural components such as columns, arches, rings, and plates has occupied the attention of mathematicians and engineers for many years, the column investigations of Euler in 1759 being among the earliest recorded studies. Until rela- tively recent times, however, the investigations have centered on a static analysis of the problem. In 1952, Ziegler [41] formulated a set of criteria for the solution of stability problems which expresses the limit of appli- cability of static analysis. This work was further de- velOped by Ziegler [42] when he reviewed the classification of force systems and proposed in particular that a non- conservative force system requires special treatment, that of dynamic analysis. Since this work appeared, inter- est in dynamic stability has steadily increased. Many Russian scientists have made contributions to the dynamic analysis of structures, perhaps the most nota- ble among them being Bolotin [l], [2], who has collected an extensive bibliography of the efforts of the Russian scientists in [l]. The work of Koiter [30], which is reviewed in English in a later article by Koiter [29], firmly estab— lished the need for consideration of imperfections in structural analysis, and provides a method for determining the influence of these effects. Among other recent contributions are the works of Hermann and Bungay [25] and Herrmann and Jong [26] in investigating the class of problems in which structures are subjected to nonconservative forces, and the work of Budiansky and Hutchinson [18] in investigating impulse loading of arches. In addition, the work of Hsu [27] in defining in a mathematical sense the concept of dynamic stability, and Humphreys [28] should be noted. Particularly pertinent to this present study are the papers by Newmark [35] and Eppink and Veletsos [24]. Newmark has established a method of analyzing dynamic equations, of the type to be considered here, for linear systems of equations, which is particularly well suited to the modern, high speed digital computer, and Eppink and Veletsos have applied this method to the solution of a circular arch problem, which is of a considerably more restrictive nature than some of the problems to be con- sidered here. The dynamic analysis of structures and structural components, in particular those studies relating to buckling, has progressed to such an extent that a defini- tion of the stability of a structure has become difficult. Stability in a mathematical sense, and stability in a physical sense for a particular structural component may not be identical concepts. For example, as pointed out by Koiter [30] in his paper dealing with imperfections, a structure may be stable in a given equilibrium position, yet neighboring configurations achieved by imperfections in the structure may not be stable. In particular, there may be some lower loading level for which the structure tends to buckle. Thus, even though a particular solution to a problem is stable in a mathematical sense, the physi- cal structure may exhibit tendencies to behave in an unsta- ble manner due to imperfections. Other papers of interest on the subject of elastic stability include, for rings, the work of Ball [13], Boresi [l6], DenHartog [21], Lind [32], Schreyer and Masur [36], and Stuiver [40]. Those additional papers consider— ing principally beams and columns include works by Bailey [12], Beatty [l4], Beatty and Hook [15], Burgreen [19], Dengler and Goland [22], Lindberg [33], Morris [34], Seames and Conway [37], and Sevin [38]. The majority of the works referenced here and listed in the bibliography approach the problem as a continuum rather than by the method of segmenting into systems of bars and masses as is done here. The usual method of solution in these works is other than numerical. The survey article by Eisley [23] provides addi- tional historical comment and affords a more complete bib- , liography than that undertaken here. 1.3 Objective For the purpose of this investigation, the stability of a structural component will be considered assured if the structure, when displaced to any adjacent configuration. by any means whatever, is not subject to a catastrophic failure. Hence, an arch does not "snap through" if it is displaced to an adjacent configuration with no change in the loading of the arch. This investigation is then con- cerned with determining limiting conditions of stable be- havior in-a physical sense, for a particular structural configuration, in addition to determining the stable re- sponse of the structural component to dynamic loading con- ditions. Additionally, in some problems, static behavior will be investigated. No attempt is made to define sta- bility, and no claim to a definition of stability is to be implied. The purpose of the study, the results of which are given in this thesis, was to develop and evaluate a simple model for study of planar structural components such as beams, rings, columns,and arches subjected to conservative and nonconservative static and dynamic loadings. In study- ing the static and dynamic response of such structures and structural components, two alternate approaches are: a) To idealize the component, write the differential equations of motion or equilibrium,and solve these- equations. In many cases, closed form solutions are not easy, or in fact possible, to find with the mathematical skills presently available. Nu- merical methods, such as finite differences, are often used, satisfying an approximation to the governing differential equations at a finite number of points. b) To model the physical problem by a finite number of discrete elements, and obtain exact solutions for this model. Which of these two methods is preferable is at least partially dictated by the preference of the problem solver, and in some cases the two methods are entirely equivalent. It is possible, however, that it is easier to intuitively estimate the degree to which the discrete model represents the physical problem than it is to estimate the nearness of an approximate solution to the exact solution of the differential equations. The method of approach used here is the second of those listed, to represent the structural component as a system of bars and masses connected by torsion springs, to formulate the mathematical conditions describing the subsequent response of this system of bars and masses to an external excitation, and to determine a solution to this system of equations which is consistent with the constraints imposed by the geometry of the struc— tural component and its method of support. The interpre- tation of the solutions so arrived at will provide insight into the elastic stability of the structural component and its fundamental natural frequency when dynamic analysis is indicated. Specific structural components studied and their method of loading are: 1. Straight columns subjected to conservative loads and the nonconservative follower type of load. 2. Circular rings subjected to normal pressure and central pressure. 3. Beams subjected to initial displacements but with- out external loading. 4. Arches subjected to symmetric and unsymmetric geometry and load conditions and including pure impulse loadings. The method of solution used is evaluated by com— parison of the solutions obtained to known results, where these results are available. It was initially intended to study the effect of transverse shear deformation, axial deformations, bending deformations, and rotary inertia. However, it was not found possible to include transverse shear and rotary inertia effects with the specific solution method employed. CHAPTER II DEVELOPMENT OF THE MODEL 2.1 Method of Approach The structural member to be investigated is first replaced by a model which consists of a number of weight- less bars connected together by frictionless hinges at the nodal points where the mass of the structure is considered to be concentrated. Then the motion of each mass is con- sidered, employing the necessary conditions of compatibility to insure a continuous structure. As the number of bars be- comes large, the condition of the true mass distribution is approached. 2.2 Formation of the Physical Model The bars are considered to have axial elastic properties, but are incapable of bending. In allowing for the influence of bending of the structural component, it is assumed that torsion springs of constant ki are at- tached across each nodal point of the model so that as adjacent bars rotate, these springs are deformed. The model is illustrated in Figure 2.1. 10 Figure 2.l--Typical mass point arrangement. Each mass point is considered to be acted upon by a system of externally applied forces, in addition to the forces and moments in the bars attached to it. The result- ant of these external forces is denoted by the force vector Pi in Figure 2.1. Also in Figure 2.1, the angle ¢i denotes the angle the ith bar forms with the horizontal when the structure is in the original undeformed position, the angle th wi denotes the total rotation of the i bar into its cur- rent position, Ni denotes the axial force, Si the shear th force, and Mi the bending moment attributed to the i bar. 2.3 Equations of Motion Using the model shown in Figure 2.1, the dynamic equations of motion are written for the coordinate direc- tions shown. m.u. = Xi and mivi = Y. (2.1) 11 In equations (2.1) the dot is used to indicate differentia- tion with respect to time. The forces Xi and Yi are ex- pressed by the equations (2.2). Xi = ‘N1 °°S(¢i + u’1) + N1+1 °°S‘¢i+1 + 1"14.1) Si Sin(¢i + wi) + Si+l Sin(4’1“ + 1"1+1) + Pi sinei (2.2) Yi = -N1 sin(¢i + wi) + Ni+l sin(¢i+l + wi+l) + Si Coswi + 1pi) ' S1+1 C°s(¢i+1 + wi+1) P. cose. 1 1 2.4 Shear Force, Begging Moment, and Normal Force For a typical bar loaded as illustrated in Figure 2.2, the shear force-bending moment relation is derived, using a statical moment equation. The use of a statical moment equation is justified by the assumption that the bar is weightless, and thus its dynamic resistance to rota— tion can be ignored. Summing moments about the end A of the bar in Figure 2.2 gives (2.3) In equation (2.3), £1 is the original length of the bar, and 6i is the current axial elongation or contraction of 12 Figure 2.2--Force and moment distribution on a typical bar. 13 the bar due to a change in position of the ends of the bar. The bending moment is related to the rotation of the bar by the relation Mi =-ki(wi - win) , (2.4) where the spring constant ki is given by k. = 1 i . (2.5) In equation (2.5), E1 is Young's modulus of elasticity of the bar, and I1 is the area moment of inertia of the bar cross-section measured about an axis perpendicular to the x-y plane and passing through the centroidal bending axis of the bar. The axial force Ni is formed using the relation N. = a. 6. , (2.6) where ai is the axial stiffness of the bar, expressed as a = 1,1 . (2.7) 14 2.5 Axial Deformation and Rotation of a Bar There remains the determination of the extension 6i and the rotation $1 of a typical bar which results from a change in position of the end points of the bar. Figure 2.3 shows a typical bar in its original and displaced configurations, where for ease of illustration,‘ the displacements are greatly exaggerated. From Figure 2.3 we see that ‘1 cos¢i u. + (ILi + Si) cos(¢i + 01) - u. 1-1 1 (ui-l - ui) + (ii + 5i) cos(¢i + $1) , from which 6 = £1 cos¢i - (ui-l - ui) _ 2 i cos(¢i + 01) i But we can also write 2. cos¢. + (u.— u._ ) C°S(¢i+'wi) = 1 1 1 1 1 . “\/[5Licos¢i+(ui-ui_l)12+[£isin¢i+(vi-vi_l)]2 Thus we have for the change in length of the bar (2.8) __ 2 . 2‘. 51: -11 + Wszicos¢i+(ui-ui_l)] +[zis1n¢i+(vi-vi_ln Also from Figure 2.3 we have 15 tan¢. + tanw. 2.sin¢. + (v. 1 1 = 1 1 1 l-tan¢i tantpi llcos¢i + (ui - ui_i) Figure 2.3-—Displacement of a typical bar. Thus, after some computation, we can express the bar rota- tion wi in terms of the initial geometry of the bar as, w = Arctan (Vi - vi_l)cos¢i - (ui - ui_l)s1n¢i _ _ . (2.9) 1 2i + (ui ui_1)cos¢i + (vi vi_1)81n¢i Expanding equations (2.8) and (2.9) in series, and retain- ing only up to cubic terms in the relative displacement of the ends of the bar there results for the bar elonga- tion, 16 . 2 . 2 Sin ¢i 3 cos¢is1n ¢i 51:‘ui'ui-1’°°S¢i+(ui’ui-1) ‘II;“ '(“i'ui-l) 2, 2 i . 2 coszcbi 3 cosz¢isin¢i +(Vi'Vi-1’Sln¢i+‘vi‘Vi-1’ ‘EII" -(Vi-Vi-l) 2i: i cos¢i sin¢i ‘(ui‘ui-1)(Vi'vi-1) mi . 3 2 . 2 Sln oi- 3 cos ¢i s1n¢i -(u.-u. ) (v.-v. ) 1 1-1 1 1-1 2* 22. 1 2 cos3¢i- 3 cos¢i sin2¢i -(ui-ui_l)(vi-vi_l) 22 2 + ... (2.10) i Equation (2.10) has been written so that a shortening of the bar will result in a negative value for Si. The equations (2.8) and (2.9) are in a form well suited for use with a high speed digital computer, whereas equations (2.10) and (2.11) are somewhat more convenient for use with a desk calculating machine, or for hand cal- culations, in completing the problem solution by some numerical scheme. 17 For the rotation of a bar there results, sin . cos . 'n . $1 2 ¢1 Sl ¢1 wi = ’ (“i ‘ ui-l) “I;’ ' (“i ‘ ”1-1) £ 2 i 3 sin3¢i - 3 coszcbi sin¢i + (“i ' ui-l) 3 32. 1 cos¢. cos¢. sin¢. 1 2 1 1 + (Vi ' Vi-l) “I;’ + (Vi ‘ Vi-l) 2 2 i 3 cos3¢i - 3 cos¢i sin2¢i ‘ (Vi ’ Vi-l) 3 . 32. 1 cos2¢i - sin2¢i + (“i ‘ “1-1)(V1 Vi 1) 2 2 i 2 cos3¢i - 3 cosoi sin2¢i + (ui - u1 1) (v1 - Vl-l) 2 3 i 2 sin3¢i-3cosz¢isin¢i + (ui-ui_l)(vi-vi_l) 2 3 + .. (2.11) i The formulation of the problem is now complete. Using the above expressions for the axial deformation and rotation of each bar of the structural components being analyzed, the equations (2.1), with equations (2.2) through (2.7) now become a system of simultaneous differential equa- tions in the displacements of the mass points of the seg- mented model. For each mass point two equations will exist, thus, for a model containing n members and n - 1 mass points 18 there will be 2(n - 1) equations to be solved in 2(n - 1) unknown displacements. Of course, using the model de- scribed requires that n > 1. Such a system of equations is very difficult to solve formally in such a way that the resulting solution would be applicable to several classes of problems such as rings, columns, arches, and beams. Therefore a numeri- cal solution of the system is indicated. 2.6 Reduced Form of the Equations In order that the numerical solutions obtained have greater general applicability, the equations of the problem are next written in terms of the elastic constants and cer- tain geometric properties of the segmented model. Introducing the quantities, ' a' I k. ' mo _ 1 _ 1-1 _ 1 “i ‘ ET ' ki ’ k. ' mi ‘ ET ' (2'12) 1 l 1 we have for the new coefficients in the equations of the problem, I l -(‘Pi-¢i_l)+ki (wi_l-wi) N. - l:- = Ni ‘ ET “1 51 ' Si L.+67’ ' 1 l l (2.13) 19 where “i is the mass per unit volume of the structural com- ponent. In equations (2.12) and (2.13), primes denote the new coefficients used in formulating the problem. Finally then, combining equations (2.1), (2.2), (2.12), and (2.13) there results, I 8 “i = [' Ni °°S‘¢i+wi) + N1+1 C°S(¢i+1+wi+1’ v ' Ei ‘ Si Sln‘¢i+¢i’ + Si+1 Sln‘¢i+1+¢i+1’1 u—a'Tz—z i i 1 Pi + $7 Slnei (2.14) J. and I I Vi = [’ Ni Sln(¢i+wi) + N1+1 Sln(‘bi+1““’i+1) I 1 Bi + Si Sln<¢i+wi) ‘ Si+l °°S(¢i+1+¢i+1’] "“T“§ u.a. l. l l 1 Pi ‘ - a: COSEi . (2.15) Equations (2.14) and (2.15) are the equations of the prob- lem in a "dimensionally reduced" form. Note that in the event no forces Pi are applied directly to the structural component being analyzed, the mass of each nodal point is not required since the ratio Pi/mi becomes the local gravi- tational constant. Thus for problems where the motion is 20 initiated by an initial displacement or velocity, but where no forces are applied, except possibly forces having a known magnitude relative to the total weight of the structural component, the input constants consist of the elastic con- stants in axial stiffness and bending for each segment of the model, the ratio of the bending stiffness for adjacent bars of the model, the mass per unit volume of the material of the structural component, and Young's modulus of elas- ticity. Also needed are the conditions of loading, and the conditions for initiating motion for the problem being considered. It is also possible, by the introduction of the dimensionless ratios, ui _ vi u. = —— t v. = — 1 2. ' ' 1 1. ' l 1 to write equations (2.1) and (2.2) in the non—dimensional form, 8 - - fi' c (¢ +1 ) + E" ki+1 1 0 (¢ +w ) 1 ‘ 1 03 1 1 1+1 k. 2. °5 1+1 1+1 1 1+1 - S sin(¢ +0 ) + S ki+l £1 in(¢ +0 ) 1 1 1 1+1 ki £1+1 5 1+1 1+1 + PK 'n . 1 $1 81 21 and, u- - k1+1 11 V1 ="Ni Slnwiwi) + N1+1 k. 9.. 51“ (¢1+1+"’1+1) 1 1+1 k. 2. — — +1 1 +5. cos(¢.+1p.) - s. 1 cos(d>. +41. ) 1 1 1 1+1 ki 2i+l 1+1 1+1 -»Pi cosei . In these non-dimensional equations, use has been made of the notations, N 2 8.2. P.2. fi' = 1 1 § = 1 1 P = 1 1 1 ki ' 1 E1 ' 1 Ei ' 2.7 The Column Problem To solve a column problem it is first necessary to consider in detail the constraints that are implicit in equations (2.8) and (2.9). To apply these equations re- quires that the axial deformation of a bar be dependent on the position of the ends of the bar. If the structural component being considered does not possess geometric con- straints necessary to insure that the displacements can be calculated, equations (2.8) and (2.9) do not apply. Since a column lacks these constraints, some additional constraint must be imposed on the column problem formulation. One method of introducing the necessary condition is the construction of an equation of constraint which 22 would relate the column geometry and the applied external forces in a manner independent of the relations expressed previously. A second method of approach is to introduce assumptions into the problem formulation without simulta- neously introducing errors of unacceptable magnitude. The latter of these methods is employed here. Starting with equations (2.2) where £1: 0--the structural component is initially horizontal--we are X II I Nicoswi + Ni+lcoswi+l - Sis1nwi + . ‘ . + . ’ . Sl+151nw1+1 P151ne1 (2.16) Y. = - Nis1nwi + Ni+ls1nxpi+l + Sicoswi — Si+lcoswi+l - Picosei , where Ni and Si are as defined before. Now if the constraint introduced is that the bars are :inextensible,, and all of the same length, the distance between mass points does not change and we have for the bar next to the support in 9. a - £7 fi‘JLul uzquvl Figure 2.4--The column model. 23 Figure 2.4, 2 - ul= 2 - v1 or u1 = 2 - 2 - vl , also 22 - u2 - (2 — ul) =1/22 - (v2 - v1) , or and in general _ _ 2 _ y 2‘ u. — 1 + ui_l -\/2 (vi vi+1) . (2.17) Considering the rotation of the bars, and using equation (2.17) there results, -1 v. — v._1 I('i"'t""‘1"‘ z 1 1 2] ' Va ‘ (V1 ‘ V1-1) J which replaces equation (2.9). (2.18) To determine the normal forces Ni in the bars now that no axial deformation is allowed, an additional assump- tion is useful. Since the masses are numerically quite small for any real problem, and since by the very nature of a buckling problem, the accelerations of the mass points will be small in the vicinity of the actual buckling con- figuration, since then the period of oscillation tends to 24 infinity, the inertia forces in the horizontal direction will be ignored. The validity of this assumption can only be borneout by the results of the computations which appear later in Section 3.4.4 Application to a Column. Using this assumption and the first of equations (2.16) there results _ 1 Ni _ cosxpi [Picos(Bwi) + Ni+lcoswi+l - 8131nwi + Si+ls1nwi+ll . (2.19) In equation (2.19) the angle Bwi replaces the angle ei in equations (2.16) since it is anticipated that the column load can be represented as a function of the angle of inclination of the adjacent bar. Note that if B = 1, the problem is that of the nonconservative follower load problem, while if B = 0, we have the conservative load problem, or the Euler column problem. To summarize then, the equations governing the problem represented in this section are, 25 miv. = - Nis1n1pi + Ni+ls1n1pi+l + Sicoswi - Si+1coswi+l + Pis1ani , -1 v1 ' v1-1 w. = tan ' (2.20) 1 V22 - (v. - v. )2 1 1-1 N = ——£—— [P cost + N cosw i cosxpi i i i+l i+l - Si31nwi + S1+131n¢1+11 and s=-’i(w -2w+w) i 2 1+1 i i-l ' 2.8 Linearized Form of the Equations As a first attempt to establish the validity‘of the use of the model which results in equations (2.2), consider the following reduction of these equations to linear form. For a straight bar, ¢i = 0. Further let 8i = 90° and hence only axial forces are applied. This gives the equations (2.21). -Nicoswi+Ni+1coswi+l-Sis1nwi+si+ls1nwi+1+Pi = Xi (2.21) -Nis1nwi+Ni+131nwi+1+Sicoswi-Si+lcoswi+l = Yi 26 Now from equations (2.3) through (2.7) where the bar is of constant cross-section and bending stiffness, we have EA N1 - 7r 61 ' and ‘ (2.22) _ EI _ _ _ S1 ' 2(2 + 51) [ (21 w1+1) + (wi_1 wi)] . Placing equations (2.22) in (2.21) there results, with2.= h. (w. 'W- )-(W-'W- ) EA EI 1+1 1+2 1 1+1 . TF‘51+1°°3¢1+1 ' 51°°Sw1’ ' I?“ [ h + 51+1 )Slnw1+1 (w- ‘ W- ) - (¢-_ ‘ W-) _ [ 1 1+1h + 8.1 1 1 ] Sinwi + Pi = xi I 1 . and (2.23) - (¢-'¢- )-(¢- '¢-) EA . . EI 1 1+1 1-1 1 TF‘51+181n¢1+1‘5151nW1)‘TT'[i h + 51.“ J °°Sw1 [(w1+1 ' w1+2) ' (W1 ' w1+1)] h + 51+1 c°S“’1+1 = Y' ° . Next, assume small mot1on so that 61 = 61+1 = 6 and introduce only the linear parts of equations (2.10) and (2.11) to form, 27 EA E1 2 _ IT’QO) "ir [0‘w1 ’1 + P1 ‘ x1 ' or Pi = Xi , (2.24) and Nh(h + 6)(vi+1-2vi+vi_l) - EI(vi+2-4vi+l+6vi-4vi_l + Vi-Z) = Yi . (2.25) Then if the horizontal motions are neglected we have, with N =._ P, (V1-1'2V1+V1+1) _ BI (V1-2'4V1-1+6V1’4V1+1+V1+2’ 21 - P h2 h (2.26) The left-hand side of equation (2.26) is the finite difference approximation-to the first two terms of P 32” + EI 34V + 91 32V = 0 - (2 27) 3x2, 3x 9 at2 Equation (2.27) represents, in the sign convention used here, the equation of motion for the vibration of a' beam acted upon by axial forces as derived in Timoshenko [10], pages 374-375. CHAPTER III APPLICATIONS AND RESULTS In this chapter, the theory developed in the pre- ceding chapter is to be applied to a variety of problems, both static and dynamic. The primary objectives of this chapter are to illustrate the application of the theory, to compare results to known solutions so as to verify the ac- curacy of the equations developed from the theory and to extend the applications to problem areas where no known results are available. The problem solutions presented are of two types. First, certain problems having known so- lutions will be solved using sufficiently simplified systems of bars and masses so that hand calculations can be presented. These solutions will provide some insight into the conver- gence of the system of equations to known solutions for the various problems. Finally, problems of greater complexity will be solved using high speed digital computers as tools. These problems will be predominantly problems which require dynamic analysis, and will also include nonlinear effects which will be ignored in the first category of problems. In this way will be indicated some idea of the range of applic~ ability of the model as well as the accuracy that can be generally expected in solving problems of the type considered here. 28 29 3.1 Axial Strain Effect on the Buckling of a Straight, Pin-Ended Bar This problem will be approached in two ways so that bifurcation and post buckling behavior can both be observed. Results are compared to the work of Huddleston [43] wherein an exact solution to this problem is presented. 3.1.1 Bifurcation Values Starting with equations (2.24) and (2.25) where we set nv. _ _ _ 1 _ _ P2 _ 2 x1 ‘ Y1 ‘ 0 ' V1 ‘ ‘2’ ' 5 ‘ nAE ' h ‘ n ' there results Pi = 0 , and P A22 P ' g2; (T (1 " A‘E) (V1+1 ‘ 2"1 + V1-1) + V1+2 ' 4"1+1 + 6"1 Next introduce the notation I l P c=——2. , 1 =-§3—(-5)(1-X§) . A2 n AE and we have Mvi+1 - 2vi + Vi-l) + vi+2 - 4V1+1 + 6vi - 4Vi_l + V1-2 = 0 (3.1) 30 Using equation (3.1), A can be determined for any selected displaced configuration. We note that in general we have EI a1 8 P = a and = B , or a = -' , 27 A22 C we also have, P _ -A_'E- _ B 0 Using this notation, the equation relating B and.) becomes 2 2 B - B + Cn A = 0 . (3.2) which has the roots 2 B _ 1 l _ N l E'IEfl/E'z T:- - ‘3'” Then with 1 determined from equation (3.1), the buckled con- figuration equation, g'= a, can be determined from equation (3.3). So finally we have _P_ Por 8— d (3.4) a cue where d is a factor defining buckling of the segmented model whose solution is sought. For example, with n = 2, d = 8. Table 3.1 displays the results of applying equations (3.1), (3.3), and (3.4) for several values of n, where the value of C is taken to be C = 0.02. 31 Table 3.l.--Numerical values of the ratio of applied load to critical buckling load for a pin-ended column considering the influence of axial strain. C = 0.02. n a P/Pcr 2 10.0000 40.0000 1.2500 5.0000 3 11.7712 38.2288 1.3079 4.2476 4 12.4951 37.5049 1.3331 4.0015 5 12.8533 37.1467 1.3464 3.8911 6 13.0555 36.9446 1.3534 3.8298 7 13.1783 36.8217 1.3579 3.7951 Comparison to Huddleston [43], where the exact solution is given to be 1.371 and 3.695 for the two values of the ratio at the bifurcation point, the values in Table 3.1 show that the results for n = 7 are in error by -0.948% and 2.679% respectively. Using the values of Table 3.1 and em- ploying three point extrapolation for n = 5, 6, and 7, there results 1.3715 and 3.7002 respectively for the bifurcation values. The extrapolation scheme used is that referred to as Richardson's extrapolation. Two schemes will be used in this thesis, the two point extrapolation, and the three point extrapolation, 4 4 “f _ A n1 f (n? - nian? - no?) 1 (an - nianiz - no?) n 4 C + A ' I c (nfz _ n02)‘n12 _ nCZ) A=A 32 where the subscripts f, i, and c refer to the fine, inter- mediate and coarse grid sizes, respectively. 3.1.2 Deformed Shape of a Pin- Ended’Bar The shape of the deformed bar will be characterized by two values, A/2 and 6/2, which are illustrated in Figure 3.1, the two bar model. fly P 1P1 (P1 P __|g.._ Figure 3.1.--The two-bar model of a deformed, pinrended bar. For the bar of Figure 3.1, using equations (2.23) in which we let Xi = Yi = 0 at point 1, there results, -H (61 coswl - 6l coswl) - ‘H h + 51 sin(-wl) " h+ 5 ] sinlpl + 0 = 0 l and 33 (1D “(w-(w -w) (61 sin(-wl) — 61 sinwl) - §% 1 i + &l l l cosw1 11 H ("Wl + I("1) " (191 + 1111) - h + 31 cost!)l = 0 . In these equations, the conditions of symmetry dictate that 61 = 62,01 = -w2 and thus the first equation is an identity. The second equation reduces, with the help of the relations Nh 5 = AE and P cosw1 + N = 0 , to P _ zwl f EI - sinw P cosw ' 1 h2 1 _ 1 AE Further introducing 2 h=£ ’ P =TTEI ’ i=2... I C=_:£2. ' 2 cr 22 1T2 Pcr A2 there results a2 _ a + lswl = 0 Ccos$i Csin2¢l ’ which has the roots a = ————l—— 1 1 /1 - 32c¢ cotw . (3.5) 2Ccosw1 1 1 Additionally, from the geometry of Figure 3.1 we can write sinwl ——§—- (1 - Ca coswl) and (3.6) A 2 nu» I - 1 - coswl (l - Ca coswl) . 34 Proceeding in a similar way for the three element approximation to the bar, there results in place of equations (3.5) and (3.6) the following set of equations: _ 1 _ , a - 2 C COSWI [l 1 V1 36VC 01 cdtwl] A sintp1 i- = 3 (1 " CC). 005191) I (3’7) %-= % [2 - 2 coswl + Ca(1 + 2 cos2 wl)] . Some values for the ratios A/2 and 6/2 of equations (3.6) and (3.7) are given for specific values of P/Pcr' with C = 0.02, in Tables 3.2 and 3.3 and presented graphically in Figures 3.2 and 3.3. Again it is seen that convergence is rapid and results compare favorably to those of Huddleston. Additional references to this problem are the works of Beatty [l4] and Beatty and Hook [15]. In [15] experimental results of tests on rubber bars are presented. Table 3.2.--Numerical values of deformed shape of the two element pin—ended column considering axial strain. C = 0.02- 11 P/Pcr A/2 6/2 0 1.2500 5.0000 0.0000 0.0000 .2000. .8000 .2 1.2528 5.1243 .0798 .0195 .2125 .8075 .5 1.2691 5.8528 .1970 .0427 .2788 .8436 1.0 1.3447 10.2230 .3718 .0489 .5225 .9372 1.5 1.5301 86.8216 .4901 .0086 .9305 .9988 2.0 1.9470 - .5136 - 1.4701 .9460 2.5 3.0134 - .4148 — 2.1106 .6905 35 1000‘] I 0 I : 1 1L 1 1 0 0.1 0.2 0.3 0.4 0.5 A 2 2 element 4 element 0 D 3 element Huddleston 13 ————— o Figure 3.2.--Bar with axial deformations considered, post buckling. 36 7.0-» ‘f . 6.0‘)’ . I l 5.0"- p 7 ___ _. .” Pcr 4.d+ ./ 3.00— ‘1 Q I O 2.0-”- . __ A 2 element-Lo 4 elementfl l'O'L 3 1 H ddl e ementA _1_1_ _es_tgn_o 1 1 1 i f i 0 0.5 1.0 1.5 2.0 2.5 Figure 3.3.--Bar with axial deformations considered, bifurcation point and post buckling. 37 Table 3.3.--Numerica1 values of deformed shape of the three element pin-ended column considering axial strain. c = 0.02. wl P/Pcr A/2 6/2 0 1.3079 4.2476 0.0000 0.0000 .2354 .7646 .2 1.3090 4.3560 .0509 .0153 .2427 .7774 .5 1.3169 5.0137 .1266 .0332 .2823 .8458 1.0 1.3713 8.9111 .2431 .0374 .4368 1.1533 1.5 1.5336 - .3260 - .7124 - 2.0 1.9226 - .3468 - 1.0994 - 2.5 2.9350 - .2839 - 1.6029 - A careful examination of the value C = 0.02 used here may lead to some confusion. If it is assumed that the bar is of square cross-section, we are led to the result that the width of the bar approximates one-half the bar length. If, however, a wide flange section is considered, with the area equally distributed between the flanges, and with h the depth of the section, we find the ratio of h to 2 to be about 0.283, which is on the fringe of practicality. Also, the values of the ratios 6/2 and A/2 are permitted to grow to quite large values. Of course when 6/2 becomes greater than unity, the bar is subjected to tension, and hence with axial deformations considered, this ratio will increase. When these large deformation states are reached, the cross- section distortion will be appreciable and will contribute to the load values. These effects are not considered here, and so the problem takes on a flavor of the academic problem for displacements of such magnitude. 38 If stress considerations would be included in the analysis, and if plastic deformations were accounted for, these large motions would be viewed in a more realistic environment. 3.2 Static Analysis of a Circular Ring The circular ring problem will be investigated under the assumptions of small deformation and inextensibility of the ring. Two types of loading will be considered. First, a uniform external pressure will be assumed to be always di- rected toward the center of the ring, and second, this uni- form external pressure will be considered to remain normal to the ring in all deformed and undeformed configurations. Each of these problems has been given considerable attention and solutions to them are known, as for example the work of Boresi [16]. 3.2.1 Central Loading of a CIrcular Ring The problem of the circular ring under uniform cen- tral external pressure will be considered in detail for the six bar approximation as illustrated in Figure 3.4(b). The dotted lines represent the deformed configuration, assumed known beforehand. The equilibrium of each joint is then considered in order to derive the equations of equilibrium for this deformed configuration. 39 Figure 3.4.--Four and six bar approximations to a circular ring. It is equivalent to use equations (2.2) for vanishing acceleration, since the equations (2.2) are the equations of equilibrium for each node point. We have at node a F + 2Nl cos(60 + 0) - 33% (30) sin(60 + 0) = 0 , _(3.8) where N1 is the force in bar ab and the ring is assumed to be of constant bending stiffness. At the node b there results 1 and (3.9) N2 + F cosLaob - N cos(60 + (p) + E;- (3111) sin(60 + 11)) = 0 2 -N1 sin(60 + w) — F sin [_aob - Eé-(3W) cos(60 + w) = 0 , 2 where N2 is the force in the bar bc. 40 Combining these three equations so as to eliminate N and N2, and considering w to be so small that sin 0 = w l and cos 0 = 1, there results 4EI F = . (3.10) £2 3.2.2 Methods of Pressure Distribution There are two ways to determine the force F at each joint that are considered here. The force is considered to be a function of the externally applied pressure on the ring in each method. In the first of these methods, it is con- sidered that each joint carries 1/6 of the total pressure for the six element ring problem without regard to the in- clination of the bars adjacent to the joint. In the second method, the angle of the bars is considered as in Figure 3.5, where a typical joint is illustrated, and where symmetry is assumed at a joint. The results of these two approximations being used along with equation (3.10) gives two values of pressure required to maintain the assumed equilibrium position, _ EI _ EI PCI' - 3.0558 g? and PCI' - 3.6950 3 I (3011) where R is the radius of the undeformed ring. Before comments about the accuracy of this result are to be made, it will first be established whether this method of approach tends to converge. This requires the analysis of the same problem using different numbers of bars to approximate the circular ring. However, before establishing 41 P p; 91¢ I p2cos[90-(6+0H. 6+0 ‘5', 9% Figure 3.5.--Method of force determination at a joint for a uniformly loaded ring. the convergence of this model, we first consider the same six bar approximation to a circular ring under the influ- ence of normal forces. 3.2.3 Normal Loading of a Circular Ring To approximate normal uniform external pressure loading of the circular ring, it is assumed that the force at each joint always acts along the bisector of that joint. Using this approximation there results at joints a and b the following equations, derived again from equilibrium considerations. At the node or joint a F + 2Nl cos(60 + 0) - g‘E—Jiwsimeo + 0) = o , (3.12) at b 42 3EI F cos(60 - £0) - N cos(60 + 0) + N + ———U)sin(60 + 0) = 0 2 l 2 22 and (3.13) . 1 . 3EI _ -F 31n(60 - 50) - Nl 51n(60 + 0) - ——§0)cos(60 + 0) - 0 . 2 Once again, eliminating N1 and N2, there results 9 22 When the two methods of pressure distribution are used again there results _ EI _ EI Pcr — 2.5465 1:3- and Pcr — 3.0793 ? . (3.15) These typical calculations are repeated for each approximation, the foregoing being merely illustrations of the method used. Table 3.4 lists the critical pressure co- efficients that result from applying this method to 4, 6, 8, and 10 bar approximations, as well as the results of extrapolations of the data gathered. Table 3.4.--Critical pressure coefficients for a statically loaded uniform circular ring under uniform ex- ternal pressure. Number of Bars Central Pressure Normal Pressure 4 1.8006 2.8284 1.8006 2.8284 6 3.0558 3.6950 2.5465 3.0793 8 3.6398 4.0428 2.7563 3.0615 10 3.9457 4.2177 2.8526 3.0493 Extrapolation Results 4,6 4.0600 4.3883 3.1432 3.2800 6,8 4.3906 4.4900 3.0260 3.0386 4,6,8 4.5009 4.5239 2.9861 2.9582 6,8,10 4.5451 4.5504 3.0225 3.0214 43 The solutions that are accepted as correct for the two problems considered here are, as found in the paper by Boresi [16], for central pressure, 4.5, and for normal pres- sure, 3.0. The above extrapolated results differ from these by 1.002%, l.1200%, 0.75%, and 0.7133% for the 6, 8, 10 bar extrapolations respectively. Additional information perti- nent to the validity of this model will be reported in a later section in dynamic, nonlinear analysis of the circular ring problem, section 3.11 3.3 Follower Force Problem The first dynamic analysis problem to be presented is the follower force problem. Historically this problem is of considerable interest. The original prOposal of this problem led to a paradox. Static analysis yielded the re- sult that a column subjected to a force which always remained tangent to the end of the column would not collapse under any force, regardless of the magnitude of the force. Subsequently, it was shown that dynamic analysis produces a solution, as given in references [2], [41], [42], which has been further amplified in later works, [25] and [26], to show the signifi— cance of dynamic analysis. This problem is the historical beginning of the modern day dynamic buckling studies. Consider the two bar, two mass point model of a column illustrated in Figure 3.6. This is the same model, exclusive of damping, as that considered by Herrmann and 44 Jong [26]. In this problem, the mass is not considered to be constant along the length of the column. p 02 m 01 2m 3 11. 1h h} / Figure 3.6.--Two mass point approximation of the follower load problem. Applying equations (2.20), where small motion and the illustrated mass distribution is used, requires m = j-wh and N1 = N2 = -P . There results P(w - w ) + EI(20 -3w ) = 2mh0 l 2 E2 ' 2 1 l and P0 + El (0 - w > - P0 = mh (0 + 0 ) 231 2 2 1 2 ' which become, after introducing the assumption of harmonic motion in the form, and the notation 2 3 2 _ Ph 2 _ mh w a 7 21— ' B - -—§1— . (3°17) 45 the state equations for the problem, (282 + 3 - a) (-2 + 0) Al (82 - 1) (82 + 1) A2 0 Setting the determinant of the coefficient matrix to zero results in the quadratic in 82: 284 + (7 - 2a)82 + 1 = o , which is the result arrived at by Herrmann and Jong [26]. This quadratic has the solution, 82 = "7 2.20” i 1 ‘\/(7 - 20.)2 -8 . (3.18) The value of a for which 8 is a pure imaginary number is the desired solution for the critical load, since for any other value of a in equation (3.18), there results a value of 8 having a non—zero real part and the response 0i of equation (3.16) then becomes either a damped oscillation or an ever increasing oscillation, depending on the sign of the real part of 8. Hence, the ratio of the applied load to the actual critical load becomes _P_. CI = 3.3813 . Using equations (2.2), (2.3), (2.4), and (2.5) ig- noring axial deformations of the bar segments, and consider- ing as our example the three mass point approximation, there results the set of equations in the slope angles 01, 46 (4 - (1)0l + (-3 + (1)02 + 03 = -801 -01 + (3 - a)02 + (-2 + a)03 -8I I3. ‘ I U m I I 0" l> 5~~ . I 16’ I’A/I I I fall; ;‘;-'-"-""" I ' I 5%" If. I | HI I-IJ l u c | ' l xl l0 I 1 Fill U l i I 1 i 1' [tr 0 5 8.0011)12.402 15.088 P2 EI Figure 3.7.--Inf1uence of mass points added for follower force column problem. 48 is prOportional to the value of the force applied at buckling, or when the oscillation frequency is zero. Table 3.5.--Buckling and natural frequency data for follower force problem. Number EI of Bars acr w/ EAEZ for a = O 2 13:3: 3 12.40 1333; 4 15.09 13:2: Extrapolation, using Richardson's three point for- mula gives the value for the buckling parameter of a = 19.414 , cr which compares to the value “or = 20.05 given in [2]. One other extrapolated result is of interest. Using the smaller values listed in the third column of Table 3.5 for each problem, extrapolation gives the number 3.50 which compares to the value 3.52 as given in reference [10]. This problem is easily solved using three other in- dependent methods, all of which give the same results. These methods are finite differences using equation (2.26) and the appropriate diSplacement boundary conditions for a clamped- free bar, and application of the Lagrange equations of 49 dynamics to the discrete model in terms of either slopes or displacements. The finite difference method yields, for the three mass point model used before, the equations 7vl - 4v2 + v3 + 0I(v2 - 2v1) - 8 v1 2 -4vl + 5v2 - 2v3 + 0I(vl - 2v2 + v3) - 8 v2 (3.21) 2v - 4v + 2v = 82v 1 2 3 3 ’ which, when solved for the characteristic equation again produces equation (3.19). Expressing kinetic and potential energy in terms of the lepes, and using virtual work to derive expressions for the generalized forces gives for the three mass point ap- proximation to the problem, the kinetic energy, with w representing the mass per unit length. _ l - 2 - 2 . 2 T — -2--Ah[vl + v2 + v3'] or, _ l 3 - 2 - - 2 1 ° ' ° 2 the potential energy, 2 l v = EI/h[0l + 2- (02 - 02)2 + % (03 - 02)2] . and for the generalized forces, 01 = Ph(¢1 — $3) I 02 = Ph(w2 - $3) I Q3 = 0 0 Then employing the Lagrange equations and assuming harmonic motion, equation (3.19) is again produced. 50 The energy method can also be used with displacements as generalized coordinates rather than rotations, and again the characteristic equation (3.19) results. 3.4 Computer Solutions of Nonlinear Problems The equations developed in Chapter II will be solved by using a step-by-step numerical iteration form of integra- tion as described by Newmark [35]. A brief discussion of this method follows. It is initially assumed that the acceleration, velo- city, and displacement for each mass point is known at some time tk' It is desired to determine the values of the acceleration, velocity, and displacement for each mass point at some subsequent time t l' where (t - tk) W111 k+ k+l in the future be represented by the time increment, At. The solution is begun on the assumption that a linear variation in the acceleration will exist during the time interval At. This assumption provides a set of equa- tions from which, for known values of the acceleration at the times tk and tk+1’ computed. Using a dot to denote differentiation with re- the velocity and displacement may be spect to time, there results from these assumptions, within a time interval, X = ik + aT + bT2 and x = a + ZbT , 51 but we have the following conditions to satisfy: T = 0 , x = x k T = At , x = xk+l which give, x x _ k+l k a - xk and b 2At , therefore, x 2 . _ . ” k+l k 2 , or :2 =$<+.1—(x +$E)At (322) k+l k 2 k+l k ' ‘ However, we also can write At At x - x - . " k+l k 2 x] 1 xk’ -[0 xdr +‘[0 [Xk + xk + ( 2At )1 ] dT . which gives At2 At2 xk+l = xk + xkAt + xk ——3 + xk+l ‘6‘ . (3.23) Equations (3.22) and (3.23) are applicable to each mass point of the segmented structural component, and may be applied to the motion in either of the coordinate directions. Using equations (3.22) and (3.23), the trial velo- cities and displacements :flmr each mass point are computed in their coordinate directions for the end of the time in- terval. Using the trial displacements thus computed, the bar rotations and length changes are next computed. 52 q—-— ‘p———-— Figure 3.8.--Acceleration, velocity and displacement vs. time curves for a representative mass point. 53 The axial force and shear force are then determined, and from these trial values, plus the additional active forces applied to the structural component, the resultant force is determined for each mass point. These forces are used to compute the acceleration which should exist for the assumed trial displacements. The derived accelerations are compared to the assumed accelerations. If the two values are in satisfactory agree- ment, the problem is re-cycled to compute the values of acceleration, velocity, and displacement at the end of the next time interval. If, however, the two accelerations are not in satisfactory agreement, the problem is repeated in the same time interval with a new assumed value of accelera- tion. A satisfactory method of selecting an improved trial acceleration is to use the latest derived value of acceleration. The iteration technique described is repeated as often as necessary to bring the difference between the trial acceleration and the latest computed value of the accelera- tion within a predetermined tolerance. The number of itera- tions required depends largely on the accuracy desired, as well as the accuracy of the initially assumed acceleration, and varies within a problem from one time interval to another. The method converges rapidly enough to be entirely practical in use with a high speed digital computer such as the CDC 3600 machine. It is also possible to use this scheme on the IBM 1620 II computer, but there is a significant time 54 difference in the solution of a problem between the two machines. The flow chart for the computer program and the program listing for the arch problem, appears in the Appendix. If the length of the time interval is not selected sufficiently small, the method described here will not con- verge. Newmark [35] has given a method for determining a suitable time interval and has shown that for such an inter- val the method described will converge and provide stable solutions. However, his criteriatue derived for linear problems and is therefore not entirely satisfactory here. Newmark's method does, however, give an estimate beyond which the computer program itself can be made to seek a suitable time increment. Newmark's approximation is At < 0.389T , (3.24) where T is the smallest period of oscillation of the segmented structural component being used. 3.4.1 Application to the Period of Vibration of Beams The problem considered here was a bar whose modulus of elasticity of 10 x 106 psi. approximates that of aluminum. The beam selected had a length of 8.0 inches and cross sec— tional dimensions such that the ratio of the beam axial rigidity to bending rigidity is 3.0 in-z. The mass per unit volume was selected as 0.000246 lb—secz-in-4, again approx— imating that of an aluminum bar. These numerical data were 55 not selected to represent a particular cross section of a beam. The computer program was designed to seek out a suitable increment of time after an initial estimate was made, using the approximation given by equation (3.24). The time increment used for the five mass point approximation was 1.0 x 10-7 seconds. Data were then printed out for each 100 time increments. Typical data output for a five mass point approxima— tion of a simply supported beam appears in Table 3.6.~ The data shown are for the time spanning one quarter period of oscillation of the beam. Only three of the five mass points are recorded since symmetry makes the remaining data redun- dant. The motion was initiated by displacing the beam into an approximate half sine wave shape with the center displace- ment being 0.001255 inches from the original equilibrium position. No sag due to the weight of the beam was taken into account. Table 3.7 lists the periods of oscillation measured with this model for a simply supported beam for several values of n, the number of mass points, along with the percent deviation from the known linear solution to the continuum problem. The final data refinement is by linear interpolation between the listed values of the print out data. The results of performing this interpolation for the print out data as shown are not significantly different from the results obtained by using smaller print out time incre- ments. Table 3.7 also includes similar data for a seven mass point approximation of a beam clamped at both ends. Table 3.6.--Typical computer output. 56 Five mass point beam oscillating in its fundamental mode in the re- gion of the one quarter period of oscillation time. Time Displacement Velocity Acceleration V1 9.41389-05 -2.12268+01 -l.02727+06 9.000-05 V2 9.28613-05 5.4Sl6l+00 l.7l364+06 V3 1.88278—04 -4.24537+01 -2.05454+06 Vl - 7.26654-05 -5.21242+00 2.77043+06 1.000-04 V2 6.66924—05 ' -2.2539l+01 -4.78336+06 V3 — 1.45331-04 -l.04248+01 5.54087+06 V1 - 5.97098-05 -l.54737+00 —2.35551+06 1.100-04 V2 - 2.69315-04 -2.83320+01 4.17547+06 V3 - 1.19420-04 -3.09475+00 -4.71102+06 V1 - 1.98496-04 -2.15425+01 2.05975+05 1.200-04 V2 - 3.32956-04 7.64956+00 -l.83159+05 V3 - 3.96993-04 -4.30850+01 4.11949+05 Table 3.7.-—Results of computations of periods of oscilla- tion for transversely loaded beams. Number of Mass Computed Period Percent Points, N (Seconds) Deviation Ends Pinned 3 .0004280 21.90 5 .0003910 11.41 7 .0003790 7.98 9 .0003730 6.28 11 .0003680 4.84 19 .0003598 2.51 Ends Clamped 7 .000170 9.67 57 In Table 3.8 are listed the results of some extra- polations using the data of Table 3.7., Both two and three point extrapolations were made. Table 3.8.--Results of extrapolation of data from Table 3.7. Transversely loaded beam excited in the funda- mental mode. Mass Points Computed Period Percent Combined (Seconds) Deviation- 3,5 .000370 5.40 5,7 .000367 4.56 3,5,7 .000365 3.99 7,9 .000363 3.42 5,7,9 .000360 2.57 9,11 .000356 1.43 7,9;11 .000353 0.57 These results indicate that the problem as formulated and programmed is capable of producing satisfactory results for the determination of the period of oscillation of a beam. Other modes of motion can be achieved by apprOpriate selection of initial conditions and symmetry or asymmetry conditions in the program. 3.4.2 Application to a Circular ArCh In order to help establish further confidence in the physical and mathematical model used as well as the computer program, the linear response of the discrete model approximat- ing the inextensible circular arch is studied. Again the 58 period of oscillation of the arch is sought. The arch selected is a sixty degree, pin-ended arch having physical characteris- tics similar to aluminum, and with a ratio of axial rigidity to bending rigidity of 3000 in-Z. In order to be certain of excitation in the same mode as for the problem having known analytic solution, the top of the arch was restricted so no motion can occur. Table 3.9 lists the results of the computer computations with the percent deviation from the known period of oscillation--Den Hartog [21]--shown for each mass point that moves. The mass points are numbered so that the mass point at the tOp of the arch is number four, as illustrated in Figure 3.9. Table 3.9.--Inextensive circular arch of seven mass points. Fundamental period of oscillation for each mass point that moves and deviation from linear con- tinuum results. Computed Period Percent Mass Point (Seconds) Deviation 1 .00506 -6.65 2 .00583 7.93 3 .00566 4.43 The results achieved from the beam and circular arch problem compare favorably with the known, linearized solutions of these problems. Also, it is evident from.Table 3.7 that as the number of mass points increases, the accuracy of the results also increases. Thus, as the model approaches 59 Figure 3.9.--Inextensible arch of seven mass points. a continuous distribution of mass along the bar or arch, the solution tends to the continuum solution. As the motion amplitudes are allowed to increase, greater deviations from these linear results would be expected to occur since then the nonlinear effects which are included in the com- puter program would become more dominant. 3.4.3 Application to a Parabolic Arch The parabolic arch problem is further sub-classified into the following categories, 1. Single mass point, single degree of freedom. 2. Single mass point, two degrees of freedom. 60 3. Three mass points, three and six degrees of freedom. 4. Five mass points, five degrees of freedom. Since the present analysis is not intended to be limited to shallow arches, no definition of "shallow" is necessary here. The parabolic arch used was selected arbi- trarily, and its center line lies along the curve 2 4 x ) y = -9- (x - i” . (3.25) Thus we have an arch of L units of horizontal span. with a rise of L/9 units at the peak of the arch. The ori- gin of coordinates is taken to be the left, pinned, end of the arch and a Cartesian coordinate system is used through— out computation. While a symmetric arch has been selected,. the results of computation in category two show that such a condition is not necessary to obtain results.r Physically the material of the arch was again similar to aluminum, and had a ratio of axial stiffness to bending stiffness of 3000.in-2. The time period over which data was collected was that time period sufficient to determine if snap buckling had occurred. The criterhmlused was that if all mass points moved to a position below the line between the end points of the arch, snap failure was assumed. In all attempts where this criterionwas not met, the arch returned to a position and assumed velocity components indicating a trend to re- cover from the induced motions. When the criterfizlwas met, 61 snap-through was evident from phase plane type observations of the displacement vs. motion of the mass points. The time interval of computations was not the same for all problems in the categories listed. The largest possible time increment was always used, to a multiplicative factor of ten. This was done.to conserve computation time as much as possible. The time intervals of computation ranged from .0001 seconds for the single mass point, single degree of freedom problem to .000001 seconds for the three and five mass point arch problems. The computer program is so devised that when convergence tests fail, the time inter- val of computation is reduced by a factor of ten. This action has not produced a discernible discontinuity in any problem output. Further, attempts to refine data computa- tions by using a finer time grid have not led to the con- clusion that the problem solutions are time interval sensitive, so long as sufficiently small intervals are selected so that the solutions do converge.- However, it is to be expected that some variations will occur since the assumption in the solution is that the acceleration is linear over the time interval, and any variation from this assumption will be reflected in the computer output. Finally, Newmark [35] gives a proof of the accuracy of the method based on the time interval selection as shown in Equation (3.24) . Within category 1, the results from five problems are reported here. These are, 62 a. critical initial velocity (pure impulse) with a sensitivity of one inch per second. b. Critical force, applied suddenly and left on with a sensitivity of one pound. c. Critical impulsive force, sensitivity of one pound over the time interval of application. d. Maximum force applied over forty time increments, sensitivity of one pound. e. Maximum force applied over eighfiztime increments with a one pound sensitivity. These problems were selected to help show the scope and resolution of the problems that are possible with this program. Computing time is low for each of the first three problems, averaging approximately four seconds per problem using the CDC 3600 computer, and approximately five minutes per problem using the IBM 1620 Ercomputer. As more mass points are added, the computation time increases rapidly. The last two problems in category 1 were included as approximation of a static loading situation for this single degree of freedom problem. The force was increased in steps and the resultant motion was observed. The results of these computations appear in Tables 3.10 through 3.14 and in Figures 3.10 through 3.14. 63 Table 3.10.--One degree of freedom arch-~pure impulse. Time v v v v (Seconds) (Inches); (In/Sec): (Inches) (In/Sec) .0002 -.688121 -2040.03 - .688330 -204l.69 .0004 -.9l4688 — 530.32 - .915445 - 535.97 .0006 -.971352 - 123.46 - .974338 - 144.14 .0008 -.981660 4.38 - .993656 - 79.05 .0010 -.968972 142.18 -l.017283 - 193.92 .0012 -.904556 600.78 -1.098872 - 747.18 60 = 4417 in/sec 60 = 4418 in/sec Table 3.11.-~0ne degree of freedom arch--suddenly applied force. Time v 6 v v (Seconds) (Inches) (In/Sec) (Inches) (In/Sec) .0002 -.41266 -3546.30 - .4127? -3547.44 .0004 -.83807 -1oss.91 - .83852 -1058.81 .0006 -.95242 - 258.82 - .95415 - 270.67 .0008 -.97800 - 36.52 - .98490 - 84.49 .0010 -.97225 102.89 -l.00004 - 90.39 P = 3877 lbs. P = 3878 lbs. 64 Table 3.12.--One degree of freedom arch--impulsive force. Time v v v v (Seconds) (Inches) (In/Sec) (Inches) (In/Sec) .0002 -.12467 '1388.13 -.12480 -1389.55 .0004 -.40969 -1184.87 -.41021 -1187.64 .0006 —.57293 - 500.07 -.57441 - 507.98 .0008 -.63445 - 164.83 -.63888 - 189.47 .0010 -.65482 - 21.54 -.66506 - 97.85 .0012 -.64487 93.96 -.68675 - 141.98 F = 1035 lbs. F = 1036 lbs. Table 3.13.--One degree of freedom arch--1oad applied over 40 time cyles. No. of v v v v Cycles . (Inches) (In/Sec) (Inches) (In/Sec) -.01528 - 92.60 -.01529 - 92.67 -.03571 — 11.31 —.03574 — 11.32 12 -.06113 - 94.56 -.06118 - 94.65 16 -.07945 - 28.11 -.07952 - 28.13 20 -.11167 - 100.26 -.lll77 - 100.39 24 -.l3425 - 40.85 -.13439 - 40.89 28 -.17071 - 124.64 -.17087 - 124.82 40 -.32237 - 212.76 -.32291 - 213.78 48 -.46205 - 68.79 -.46601 - 80.46 52 -.47359 6.024 -.48684 - 34.53 60 -.39653 213.15 -.53577 - 159.38 'R = 1326 lbs. R = 1327 lbs. Table 3.14.-40ne degree of freedom 65 80 time cycles. arch--1oad applied over No. of v 6 v 6 CyCIes (Inches) (In/Sec) (Inches) (In/Sec) 10 -.02221 - 41.06 -.02223 - 41.06 20 -.05331 - 27.34 -.05335 - 27.39 30 —.07904 - 12.55 -.07911 - 12.55 40 -.11074 - 45.10 -.11084 - 45.10 50 -.15072 - 60.13 -.15086 - 60.21 60 -.19673 - 64.93 -.l9694 - 65.12 70 —.25349 - 84.19 -.25381 - 84.48 80 -.33943 -135.50 -.34015 -136.09 90 -.45843 - 41.07 -.46408 - 56.54 95 -.46586 7.79 -.48891 - 58.49 100 -.44865 68.27 -.55148 - 58.95 R 1336 lbs. R 1337 lbs. 5000 o g 1 \ O c . '11 . *- = 4000 > >1 0\ u u - .8 o ,2 g 1000._ v = 4418 )3 1 1 L ‘ l . 1' ' I l I . r 0 .2 .4 . .8 '1.0 Displacement 6 (inch) 0 = 4417 A Figure 3.10.--Sing1e degree of freedom arch--pure impulse loading. 66 4000~ P = 3890 P = 3878 P = 3600 3000" ’6 m m \. c I: 2000" >1 4) -a o o r-I m > 1000" L . 1 1 1 1 1 1 L l 1 r r [ I 1 I I I T l _I_ .2 .4 .6 .(1 .8 . 1.2 Displacement (inches) Figure 3.11.--Single degree of freedom arch--impulsive force 2000“ = .8 F 1040e5 0 _ ‘3 F 1000_o c :1 34 1000*? -a o o '3 > i i i ; 1 i i = 4 = '2 04 $ .6 08 1.0 Displacement (inches) , Figure 3.12.--Sing1e degree of freedom arch--suddenly applied force. 67 .7 TE .33 o .5 R= 1332 C 'H 1320 V R: u .A C: 0’ a 8 m .3 H C: m -H O .1 : : .L .L 20 4o 60 80 No. of Time Increments Figure 3.13.--Single degree of freedom arch, force increas- ing to maximum over 40 time increments. 06-0- 3? 0 .C‘. U c :1 .4.. u G i U .3 a. ~21" m -H D 4 4v 5 i i : 20 40 60 80 100 No. of Time Increments Figure 3.14r-Single degree of.freedom arch, force increasing to maximum over 80 time increments. 68 The second parabolic arch problem considered was the two degree of freedom, single mass point problem. This prob- lem was considered since it demonstrates the unsymmetrical capabilities of the program.‘ Three different problems are” reported here. These are, 1. Bars of different length, but having the same stiff- ness characteristics. The bar lengths used, compared to the standard problem, were 21/2 = 0.89436 and 22/1 = 1.10607. 2. The bars are of equal length, but have different stiffness characteristics. The axial stiffness ratios used were kl/k = 1.0000 and kz/k = 2.0000. 3. The assembly is symmetric, but a horizontal, non- symmetric force of 200 pounds is applied. In each problem, motion was initiated by providing a pure velocity impulse in the vertical direction for the mass point. Nonsymmetrical behavior resulted in a manner which might be expected in each case, as illustrated in Figure 3.15. The results from these problems appear in Tables 3.15 through 3.18, along with a symmetrical problem output where the initial impulse was the same as that for the non- symmetrical problems. This symmetrical problem is included for comparison purposes, as a control problem. 69 v . .0007 v .0004 0.05" ‘ 0.05“ 0 4- 0.00 0 -- 0.00 .0006 _0.10-_ “0002... ‘- i i i f f } -0.002 0 0.002 ‘ -0.005 0 0.005 u u a) Unequal bar length b) Unequal bar stiffness v V 0020"" .0007 0 d— 0000 .0006 0 -. 0.00 .0001 -°°°5 -0.10.. .0015 .0002 .00022 -0.20.- 1 1 L l l l I v -0.002 0 0.002 u 0 u c) Applied horizontal force d) Control problem ’ Real time values shown adjacent to points. Figure 3.15.--Unsymmetrica1 behavior of single degree of freedom arch approximation. Table 3.15.—2Two degrees of freedom arch-~unsymmetrical mass-- 70 initial velocity of lOOO'ifi/sec. Time u u v v (Seconds) (Inches) (In/Sec) (Inches) (In/Sec) .00010 .00237 10.11 -.08336 -528.81 .00020 .00207 15.02 -.09060 394.94 .00030 -.00012 -20.00 -.01557 986.08 .00040 -.00183 -38.33 .07290 598.17 .00050 -.00144 9.41 .07949 -482.95 .00060 .00016 46.49 —.00349 —998.24 .00070 .00143 14.35 -.08517 -449.56 .00079 .00221 -29.72 -.09297 335.55 Table 3.16.--Control problem--one degree of freedom-—initial velocity 1000 in/sec. Time v v (Seconds) (Inches) (In/Sec) .0001 -.09391 -817.3 .0002 -.15438 -362.2 .0003 -.16298 193.7 .0005 -.03079 983.2 .0006 .06648 885.8 .0007 .13246 368.0 .0008 .13358 -344.8 .0010 -.02743 -987.0 71 Table 3.17.--Two degrees of freedom arch--bars of unequal stiffness-initial velocity of 1000 in/sec. Time u;z u v v (Seconds) (Inches) (In/Sec) (Inches) (In/Sec) .0001 -.00732 2.85 -.08197 -504.51 .0002 -.00748 -10.93 -.08500 448.68 .0003 -.00040 - .59 -.00637 997.60 .0004 .00759 13.32 .07670 488.13 .0005 .00670 - .55 .06960 -615.47 .0006 -.00189 - .96 -.01915 -978.02 .0007 -.00802 - 5.55 -.09005 -327.89 .0008 -.00649 -15.80 -.07474 614.09 Table 3.18.—-Two degrees of freedom arch-—horizontal force of 200 pounds applied to the mass point--initia1 velocity of 1000 in/sec. Time u . .. u . v 6 (Seconds) (Inches) (In/Sec) (Inches) (In/Sec) .00010 -.00202 54.41 -.08665 -618.54 .00015 .00000 - 7.22 -.10832 -236.62 .00018 -.00117 -59.65 -.11166 14.74 .00020 -.00230 -46.47 -.10968 182.74 .00022 -.00281 — 1.15 -.10438 346.16 72 Categories 3 and 4, the three and five mass point problems, are included to show the influence of adding mass points for the parabolic arch problem, and also to show the effect of neglecting the horizontal motion of mass points in the analysis of an arch. In all of these arch problems, care has been exercised to insure that the same arch is being approximated. The only initial condition considered for these problems was a vertical velocity impulse of the cen- tral point of the arch. Since no reference to this problem is available in the literature, no conclusions about convergence to a known solution are possible. The computer program, whose listing appears in Appendix B, was devised to seek out the buckling impulse to within one inch per second of the initial velocity of the peak mass point. If an impulse velocity less than the criti- cal value is used as initial input data, the program will cause this value to increase by increments selected by the programmer until buckling occurs. Buckling is defined in the program as the time when all mass points simultaneously fall to a position below a straight line connecting the ends of the arch. Stability is determined by observation of the output, and is readily detected by noting the return of the mass points to their original positions. The program is so constructed that the programmer must decide the length of run time sufficient to insure that if the arch has not 73 buckled, it will be stable. A little experience with the program in the form of several trial runs readily establishes this time factor, but some care is necessary to avoid false indicators. No completely satisfactory scheme other than this one has been found. No attempt was made to try more than five mass points for an arch approximation. The computer time quickly becomes quite large beyond a three mass point approximation., In fact, the time required for the IBM 1620 II computer to load, execute and print the output for the three mass point problem' as given here, starting with an initial input impulse velocity of 16,000 inches per second and working to completion of the problem with a sensitivity of one inch per second, was ap- proximately six hours. During this time, twelve values of impulse velocity were examined by the computer. The five mass point problem required even more time, since more mass points were involved and also since it was necessary to ex- tend the run time for each value to insure stability, com- pared to the three mass point problem. It was not possible to obtain a sufficiently large block of computer time to permit the five mass point problem to be run on fully auto- matic as described for the three mass point problem. Of special interest in the three mass point parabolic arch problem is the influence of restricting the motion of mass points to the vertical direction, which is the problem usually found in the literature on snap buckling of shallow 74 arches. The results show quite clearly that this additional constraint on each mass point has a marked influence on the oscillation mode and frequency, as well as the critical value: of the buckling impulse. For example,.the work of Humphreys [26] is explicitly limited so that only motion perpendicular to the straight line connecting the ends of the arch is con- sidered. It would appear that this constraint may somehow need to be included in any definition of shallowness of an arch. It is possible that this constraint provides the ori- terionupon which shallowness might be judged. For example, it is possible that any arch for which inclusion of this constraint would not appreciably alter the buckling load may be considered shallow, otherwise not.- In any event, the program used here is capable of providing results without consideration of any definition of shallowness of an arch. Figures 3.16 through 3.19 show the results of the _ problem solution in the vicinity of the critical values for the initial velocity imparted to the peak mass point of the arch. Only the vertical motion of each mass point is plotted. The figures plotted represent all the output from the computer program, but do not include all the points that were calculated by the computer. To conserve computer time, and to avoid the problem of becoming inundated with data, the computer was programmed to print out data only after one hundred time increments had been processed. Thus each set of points plotted represents one hundred sets of calculated 75 2.01. vo=l6229 v0 = 16230 1 ml 5 m1 106‘” ‘g o a 2 " o A 5 1.2+ ‘ 4.} a m E 8 I I 3 0.8“ 3‘ ° A -H o 0.4-- . , : ; \ 4 60 80 \ 100 10"5 secs. Y 2 - m m1 2 m1 i, % I l X 2 4 6 8 10 Figure 3.16.--Three mass point arch, restricted motion. Displacement (inches) 76 _0I4d (Displacements of m to plot here.) 1 Time x 10_5 secs. for t < 50 x 10"5 are too small Figure 3.17.—-Five m m L 10 mass point arch, restricted motion. X 77 Displacement (inches) l I A 1 j l I l I J. l l db fih - Id Time x 10"6 secs. v = 1900 v = 2000 v = 2100 m1 % “‘1 43 m1 ----A m2 0 m2 —5 m2 ----——0 Figure 3.18.--Three mass point arch, mass points free. 78 1.5x10'4 a) Initial Velocity 16000 in/sec. 4 T = 2x10- b) Initial Velocity 18000 in/sec. Figure 3.19.--Disp1acement response of three and five mass point arch approximations. 79 data. Therefore, the seemingly large jumps in motion made between points did not in fact actually exist during calcula- tion. This time increment of print out is a variable that the programmer may select for each problem run. . In all of the problems reported here, the same, or similar statements can be made. The output of the computer, and hence the points plotted do not represent all the data points calculated. The frequency of print out was varied to suit the problem and was selected, by trial, so that a bal- ance between having a sufficiently large number of points to plot and having a large amount of unnecessary data was achieved. Table 3.19 shows the gnumerical output for selected time values for the three mass point problem with the mass points free to move horizontally and vertically, for one value of the initial velocity.‘ A study of a physically symmetrical, symmetrically loaded arch was performed previously, using the same computer program used to collect these data, where it was not assumed in inserting data to the program that symmetry did exist.,' The mass point at the peak of the arch exhibited horizontal, nonsymmetric motion only to within the ability of the com- puter to approximate zero through the internal errors of computation which occur. These are errors such as trunca- tion, orround-off errors, depending on the particular model of computer used. Additionally, symmetry was evident from 80 mmoo. wm.m~oa vmmmm.l No.5mHI havHH.I mm.Hm~I mwmwo. em.HmH I mmHHQ.HI mm.mmh H¢m~m.| Hm.Oh.I mmmmo. «moo. mm.vm~ I NvaH.HI mm.mmm mmHvM.I mm.mm L amuse. emoo. Hm.om ommmH.HI hm.om I mhbvm.l mm.m mambo. omoo. mm.mm~ bommH.HI mm.m~vI mmmam.I m¢.om I hammo. maoo. Hw.ama mmmmo.HI hm.mmmI hvhm~.I vH.mNHI mwamo. mace. mo.a>~ I emmmm.| mm.vmmI mvme.I m~.mva mmowo. mooo. mm.vvaI moqmm.I nm.o>m mvmmo.I mm.ov I mmmoo. vooo. on.avaI hmmmm.I o~.hma vmqvm.l mm.m I Hmvao.I mooo. om.mmwal ~¢hhH.I m>.m~¢ mmqao.l mm.mmHI mHoHo.I Hooo. Aomm\ch Ammnoch Aomm\cHV Amococh Aoom\ch Amococh Amosoommv mé N> H¢ H> Ha do mafia LI .Ucoomm mom monocfl coma mo ucfl0n mmmE mo» um mmaamfifl HMOfluHm>IInoum Eocmmum mo mmummp xwm .ucwom mmmfi ochBII.mH.m manna 81 the observation of the motion of symmetrically placed mass points. In all subsequent work reported here, to conserve computation time, symmetry conditions were inserted at the peak mass point of the arch.. The computer program does not rely on this compatibility, however, and unsymmetrical prob- lems can be worked using it. Therefore, in this problem the top mass point did not move horizontally even though it was free to do so within the compatibility limits described. The other mass points did move horizontally, with motion amplitudes of from one-half to one order of magnitude smaller than the vertical amplitudes. 3.4.4 Application to a Column Equations (2.20) have been used to determine the buckling load for a column. The specific model investi- gated has physical characteristics similar to aluminum and is approximated by a two bar, two mass point model of non- constant mass distribution. The mass distribution used was that of Figure 3.6. The principal difficulty encountered with this prob- lem was, as in any dynamically programmed problems, the simulation of static loading. To approximate static loads, the column load was incremented using a time dependent step function. After completing the calculations for each time increment the load was increased up to a predetermined maxi- mum value. After full load was achieved, the motion of the 82 column was monitored until it was evident that the column had begun to oscillate. The selection of a total time value over which full load was applied was critical. Ex- perience with the arch problem reported in section 3.4.3 indicated a starting point, and longer times were selected until the final result that is presented here. It is prob- able that even longer time intervals for the load applica- tion would permit larger loads to be applied without buckling being indicated. This is the primary difficulty of static load simulation using the Newmark scheme of iteration. The establishment of a load cycle time for a particular column model by comparison of the results obtained to a known solu- tion should provide a satisfactory basis for investigation of additional load configurations for that column model, and it should also aid in the establishment of a confidence level for the results obtained from further investigations of that particular column. For the particular problem selected as the final problem in this project, as illustrated in Figure 3.6, the following results were obtained. When the load was forced to maintain its original direction as the column displaced, the ratio of applied load to the known Euler load for the continuum column was in the range of 1.0117 < P/PCr < 1.0376 . (3.26) At the upper value of the range, the column collapsed before full load was applied. At the lower value of the range, the 83 column clearly oscillated when its motion was observed for a time increment of 2-1/2 times the time increment over which the load was applied. During this additional time span, more than one full quarter of the period of oscillation of the loaded column was observed. When the same column was subjected to a nonconserva- tive load, specifically the follower load problem, the range of load ratios was 1.0461 < P/Pcr < 1.0620 . (3.27) Again at the upper end of this range, the column collapsed before full load application, and again, more than one full quarter of the oscillation period was monitored before the problem was terminated when using the smaller value of load ratio given in equation (3.27). The load was applied in time increments of .007143 times the fundamental period of oscillation of the continuum problem--0.001 seconds of real time—-and was fully applied at .7143 times the fundamental period, for the conservative load problem. Correspondingly, the values were .007143 and 1.2664 times the fundamental period of the continuum column for the follower load problem. CHAPTER IV DISCUSSION, CONCLUSIONS, AND FUTURE STUDIES 4.1‘ The Problem Summary The elastic behavior of circular rings, beams, columns, and arches has been investigated using a lumped mass physical model to simulate the continuum problem. Both static and dynamic response to a variety of loads hamebeen considered. The influence of bending and axial strain have been included in the analysis, separately and together. Nonconservative as well as conservative load- ing conditions have been considered. Unsymmetric capa- bilities of the model have been demonstrated. The problem of the buckling of an arch has been considered where the arch geometry has not been limited by shallowness or other Special geometric constraints. Only the constraint of elastic action has been employed, and that by implication, since no bounds for elastic action have been established. Attempts to include the effects of transverse shear and rotary inertia in the beam and arch problems have not proven successful. This lack of success is attributed to the nature of the iteration scheme used in the computer program, specifically the Newmark scheme described in Chapter III. 84 85 Thus, the problem considered here has been the determination of the fundamental mode response to static and dynamic loads of the undamped, one dimensional con- tinuum elastic bar exclusive of transverse shear and rotary inertia effects. 4.2 Results of the Analysis The physical and mathematical models used to in- vestigate this problem have been shown to be reasonable models. Within each category of problem considered it has been shown that the models used do converge toward the known solutions for the continuum problem being ap- proximated, when such solutions are available. The re- sults obtained for the static problems are in.each case in good agreement with the known continuum problem solu- tions. For the dynamic problems considered, there wa8~ good agreement for the period of oscillation of a beam as mass points were added to the physical model, as well as obtaining good agreement for the period of oscillation of the inextensible circular arch and the buckling load of the column. No results are available in the literature to confirm the validity of the model for the buckling of a parabolic arch excited by impulse load. However, the articles by Eppink and Veletsos [24] and Lind [32] lend validity to the method used. Eppink and Veletsos con- sidered the radial motion of a circular arch, using 86 Newmark's method, including axial deformations and bend- ing moments, and compared the solution obtained to the solution using a modal method of analysis. Their problem solution does not have the flexibility of the one used here, in that the one used here is not restricted to a circular arch and uses all of the nonlinear terms of the equations of the problem, and does not terminate an ap- proximating series. Lind shows the application of the method of New- mark, in a very general way, to arches and rings of arbi- trary shape, and shows the method to be capable of accurate, yet practical solutions of problems otherwise not readily yielding to solution using known mathematical methods. It has further been established here that the Newmark scheme of iteration is applicable to the various types of problems considered and is not restricted to linear analysis. The models used have not been shown to provide a means for including transverse shear and rotary inertia effects, nor has it been shown possible to include the effects of damping in the dynamic problems considered. 4.3 Conclusions The models used to approximate the beam, column, and circular ring continuum problems can be expected to yield accurate results for the period of oscillation of 87 a beam, the bifurcation load values for a column includ- ing the effects of axial strain, and the buckling pressure for a circular ring. The period of oscillation of a pin-ended beam, using three point extrapolation and the results of 7, 9, and 11 mass point models was in error by 0.57 percent com- pared to the continuum linear solution. For the column including axial strain effects, the values of the critical load ratios at the bifurcation point deviated from the continuum solution by 0.036 per- cent and 0.027 percent when three point extrapolation of the 5, 6, and 7 bar approximations was used. Results for the circular ring problem showed er— rors of 1.002 percent and 1.120 percent for central pres- sure loading and 0.750 percent and 0.713 percent for normal pressure loading when the 6, 8, and 10 bar approximations were extrapolated. Based upon these results, it may be concluded that the models used are accurate and that it is not necessary to use more than twelve bars, or eleven nodes, when three point extrapolation is employed to produce results to ap- proximately 1 percent error for the types of problem con- sidered here, and when using the techniques described in Chapter III. When the incompressible circular arch problem re- sults are compared to comparable results for the problems 88 listed above, the errors are of similar magnitude. For example, the largest error in the period of oscillation for one mass point of the seven mass point inextensible circular arch was 7.93 percent. This compares to the error in the period of oscillation of a beam with seven mass points of 7.98 percent. Thus the incompressible circular arch problem results can be expected to converge in a manner similar to that of the beam results. 4.4 Future Studies Several additional problems are suggested by the work done here. These include the determination of trans- verse shear and rotary inertia effects using the models develOped here, however, using an as yet undetermined iteration scheme for the computer solutions. The influ- ence of velocity dependent damping is a problem which might be investigated using the iteration scheme employed here. The work presented here also suggests that an at- tempt at approximation of the errors introduced by finite difference approximations of differential equations might be undertaken. In particular, the solution of equations (2.2) could be compared to the solutions of certain column buckling problems, for example, those given in references [19] and [34]. These comparisons might then give insight into finite difference errors in view of equations (2.26) and (2.27). 89 Different methods of mass distribution might be investigated to attempt to determine that method which would provide most rapid convergence for the problems considered here, provided such a mass distribution method does exist. The extension of the circular ring problem to in— clude the effect of axial strain might be attempted by utilization of the physical and mathematical models in- cluded in this presentation. The physical model presented lends itself well to the solution of problems where the bar cross-sections are not constant, or where the bar is not of constant material properties. Finally, a study of the influence of time incre— ment on the rate of convergence of computer determined solutions toward known problem solutions could be useful in establishing an optimum time increment vs. number of mass points curve for several categories of problems. APPENDICES Included in this Appendix section are the flow chart and computer program used in the foregoing analysis, showing the logic sequence used in the computations. The computer program was used to gather the data on which this thesis is based, except for the column prob— lem, where some alterations of the program were necessary to incorporate the constraints peculiar to the column problem, and those problems of Chapter III where static loading is used, unless indicated otherwise. The program was so constructed that the initial conditions for the motion could be altered whenever the structural component was found to buckle. Thus, if the initial conditions are selected incorrectly, the computer changes these conditions and proceeds to seek out the value of the critical param- eter for the problem within the predetermined range of values of the parameter that is considered acceptable. 90 91 A. Computer flow chart. 9:9 INPUT PROBLEM PARAMETERS INPUT PROBLEM CONSTANTS AND INITIAL VALUES [PRINT PROBLEM HEADINgj r V [INITIATE COUNTERS] yes TEST IF SYMMETRICAL ? ADJUST INPUT PARAMETERS FOR LACK OF SYMMETRY ASSUME TRIAL VALUES FOR ACCELERATIONS AT END OF TIME INTERVAL, COMPUTE VELOCITIES AND POSITIONSfi L COMPUTE RAMP LOADING PARAMETER © ® TEST IS PROBLEM Q INEXTENSIBLE? SET BAR DEFORMATIONS COMPUTE BAR TO ZERO DEFORMATIONS 7 [COMPUTE BAR ROTATIONsl 7 q—LEST IF SYMMETRICAL @ SET SYMMETRY CONDITIONS SET LACK OF SYMMETRY CONDITIONS 93 SET BOUNDARY CONDITIONS COMPUTE AXIA FORCES AND MOMENTS I S TEST IF U MOTION @ Ye _ IS DESIRED - SET U. = 0‘ .1 ATIONS FOR END ' [OF TIME INTERVAL_ COMPUTE ACCELER-J COMPUTE DIFFERENCES BETWEEN TRIAL VALUES AND IMPROVED VALUES OF ACCELERA- TIONS .—-(IS M COUNTER > 2? ® 1 {:59 TEST FOR CONVERGENCE) (:5) ixflg, to(:> to F 92 TEST IF ACCELERATION DIFFERENCES EXCEED TOLERANCES OMPUTE CURRENT TIM RESET INITIAL CONDITIONS TO LAST ASSUMED VALUES OF ACCELERATIONS, VELOCITIES AND DISPLACEMENTS E TEST IF TIME TO PRINT)—-® COMPUTE NEXT TIME TO PRINT TEST IF ACCELERATION DIFFERENCES ARE CHANGING PRINT ACCELERATION DIFFERENCES TEST IF RATE OF CON- Q VERGENCE IS TOO SLOW t TEST IF MCOUNT > 2 Q to PRINT 'END OF PROB- LEM' to TEST IF PROBLEM HAS @ RUN FOR PRESCRIBED Q I TIME (RICREASE TOLERANCQ to (I) COMPUTE TRIAL VALUES FOR DISPLACEMENT, PRINT TOLERANCE VELOCITY AND ACCELER- ATION L to to® © PRINT DISPLACEMENTS VELOCITIES AND ACCELERATIONS USING IMPROVED ACCELERATIONS COMPUTE VELOCITIES AND DISPLACEMENTS INCREMENT MCOUNT to® C? 93 PRINT 'ACCELERATION DIFFERENCES INCREASE' ® IS MMCOUNT < 3? REDUCE TIME INCREMENT I! INCREMENT MM COUNT I @ @ PRINT 'END OF PROBLEM' PRINT PROBLEM' 'END OF 1 Ix>() COMPUTE NEW TRIAL VALUES FOR ACCELERATION, VELOCITY AND DISPLACEMENT FOR END 'OF NEW TIME INTERVAL 1 to© O 9‘9 TEST IF INITIAL CONDITIONS EQUAL MAXIMUM VALUES -u£() TEST IF INPUT INITIAL CONDITIONS INCREMENT DECREASE INITIAL CONDITIONS INCREMENT COMPUTE REDUCED INITIAL CONDITIONS PRINT NEw INITIAL CONDITIONS PRINT 'END OF PROBLEM' @ to. COMPUTE INCREASED INITIAL CONDITIONS 7 PRINT NEW INITIAL CONDITIONS tolfl) 94 B. Computer Program Listing Listed below are the variable name definitions used in the computer program which follows. Also included is a brief resume of the interpretation of statements the computer may print out in the course of the solution of a given problem. The program listing is designed for use on an IBM 1620 II computer which uses a modified FORTRAN II compiler. This listing will also function, with very minor changes, on a machine which uses a FORTRAN IV or FORTRAN V compiler. The program is an elementary one. No attempt has been made to produce a program with a high level of sophis- tication. However, since the program is capable of solving a great variety of problem types, it may be difficult to follow its logic throughout. Some careful study and trial is needed to fully utilize its capabilities. N Number of bars in the segmented structure. J Control for symmetry, equals 1 if symmetric, equals 2 if unsymmetric. Always equals 2 for two bar, one mass point problem. LL Equals 1 if mass points have two degrees of freedom, equals 2 if mass points move in vertical direction only. KA Inextensibility control. Equals 1 if bars are extensible, equals 2 otherwise. LA Ramp load parameter. Equals 1 if load is suddenly applied, equals 2 if load is gradu- ally applied, equals 3 if load is gradually removed after sudden application. PT PPT VELI DELUAC DELVAC TOLU TOLV T CC VX VY AX AY PHI XL ETA CAY 95 Number of time intervals over which load is to be applied. Number of time intervals load is left on after full application. Velocity impulse increment. Acceleration increment, U direction. Acceleration increment, V direction. Acceleration tolerance, U direction. Acceleration tolerance, V direction. Time increment. Number of time increments program is to run. Number of computation intervals between print out cycles. Displacement in U direction. Displacement in V direction. Velocity in U direction. Velocity in V direction. Acceleration in U direction. Acceleration in V direction. Initial angle of bar from horizontal. Initial bar length. Applied horizontal force. Applied vertical force. Applied distributed pressure. Initial angle of distributed pressure vector from the vertical. Mass per unit volume of the bar. Bending stiffness ratio for adjacent bars. 96 ALPHA Ratio of axial stiffness to bending stiff- ness of a bar. YM Young's modulus of elasticity. WT Weight of mass point. VC Height of mass point above straight line between ends of the arch. The program is constructed so that the maximum velocity impulse acceptable is 29000. If a different value is desired, statement number 883 must be altered. If the statement ACCELERATION DIFFERENCES INCREASE is printed out, the time increment will decrease by a fac- tor of ten and the problem will continue. After this mes- sage prints out three successive times, the program termi- nates. If the tolerance chosen is too small for signifi- cance in the computer comparison statement, the computer increases the tolerance by a factor of ten, prints the new value of TOLV and continues. If this trouble persists more than three times in any one problem, the program terminates. If the acceleration differences UDIF and VDIF can- not be improved by successive iterations, the computer prints the values obtained and continues, assuming conver- gence has been achieved, whether the tolerance is satisfied or not. 1 2 800 100 101 103 301 302 801 97 Prpgram Listing DYNAMIC ANALYSIS OF BEAMS AND ARCHES DIMENSION X(ll),VX(ll),AX(11),VY(11),Y(ll),AY(ll), 2UVEL(ll),VVEL(11),U(1l),V(ll),DEL(11),A(11),B(11), 3C(ll),D(ll),PSI(ll),PHI(11),UACC(11),THETA(11),XMV 4(11),YY(11),VC(1l),VACC(1l),WT(ll),CAY(11),XN(11), 58(11),YM(11),ALPHA(11) DIMENSION E(11),F(11),UDIF(30,11»VDIF(30,11), 2H(1l),ETA(11),XL(11),XX(11) READ 100,N,J,LL,KA,LA,PT,PPT,VELI PRINT 205 READ 101,DELUAC,DELVAC,TOLU,TOLV,T,CC,P CCT=CC*T DO 3 I =1,N READ 103,X(I),VX(I),AX(I),Y(I),VY(I),AY(I) READ 103,PHI(I),XL(I),E(I),F(I),H(I),ETA(I) READ 103,XMV(I),CAY(I),ALPHA(I),YM(I),WT(I),VC(I) VEL=VY(N) ICOUNT=0 FORMAT (515,3F10.0) FORMAT (7F10.0) FORMAT (6F10.0) DO 4 I=1,N IF (ABSF(PHI(I))-.1E-5) 301,301,302 C(I)=1. D(I)=0. GO To 4 C(I)=COSF(PHI(I)) D(I)=SINF(PHI(I)) CONTINUE MCOUNT=0 VY(N)=VEL SET COUNTERS AND CONTROLS TO INITIAL VALUES TT=0. AT=.5*T PT=PT*T PPT=PPT*T+-PT TTT=P*T MM:1 GO TO (8,5),J N=N-l GO TO 8 DELVAC=.1*VACC(1) DELUAC=.1*UACC(1) M=o COMPUTE APPROXIMATE DISPLACEMENTS DO 7 I=1,N VACC(I)=AY(I)4-DELVAC 14 15 16 18 19 20 21 211 22 23 24 25 251 252 253 26 C 36 37 39 98 VVEL(I)=VY(I)-I.5*(AY(I)I-VACC(I))*T V(I)=Y(I)-+(VY(I)-+((AY(I)-+.5*VACC(I))/3.)*T)*T UACC(I)=AX(I)-+DELUAC ‘ UVEL(I)=VX(I) + .5*(AX(I) +UACC(I))*T U(I)=X(I)+—(VX(I)-I((AX(I)-+.5*UACC(I))/3.)*T)*T COMPUTE RAMP FUNCTION PARAMETER GO TO (16,14,19),LA IF (AT-PT) 15,16,16 BBT=TT/PT GO TO 21 IF (AT-PPT) 17,18,18 BBT=1. GO TO 21 BBT=0. GO TO 21 BBT=1.-TT/PT IF (BBT) 20,21,21 BBT=0. M=M-+1 COMPUTE BAR EXTENSION AND ROTATION A(1)=U(1) 'B(l)=V(1) IF (N-l) 1,24,22 DO 23 I=2,N A(I) = U(I) - U(I-1) B(I)=V(I)-V(I-1) GO TO (25,24),J N=N+ 1 A(N) = -U(N-l) B(N)= avIN-I) DO 26 I=1,N GO TO (252,251),KA DEL(I)=0. GO TO 253 DEL(I) = -XL(I) + SQRTF((XL(I)*C(I) - A(I))**2 + 2 (XL(I) *D(I) - B(I) ) **2) PSI(I)=ATANF((B(I)*C(I)-A(I)*D(I))/ 2 (XL(I)-A(I) *C(I)-B(I) *D(I) )) CONTINUE COMPUTE SHEAR AND NORMAL FORCES DO 36 I=1,N THETA(I)=PHI(I)-PSI(I) XN(1)=ALPHA(1) *DEL(1) S(l)=(PSI(1) -PSI(2) + 2.*CAY(1) *PSI(1) )/(XL(1) + 2DEL (1)) GO TO ( 39,41) ,J PSI (N+ l) = -PSI (N) L=N GO TO 43 99 41 IF (CAY(N)) 1,42,39 42 S(N)==(PSI(N)-PSI(N-l))/(XL(N)*DEL(N)) XN(N)==ALPHA(N)*DEL(N) L= N-l IF (N-2) 1,54,43 43 DO 44 I =2,L S(I) = (PSI(I)-PSI(I +1)-CAY(I)*(PSI(I-1)-PSI(I))) 2/(XL(I) + DEL(I)) 44 XN(I) = ALPHA(I) *DEL(I) SET COMPATABILITY CONDITIONS 55 GO TO (56,57),J 56 THETA(N + 1) = -THETA(N) XN(N-+1)==XN(N) S(N+l)= -S(N) GO TO 58 57 N= N-l 58 L=M-1 COMPUTE NEw ACCELERATIONS 59 DO 61 I=1,N BBB=ETA(I)-+.5*(PSI(I)-+PSI(I-Il)) IF (ABSF(BBB-.1E-5)) 592,592,591 591 EA=(E(I)4-H(I)*SINF(BBB))*BBT GO TO 593 592 EA=E(I)*BBT 593 XX(I)=(XN(I)*COSF(THETA(I))-XN(I+-1)*COSF(THETA 2(I+-1))I-S(I)*SINF(THETA(I))-S(I+-1)*SINF(THETA 3(I+ 1)))*YM(I)/(ALPHA(I)*XMV(I)*XL(I)*XL(I)) + 4(EA*386.4)/WT(I) UDIF(M,I)=ABSF(UACC(I)-XX(I)) IF (M-2) 61,61,60 60 IF (UDIF(M,I)-UDIF(L,I)) 61,61,72 61 CONTINUE GO To 94 62 LB=N+—1 DO 93 I=1,LB 93 XX(I)=0. 94 DO 71 I=1,N IF (BBT) 64,64,63 63 =ETA(I)4-.5*(PSI(I)-+PSI(I-+1)) IF (ABSF(AAA-.lE-5)) 632,632,631 631 FA=(F(I) + H(I)*COSF(AAA))*BBT GO TO 67 632 FA= (F(I) +H(I))*BBT GO TO 67 64 FA=WT(I) 67 YY(I)=(XN(I)*SINF(THETA(I))-XN(I-+1)*SINF(THETA 2(I4-l)) -S(I)*COSF(THETA(I)) +S(I-+1)*COSF(THETA 3(I+l)))*YM(I)/(ALPHA(I)*XMV(I)*XL(I)*XL(I)+FA* 4386. 4/WT(I) 70 71 72 73 74 742 75 741 751 752 999 76 761 77 78 762 763 764 765 998 100 CHECK CONVERGENCE AND TOLERANCE VDIF(M,I)=ABSF(VACC(I)-YY(I)) IF (M-2) 71,71,70 IF (VDIF(M,I)-VDIF(L,I)) 71,71,72 CONTINUE GO TO (76,64),LL PRINT 204 IF (MM-3) 73,89 ,89 T=.1*T MM=MMI-1 GO TO 6 IF (L-l) 902,742,742 DO 75 I=1,N IF (VDIF(M,I)-TOLV) 75,75,741 CONTINUE GO TO 79 DO 752 I=1,N IF (VDIF(M,I)-VDIF(L,I)) 90,751,90 PRINT 999,VDIF(M,I) CONTINUE FORMAT (30X, llHVDIF(M,I)=PElO.3) GO TO 79 IF (L-l) 902,761,761 DO 78 I=1,N IF (UDIF(M,I)-TOLU) 77,77,762 IF (VDIF(M,I)-TOLV) 78,78,762 CONTINUE GO TO 79 DO 765 I=1,N IF (VDIF(M,I)-VDIF(L,I)) 763,764,763 IF (UDIF(M,I)-UDIF(L,I)) 90,764,90 PRINT 998,VDIF(M,I),UDIF(M,I) CONTINUE FORMAT (30x, llHVDIF(M,I)==PE10.3, llHUDIF(M,I) 2 =PE10.3) C CHECK IF TIME TO PRINT OUTPUT,PRINT OR RECYCLE 79 80 81 82 821 TT=TT+T AT=TTI-.5*T GO TO (84,80),LL DO 81 I=1,N AY(I)=VACC(I) VY(I)=VVEL(I) Y(I)=V(I) IF(TTT-AT) 82,82,6 IF (MM-2) 824,821,823 T=10.*T MM=1 GO TO (86,825),LL 823 824 825 83 807 806 804 805 84 85 86 87 88 881 882 883 891 892 893 894 89 101 T=100.*T MM=1 GO TO (86,825),LL TTT=TTT+-P*T PRINT 201,TT DO 83 I=1,N PRINT 203,Y(I),VY(I),AY(I) GO TO 88 VEL=VELI-VELI ICOUNT=1-+ICOUNT IF (VELI-.5) 89,89,806 DO 804 I=1,N AY(I)=0. AX(I)=0. vx(I)=o. X(I)=0. VY(I)=0. Y(I)=0. PRINT 805,VEL,ICOUNT FORMAT (16X,14HNEW VELOCITY==,PE 12.4,10X,IS) GO TO 801 D0 85 I=1,N AX(I)=UACC(I) AY(I)=VACC(I) VX(I)=UVEL(I) VY(I)=VVEL(I) x(I)=U(I) Y(I)=V(I) IF (TTT-AT) 82,82,6 TTT=TTT+P*T PRINT 201,TT DO 87 I=1,N PRINT 202,X(I),VX(I),AX(I) PRINT 203,Y(I),VY(I),AY(I) DO 881 I=1,N IF (VC(I) -V(I)) 881,882 ,882 CONTINUE GO TO 891 IF (TT-CCT) 6,883,883 IF (VEL-29000.) 807,89,89 IF (VELI-10.) 892,893,893 IF (VELI-1.1) 89,89,894 VEL=VEL-.8*VELI VELI=.2*VELI GO TO 806 VEL=VEL-VELI+-1. VELI=1. GO TO 806 PRINT 206 GO TO 1 90 901 903 904 902 91 911 912 913 201 202 203 204 205 206 300 102 IF (MCOUNT-1) 901,903,903 IF (M-20) 902,902,911 IF (M-20) 902,904,904 M==0 DO 91 I=1,N UACC(I)=XX(I) UVEL(I)=VX(I)-+.5*(AX(I)4-UACC(I))*T U(I)=X(I)-+(VX(I)-I((AX(I)I-.5*UACC(I))/3.)*T)*T VACC(I)=YY(I) VVEL(I)=VY(I) +.5*(AY(I)-rVACC(I))*T V(I)=Y(I)4-(VY(I)+-((AY(I)+ .5*VACC(I))/3.)*T)*T GO TO 21 MCOUNT=MCOUNT+ 1 TOLV=10.*TOLV TOLU=10.*TOLU PRINT 300,TOLV IF (MCOUNT-3) 913,912,912 TOLV=.01*TOLV TOLU=.01*TOLU PRINT 300,TOLV GO TO 89 GO TO (76,74),LL FORMAT (1X,E10.3/) FORMAT (16X,3HU 1PE16.7,4X,1PE16.7,4X,1PE16.7) FORMAT (16X,3HV 1PE16.7,4X,1PE16.7,4X,1PE16.7/) FORMAT (40H ACCELERATION DIFFERENCES INCRE 2ASE/ ) FORMAT (1H1 , 4X , 4HTIME , 12X, DISPLACEMENT , 13X , 8HVELOC ZITY , 10X, lZHACCELERATION) FORMAT (4 0X , 14HEND OF PROBLEM//) FORMAT (4OX,6HTOLV = 1PE10.3) END 10. 11. BIBLIOGRAPHY Books Bolotin, V. V. Dynamic Stability of Elastic Systems. San Francisco: Holden - Day, 1964. Bolotin, V. V. Nonconservative Proplems of the Theory of Elastic Stability. Oxford: Pergamon Press, 1963. Church, A. H. Mechanical Vibrations. New York and London: John Wiley and Sons, Inc., 1957. Davis, H. T. Introduction to Nonlinear Differential and Integrai Equations. New York: Dover PuBIi- cations, Inc., 1962. Den Hartog, J. P. Mechanical Vibrations. 4th ed. New York, Toronto and London: McCraw-Hill Book Company, Inc., 1956. Jacobsen, L. S. and Ayre, R. S. Epgineerinngibra- tions. New York, Toronto and London: McGraw-Hill Book Company, Inc., 1958. Kauderer, H. Nichtlinegre Mechanik. Berlin, Heidel- berg and New York: Springer Verlag, 1958. Langhaar, H. L. Energy Methods in Applied Mechanics. New York and London: John Wiley and Sons, Inc., 1962. Stoker, J. J. Nonlinear Vibrations. New York and London: Interscience PuBlishers, Inc., 1950, Fifth Printing, 1963. Timoshenko, S. P. and Young, D. H. Vibration Problems in Engineeripg. 3rd ed. New York, Princeton, Toronto and London: D. Van Nostrand Company, Inc., 1955, reprinted, Jan. 1966. Timoshenko, S. P. and Gere, J. M. Theory of Elastic Stability. 2nd ed. New York, Toronto and London: McGraw-Hill Book Company, Inc., 1961. 103 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 104 Bailey, H. R. "Dynamic Buckling of Elastic Columns," J. Structural Division, A.S.C.E., Vol. 89, ST 4 (Aug., 1963), 95-114. Ball,R. E. "Dynamic Analysis of Rings by Finite Dif- ferences," J. Engng. Mech. Div., A.S.C.E., Vol. 93, EM 1, (Feb., 1967), 1-10. Beatty, M. F. "Theory of Elastic Stability for In- compressible, Hyperelastic Bodies," Int. J. Solids and Structures, Vol. 3, (Jan., 1967), 23-37. Beatty, M. F. and Hook, D. E. "Some Experiments on the Stability of Circular Rubber Bars Under End Thrust," Int. J. Solids and Structures, Vol. 4, (July, 1968), 623-635. Boresi, A. P. "A Refinement of the Theory of Buckling of Rings Under Uniform Pressure," J. Appl. Mech., Vol. 22, (March, 1955), 95-102. Budiansky, B. "Dynamic Buckling of Elastic Structures: Criteria and Estimates," Paper presented at lst International Conference on Dynamic Stability of Structures, Northwestern University (1965). Budiansky, B. and Hutchinson. "A Survey of Some Buckling'Problems," J. A.I.A.A., Vol. 4, No. 9 (Sept. 1966), 1505-1510. Burgreen, D., "Free Vibrations of a Pin-Ended Column with Constant Distance Between Pin Ends," J. Appl. Mech., Vol. 18 (June, 1951), 135-139. Cowper, G. R. "Shear Coefficients in Timoshenko's Beam Theory," J. Appl. Mech., Vol. 33, No. 2, (June, 1966), 335-340. Den Hartog, J. P. "Vibration of Frames of Electric MaChineS," Trans. ' AOSOMoEo I V01. 50' APM 50-6 (1928). Dengler, M. A., and Goland, M. "Transverse Impact of Long Beams, Including Rotatory Inertia and Shear Effects," Proceedings, lst U.S. Natl. Congress of Appl. Mech. (1951) 179-186. Eisley, J. G. "Nonlinear Deformation of Elastic Beams, Rings, and Strings," Appl. Mech. Reviews, Vol; 16, NO. 9, (Sept, 1963) 677-680. Also, Applied Mechanics Surveys. Washington, D.C.: Spartan Books, 1966, 285-290. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 105 Eppink, R. T. and Veletsos, A. M. "Dynamic Analysis of Circular Elastic Arches," Proceedings, 2nd Conf. on Electronic Computation, A.S.C.E. (1960) 477-502. Herrmann, G. and Bungay, R. W. "On the Stability of Elastic Systems Subjected to Nonconservative Forces," J. Appl. Mech., 31 (Sept., 1964) 435-440. Herrmann, G. and Jong, I. C. "On the Destabilizing Effect of Damping in Nonconservative Elastic Sys- tems," Trans. A.S.M.E. (Sept., 1965) 592-597. Hsu, C. S. "On Dynamic Stability of Elastic Bodies with Prescribed Initial Conditions," International J. of Engng. Sci., Vol. 4, No. 1 (March, 1966) 1-21. Humphreys, J. S. "On Dynamic Snap Buckling of Shallow Arches," J.A.I.A.A., Vol. 4, No. 5 (May, 1966) 878-886. Koiter, W. T. "Elastic Stability and Post-Buckling Behavior," Proceedings, Symposium on Nonlinear Problems, edited by R. E. Langer, University of Wisconsin Press (1963). Koiter, W. T. "Over de Stabiliteit von het Elastich Evenwicht," H. J. Paris (1945). Levy, S. and Kroll, W. D. "Errors Introduced by Finite Space and Time Increments in Dynamic Re— sponse Computation," Proceedings, lst U.S. Natl. Congress of Appl. Mech. (1951) 1-8. Lind, N. C. "Numerical Buckling Analysis of Arches and Rings," The Structural Engineer, Vol. 44, No. 7 (July, 1966) 245-248. Lindberg, H. E. "Impact Buckling of a Thin Bar," J. Appl. Mech., Vol. 32, No. 2 (June, 1965) 315-322. Morris, N. F. "Dynamic Stability of Beam Columns with Fixed Distances Between Supports," J. of the Franklin Inst., Vol. 280, (1965) 163-173. Newmark, N. M. "A Method of Computation for Structural Dynamics," J. Engng. Mech. Div., A.S.C.E., Vol. 85 (July, 1959) 67—94. 36. 37. 38. 39. 40. 41. 42. 43. 106 Schreyer, H. L. and Masur, E. F. "Buckling of Shallow Arches," J. Engn. Mech. Div., A.S.C.E., Vol. 92, EM 4 (Aug., 1966) 1-17. Seames, A. E. and Conway, H. D. "A Numerical Pro- cedure for Calculating the Large Deflections of Straight and Curved Beams," J. Appl. Mech., Vol. 24 (June, 1957) 289-294. Sevin, E. "On the Elastic Bending of Columns Due to Dynamic Axial Forces Including Effects of Axial Inertia," J. Appl. Mech., Vol. 27, No. 1 (March, 1960) 125-131. Srinivasan, A. V. "Dynamic Stability of Beam Columns," J.A.I.A.A., Vol. 5, No. 9 (Sept., 1967) 1685. Stuiver, W. "On the Buckling of Rings Subject to Impulsive Pressures," J. Appl. Mech., Vol. 32, Ziegler, H. "Die Stabilitatskriterium der Elastome- chanik," Ingenieur-Archiv, Vol. 20 (1952). Ziegler, H. "Linear Elastic Stability," ZAMP, Vol. 4 (1953). Unpublished Papers Huddleston, J. V. "Effect of Axial Strain on Buckling and Post Buckling Behavior of Elastic Columns," paper presented at the Southwestern Conference on Theoretical and Applied Mechanics, Tulane University, Feb., 1968. 3 1293 031691714 In“ L" II‘II “ IIII'] III-Ill Al” ". u I." "