THE UNIVERSITY OF MICHIGAN Technical Report 18 A NETWORK-VARIATIONAL BASIS FOR GENERALIZED COMPUTER REPRESENTATION OF MULTIFREEDOM, CONSTRAINED, MECHANICAL SYSTEMS Milton A. Chace CONCOMP: Research in Conversational Use of Computers ORA Project 07449 F.H. Westervelt, Director supported by: DEPARTMENT OF DEFENSE ADVANCED RESEARCH PROJECTS AGENCY WASHINGTON, D.C. CONTRACT NO. DA-49-083 OSA-3050 ARPA ORDER NO. 716 administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR May 1969

,Z 0 Z- s SI -1,1

ABSTRACT A vital component of computer-aided engineering design is the base program which computes the behavior of an arbitrary design, given a minimal input of both the structural identity and the design parameters. This paper considers the computer-a: ded design of multifreedom, constrained mechanical systems (realistic machinery). Characteristics of such systems and their computational representation and graphic display output are discussed in terms of an example machine system. An outline of mathematical methods useful for generalized representation is made. Finally an organizational scheme for implementation of a generalized program is briefly discussed. A preliminary program based on these methods is presently under test. -iii

ACKNOWLEDGMENTS The author gratefully acknowledges the financial support of the Advanced Research Projects Agency, Department of Defense; and the Advanced Analytical Technology Department, Ford Motor Company. The graphic display implementation and much of the additional programming was performed by Mr. Michael Korybalski, a senior in the Department of Mechanical Engineering, University of Michigan. -V

TABLE OF CONTENTS iii ABSTRACT................. 1. INTRODUCTION................. 1 2. COMPUTER REPRESENTATION OF AN EXAMPLE DYNAMIC MECHANICAL SYSTEM. 5 3. REPRESENTATION OF CONSTRAINT BY A PARTCONTACT NETWORK..*...... 17 4. REPRESENTATION OF FORCE INTERACTION BY LAGRANGE'S EQUATION....... 27 5. REMARKS ON PROGRAM ORGANIZATION. 35 REFERENCES................... 37 -vii

LIST OF FIGURES Figure 1. An Example Multifreedom, Constrained Mechanical System. This device is actually used as a torsional vibration experiment, but serves to illustrate characteristics of realistic machinery systems generally..6 Figure 2. Sketch of the System of Figure 1, Showing Notation Used for Design Parameters and Dependent Coordinates........-... 7 Figure 3. Display of a Schematic Representation of the System of Figure 1, on an IBM 2250 Graphic Display Unit........... 11 Figure 4. A series of Photographs of the IBM 2250 Graphic Display Screen, During Simulated Start up of the System of Figure 1. Line segments fixed in the pulleys and flywheels show the angular displacements of the respective parts. The input pulley has been repositioned (lower left corner) and the belt is not shown. Simulated values of motor strength, belt stiffness and shaft stiffness are somewhat lower than actual to emphasize multifreedom effects... 13 Figure 5. General Situation of a Rigid Part, Showing Possible Rotation, Translation, and RotationTranslation Higher Pair Contact. 18 Figure 6. Schematic Representation of a Keyed Assembly, Defined as Any Group of Parts Related Entirely by Translation Contacts.. 21 Figure 7. Abstract Example Mechanical Network. Small circles are rotation contacts; small triangles -ix

are Translation Contacts. The straight lines and large polygonal shapes denote parts. Topologically, contacts arebranches; parts are nodes. ~. 21 Figure 8. Loop Structure Defined by Network Analysis for the Abstract Example of Figure 7. ~ 23 Figure 9. Abstract Representation of a Characteristic Situation in the Computation of Kinetic Energy Terms in Lagrange's Equation. The chain of parts and keyed assemblies shown runs from ground to the center of mass of part m. The chain is only a portion of a larger total mechanical system. Terms r. 29 and s.. are noted............ 13 Figure 10. Same Illustration as Figure 9, Notated to Define Terms gQ, Sap and qk..... 31 - kx

LIST OF TABLES Table 1. Distribution of Employment of Professional Engineers in the United States by Product or Service. 2 -xi

1. INTRODUCTION Effective computer-aided engineering design requires base programs which compute the behavior or arbitrary designs, given a minimal input of both the structural identity and the design parameters which define the design. This type of program is familiar in electrical engineering. Base programs such as -ECAP, SCEPTRE, CIRAN, various proprietary programs, and others have come into conventioanl use for electrical circuit design. Aeronautical, civil, and mechanical engineers have begun use of programs such as STRESS and SAMIS for structural analysis. These programs save months of design time for designs which cannot be established entirely from experience and intuitiotn, because they allow immediate construction of detailed mathematical prototypes. Since the programs accommodate automatically to changes in structural identity, a designer can, for example, change the elements and interconnections of elements in an electric circuit without having to recode the program representing the circuit. Structural or topological change is thus accomplished as easily as change in the "size" parameters for a specific design. Once such base programs are available the way is open for productive use of various applicationindependent CAD facilities such as'graphic display, remote processing, optimization subroutines, etc. A significant fraction of engineering design involves -1

-2CONSTRUCTION' ELECTRICAL and ELECTRONICS PRODUCTS SERVICES(TransportationtEducotion Communicotions, Utilitie s,Etc. ) AIRCRAFT AND SPACECRAFT 12 CHEMICAL PRODUCTS 10 MACHINERY 6 MIFING 1 TRANSPORTATION EQUIPMENT 4 FABRICATED METAL PRODUCTS 3 PRIMARY METALS 3 ALL OTHER 10,, _ II I I I.L, I 1.I I I IL PERCENT 0 10 aD Table 1. Distribution of Employment of Professional Engineers in the United States by Product or Service. realistic machinery. Examples of such systems include vehicle suspension systems, high-speed printers, canning machinery, and machine tools. Referring to Table 1, approximately 6% of all professional engineers in 1967 regarded their principal employment product or service as machinery.* Design of machines or machine-like systems is also involved in other categories as transportation equipment, aircraft and spacecraft, and chemical products. However, there are few generalized base programs available for computer-aided design of realistic machinery. Those that are available are critically *Table 1 is obtained from the results of a survey of 53,000 professionalengineers performed in 1967 by the Engineers Joint Council under contract to the National Science Foundtion. These and other results are summarized in Engineer, March-April, 1969.

-3limited in breadth of application and convenience of use. In one class are programs including DYANA which represent linear, unconstrained multifreedom dynamic systems, but these either cease to apply or require equation formulation of the user for nonlinear and constrained systems. Since machinery characteristically involves constraint, this class of programs is limited in application. In a second class are programs including KAM which represent zerofreedom systems (ideal machinery) involving detailed twoor three-dimensional constraint. But most realistic machinery exhibits multifreedom behavior in some respect. In fact, the possibility of such behavior may be the basic cause for concern in the design. Thus what is required is an approach and an associated computer program which will represent nonlinear, constrained, multifreedom dynamic systems. Ideally, the program should be fast enough for use as a base for computer-controlled optimization. In this paper an example CAD representation of a characteristic realistic machine is outlined —including system definition, network representation of constraint and dynamic conditions, integration and output by graphic display. Several distinct features of mechanical dynamic systems are then outlined, as they compare with features of electrical and structural systems. Mechanical dynamic systems in general can be represented mathematically by a part-contact network, defining rigid body constraint, and

-4by variationally derived Lagrange equations with constraint, defining force interrelations. This representation is discussed briefly in Sections 3 and 4. In Section 5 an organization for program implementation in FORTRAN IV is briefly discussed. A preliminary generalized two-dimensional program based on these methods is presently under test.

-52. COMPUTER REPRESENTATION OF AN EXAMPLE DYNAMIC MECHANICAL SYSTEM Figure 1 is a photograph of a mechanical system used as part of a torsional vibration experiment.in the Mechanical Analysis Laboratory, Department of Mechanical Engineering, University of Michigan. The system is instrumented for experimental measurement of torsional vibration of a large flywheel on a long shaft, for various input motor speeds. Although this is a somewhat less detailed system than most realistic machinery, it illustrates several important features and provides a checkout example for a generalized program. Figure 2 shows the notation used in representing this system mathematically. It is assumed that elements 3 and 7 are significantly elastic and that the motor inputs a torque which is significantly dependent on the angular displacement and velocity of element 2. This is therefore a three-degree-of-freedom system, requiring specification of at least three coordinates for complete definition of geometry. The coordinates might be 02, 04 and 08. However, the four-bar constraint arrangement represented by elements 4, 5, and 6 appears to require. expression of the dynamic characteristics of elements 5 and 6 relative to element 4. Here this is feasible, but if the constraint were more detailed it would be difficult. Also, a technique

Figure 1. An Example Multifreedom, Constrained Mechanical System. This device is actually used as a torsional vibration experiment, but serves to illustrate characteristics of realistic machinery systems generally. -6

rm 5 r6~~ r6~~~~~~~~~~~r Figure 2. Notach Of the SYStem Of Figure ISown tiorin Used for Design Paramtr n Dependent Coordinates. eer n

-8is sought that can be implemented without major difficulty in a generalized program. Therefore a Lagrange multiplier approach is chosen. The system is defined by all five coordinates 02, 0 04, 06' and 68' and is represented by five Lagrange equations plus two constraint equations: d aT aT av 2/ -- + -+. (Qi dt [0ij i aei j=lJ 3aei i 2, 4, 5, 6, 8 (2.1) ij = o j = 1, 2 (2.2) Terminology is the following: T = system kinetic energy V = system potential energy Q.= non-conservative force associated with coordinate e. 1 X.= the Lagrange multiplier associated with the constraint ~j = 0 (or ~j = 0) j = 1, 2 The constraint functions 4j simply represent the horizontal and vertical components of loop closure: 1 = (rl + + + rr i =) (2.3) (rl + r4 + r5 + r6) j = O (2.4) When fully developed, equations (2.1) and (2.2) are: Lagrange equations: 2O2 + 2k3(r282 - P404)r2 = 21 (2.5)

-92"2 [I4+m5r4 ]4+(msr4 S)[C~S(84e5)~ s+sin(e4 0) s ] 2k3(r2e2 - P484)r4 + msgrcos 84 - (r4sine4)X1 + (r4 cos 84)X2 = 0 (2.6) [Is+m5s5 ]05+l (mr4S 5) cos405)-e)4-s in(e4- es6 2] +(mgsg5)cos85 - (r sin5)ll+ (rscose5)X2 0 (2.7) 2 " [16 + m6S6 ]06 - (m6 gs6)c6s6 k7(08 - 06) - (r6sin6)X1 + (r6Cos06)X2 = 0 (2.8) I8 8 + k7(8 86) = 0 (2.9) *... Constraint equations: r4sine4 04 - rssinese5- r6sinn6 6 2'- 2 r4 4 cose4 - r505 cos O r6CO66 2 0 (2.10) r4cose44 + rsCOSe5 O5 + r6Cose6 06 2 2 e2 r 4 eS2sine 0 (2.11) r404 sine 4 - r585 sin - r6 6 6 Most of the notation in th:se equations is defined in Figure 2. The subscripts denote the element with which the corresponding symbol is associated. The term T21 is the torque exerted by the motor on element 2. Equations (2.2) (or equations (2.10) and (2.11)) have been stated in second order form so that the orders of the Lagrange equations and the constraint equations are equal.

-10This allows application of a standard numerical integration process to determine the functions 0i(t) i = 2, 4, 5, 6, 8 and Xj(t) j = 1, 2. Characteristically, equations (2.5) through (2,11) are coupled in the coefficients of the high order terms. Because the coefficients are dependent on the system coordinates, it is impossible to determine a mathematical uncoupling valid over the entire range of displacemento Therefore a simultaneous solution of equations (2.5) through (2.11) for the high order quantities 8i i = 2, 4, 5, 6, 8 and X. j = 1, 2 must be performed at every integrating step, or at least sufficiently frequently to allow accurate extrapolation. A program for the integration of equations (2.5) through (2.11) was written in FORTRAN IV, using several subroutines from the IBM System/360 Scientific Subroutine Package, including HPCG, an integration subroutine based on Hamming's modified predictor-corrector method. Program output was expanded to include the coordinates of' all points, lines, and circles of interest in the system versus time. The overall program was run on the IBM 360/67 computer under MTS (Michigan Terminal System) and the results stored on disc storage for subsequent step-by-step display on an IBM 2250 graphic display unit. Figure 3 is a photograph of the 2250 unit with the computer-generated schematic on the display tube.

Figure 3. Display of a Schematic Representation of the System of Figure 1, on an IBM 2250 Graphic Display Unit.

-12The series of photographs in Figure 4 shows successive configurations over equal time intervals, from the moment the motor begins to supply an input torque to the system. The circles can be recognized as the pulley of element 2, pulley of element 4, and flywheel of element 8. Straight lines were fixed on each of the circles to make possible recognition of their orientation. The belt was not shown, but the input, coupler, and output links of the four-bar linkage are shown. The motor torque, belt stiffness, and shaft torsional stiffness corresponding to the system in Figure 1 were too high for any multifreedom effects to easily be seen on the graphic display unit. Therefore, these parameters were substantially reduced for the photographs shown in Figure 4. Flywheel 4 is seen to lag the input from pulley 2, and flywheel 8 lags both 4 and 2. The time required to integrate the system through approximately twenty rotations of pulley 2 was approximately one-half minute. The fact that the time response of a specific example of a constrained, multifreedom dynamic mechanical system can be efficiently determined and displayed is of marginal interest. However the example does suggest the nature of realistic machinery and inspires observations on situations that will be encountered in any generalized program: 1. A multifreedom representation is reasonable for a variety of reasons, including motion input of limited

Figure 4. A Series of Photographs of the IBM 2250 Graphic Display Screen, During Simulated Start Up of the System of Figure 1. Line segments fixed in the pulleys and flywheels show the angular displacements of the respective parts. The input pulley has been repositioned (lower left corner) and the belt is not shown. Simulated values of motor strength, belt stiffness and shaft stiffness are somewhat lower than actual to emphasize multifreedom effects. -13

-.14power and presence of elastic eleme'nts. Other situations occur involving passive freedoms and.artificial freedom introduced for the purpose of computing reaction force. 2. Large, non-linear displacements occur —physically unlike any phenomena in typical electric circuits or mechanical structures. 3. The system must at least be represented as twodimensional, and it is easy to imagine a need for three-dimensional representation. If an electricalmechanical analogy is made between current and force, voltage and velocity, then force and velocity mlust each be six-element vectors-: Force has three components of translational force and three of torque; velocity has three components of translational velocity and three components of angular velocity. In contrast, current and voltage are scalars. 4. There is a constraint situation suggestive of electrical or structural networks. Conceivably every element (rigid part, elastic force, inertial effect) can be represented as a branch; every contact or point of action between branches- as a node. Representation of the present example in such a way leads to a network of 25 branches and 12 independent

-15loops in which the state variables are three-element vectors (two translational components; one rotational component). If there is a general way to avoid this complexity, it may lie in emphasizing only that constraint caused by the interconnection of rigid parts. All other force fields are included in a state function (the Lagrangian) or as generalized, non-conservative forces. This approach will lead to a larger set of differential equations than is absolutely necessary. (In the example, seven differential equations were obtained, although only three would have been required if the problem was defined entirely in terms of coordinates 2, 04, and 08 ) But there is certainly no guarantee as a practical matter that a state variable approach will lead to a minimum set of equations either, considering the eliminations involved. Also, there may be fewer numerical problems due to ill-conditioned equations, the more directly the equations are related to a manifest physical situation. 5. The only mechanical elements analogous to elecde trical capacitors (i = C dte ) are inertial fields dv (f = mdt). However, of any two parts related by an inertial field, one is always ground. Because of this simplification, a Lagrange equation repre

-16sentation of inertia via the kinetic energy term seems natural. 6. For most purposes, the sources of applied force can be considered without mass. All mass can be regarded as lumped in the interconnected rigid parts. Thus in the example the timing belt is considered a massless torsional field, exerted on element 4 by element 2. Situations will arise involving elastic elements of substantial mass, but it may be possible to represent them in lumped parameter form, again separating the mass and field effects, but introducing more system coordinates. 7. Integration times are small enough for computational representation to be feasible. With the large output displacements involved, graphic display output is almost indispensable for rapid interpretation of overall response. However, small displacement response of individual coordinates (e.g., vibration amplitudes) must be output on scaled plots to be recognizedo

3. REPRESENTATION OF CONSTRAINT BY A PART-CONTACT NETWORK General dynamic mechanical systems are now to be represented by a combination of constraint conditions based on part interconnections and Lagrange equations interrelating the effects of inertia and applied forces. The discussion here will be in terms of two-dimensional systems, although most of the methods involved have direct extension to general three-dimensional application. As a first step toward developing the constraint conditions, parts and contacts will be defined. Figure 5 shows a part in the most general situation. By definition the part is rigid and inextensible. Part orientation is defined by the angle of unit vector Xi relative to i fixed in ground. The sets of points Pim and Qin are fixed relative to part i: Pim denotes locations where part i contacts other parts such as j, k, k shown; Qin denotes locations where external forces or torques are applied on part i, including inertial, elastic, frictional, electrical, and other force effects. Except in one situation, points P. and Q. are fixed relative to part i throughout the entire range of time involved in the response. The exception is so-called "higher pair" contact, in which case one part either rolls, or rolls and slips, on another part. This is illustrated by the contact of parts i and Q. In -17

-18PatA Part k Part k..-..~i Sik Part ji A A Possible Rotation, Translation, and Rotation-Translation Higher Pair Contact.

-19this case the contact points such as Pi3 are fixed relative to their part only during the instant in which they locate the actual contact. During numerical integration for overall time response the positions of these higher pair contact points relative to their parts must be updated at each integration step so that they will always very nearly locate the actual contact. Certain unit vectors must also be defined relative to the part. In three dimensions unit vectors mij and mig would be required to define directions of relative rotation between the parts, but in two dimensions all rotations may be assigned the direction k. In either case unit vectors Sik and siQ are required to define the directions of relative translation of any parts (2 and k) involved in translational contact. Such unit vectors of slip are fixed in their respective parts over the entire time range, except again in the case of higher pair contact illustrated by parts i and Z. In this instance the orientation of the slip vector must be updated at each integrating step to maintain an orientation very nearly that of the tangent vector to the contact point. Two basic contact types are defined: rotation and translation. To simplify system topology the convention is introduced that each contact relates two and only two partso Thus, a rotation contact is either a pin joint connecting two parts such as the contact at point Pil;

-20or a rolling contact (with or without slip) such as the contact at point Pi3' A translation contact is either a keyway joint connecting two parts such as the contact at point Pi2; or a slip between two surfaces such as the contact at point Pi3. If it is necessary to represent at a single joint a combined rotation-translation or a double translation, this is accomplished by interposing an arbitrarily small part between the two contacts. Thus each contact relates two parts as required. This scheme extends readily to three-dimensional analysis in which an ideal ball joint is, for example, represented by three rotation contacts separated by two zero-size parts. Graphically, rotation contacts are represented by small circles; translation contacts by small triangles. Parts are represented by straight lines or large polygonal shapes. Now another less familiar entity important in constraint statements will be defined: the keyed assembly. A keyed assembly is simply any set of one or more parts related entirely by translation contacts. A schematic keyed assembly is shown in Figure 6. The basic reason for the usefulness of the concept of a keyed assembly is that all parts in the assembly undergo common angular displacement. A part-contact network is the topological network of nodes and branches in which every part is included as a node and every contact is included as a branch. An appropriate tree derived from the part-contact

-21Part i Port I Part k Part j Figure 6. Schematic Representation of a Keyed Assembly, Defined as Any Group of Parts Related Entirely by Translation Contacts. 9 $ Figure 7. Abstract Example Mechanical Network. Small circles are rotation contacts; small triangles are translation contacts. The straight lines and large polygonal shapes denote parts. Topologically, contacts are branches; parts are nodes.

-22network forms the basis on which the dimensional constraint equations such as (2.10) and (2.11) are constructed. This is completely analogous to the process by which the Kirchoff's Voltage Law equations are constructed from electrical networks. In both instances the only topological problem is determination of an appropriate tree from the input network. This is conveniently handled by network analysis.2 The results as concern mechanical systems can be illustrated by the hypothetical mechanical network shown in Figure 7. The branch-node incidence matrix for this example is the following: Branches 1 6 3 4 5 8 2 7 9 10 12 11 1 1 1 0 0 0 0 0 0 0 0 0 0 2 -1 0 0 0 0 1 0 0 0 0 0 0 3 0 0 0 0 0 -1 0 0 1 0 0 0 4 0 0 1 0 0 0 1 0 -1 0 0 0 ~ 5 0 0 0 0 0 0 -1 0 0 1 0 0 o 6 0 0 0 1 0 0 0 1 0 -1 0 0 = [A] 7 0 -1 -1 -1.1 0 0 0 0 0 0 0 8 0 0 0 0 -1 0 0 0 0 0 1 0 9 0 0 0 0 0 0 0 0 0 0 -1 1 10 0 0 0 0 0 0 0 -1 0 0 0 -1 (3.1) Here, precedence in the sequence of columns has been given those branches (contacts) relating parts nearest ground. To obtain the loop matrix the non-augmented branch-node incidence matrix is formed by eliminating the row corresponding to the ground node (part). Then by a process similar to Gaussian elimination a cut-set matrix is developed.2

-23The transpose of this matrix is the loop matrix, listed below for the example of Figure 7. Branches 1 6 3 4 5 8 2 7 12 9 10 11 1 1 -1 1 0 0 1 0 0 0 1 0 0 a*2 0 O -1 1 0 0 1 0 0 0 1 0 = [B] O3 O O O 1 1 0 0-1 1 0 0 1 (3.2) Note that development of a tree occurs implicitly in the construction of equation (3.2). In this case a tree results in which most tree branches directly contact the ground node, or contact only one or two nodes away from ground. This tends to define loops involving a smaller number of nodes and branches than otherwise. Figure 8 outlines the loops defined by equation (3.2) Loop Structure Defined by Network Analysis for the Abstract Example of Figure 7. Loop 3 Loop I for the Abstract Example of Figure 7.

-24For any given loop there are two useful constraint conditions for mechanical systems. One states that the sum of the relative translational displacements around a closed loop is zero; the other that the continued product of the relative angular transformations around the loop is unity: Translational constraint: m z ql = 0 (3.3) i=l Angular constraint: 2 [Fln~ll [B. ili = [1] (3.4) [IFn] n [Ei i- 1 ] i=n In equation (3.3) the qi vectors are defined as being directed from one contact point to the next around the closed loop. (Note that for each contact there are two contact points, one in each of the related parts.) An individual qi vector either spans a part, from one contact to another, or spans a contact, from one part to another. The qi vector is zero if it spans a rotation contact, but is finite if it spans a translation contact. In any case it can be expressed in polar coordinates as the product of magnitude and direction: qi= qiqi = qi[co~si i + sineij] (3.5) Consideration of equations (3.3) and (3.5) suggests that

-25the dependent variables involved in the translational constraint will be the angles of parts (more specifically their keyed assemblies) relative to ground ei, and the distances of slip at the translational contacts sk92 In equation (3.4) the [EiBi 1] terms are 3 x 3 Euler transformation matrices involving as dependent variables the angle of each keyed assembly relative to the keyed assembly preceding it in the loop. The condition (3.4) is required for three-dimensional problems* but degenerates to an unnecessary condition in two dimensions: specifically that the sum of the relative angles around a closed twodimensional loop is zero. Thus all loop constraints required for two-dimensional analysis will have the following form: m m i i q. qi Cosa. = 0 (3.6) ~i=l 1i=l m m j E qi ~ 0 = i sine = 0 (3.7) j+l i=l i=l If it is intended to include higher pair contact within the application of a general program for mechanical *Even in three-dimensional problems, direct solution of equation (3.4) can be avoided by working with its first time derivative —the statement that the sum of the relative angular velocities around a closed loop of keyed assemblies is zero: Wln +. Wi,i-1 = 0 i=n

-26systems, then the part-contact network constraints are not sufficient. Constraints defining the profiles of the contacting parts are required, as well as a condition that they do not interpenetrate.

4. REPRESENTATION OF FORCE INTERACTION BY LAGRANGE'S EQUATION The constrained form of Lagrange's equation for any mechanical system can be written as follows: d aT 9V -t X Pi- -i-+ Z + Qi (4.1) dt 9n ar n = 1 This equation compares closely with equation (2.1). Here the coordinates.i represent the angles of keyed assemblies relative to ground and the slips of the translation contacts, as discussed in Section 3. As before, T and V are the total system kinetic and potential energies respectively; Qi is the generalized force associated with the ni coordinate. There are m constraint equations fj and m Lagrange multiplier functions X (t). In the construction of a generalized program the essential task is to develop and implement a procedure for reducing the Lagrange equations represented in equation (4.1) to ordinary differential equations in standard form and to efficiently compute the numerical values of the coefficients of these equations at each integrating step. The form of equation (4.1) suggests that the greatest difficulty will be encountered in computing the coefficients of the kinetic energy terms. First the total system kinetic energy due to rotation and translation must be expressed -27

-28in terms of the system coordinates (keyed assembly angles relative to ground and slips at the translation contacts). Total kinetic energy T is the sum of the kinetic energies T due to the motion of each of the individual parts. m Figure 9 depicts the geometry for computation of the kinetic energy of part m: T 1 I 2 1 2 (4.2) m 2 m 2m m m The velocity Pm must be expressed in terms of the zeroth and first derivatives of the system coordinates intervening between ground and the center of mass of part m: t. a 1 Pm ([r. + sij ] (4.3) i=l 1 j=l 1J t. t. a 1 a 1 P E 6i[k x (r. + ij) l ij ij (4.4).) i=l j i=l j=1 13 (4.4) In equations (4.3) and (4.4), the i summation runs over all the keyed assemblies in the chain, from ground to part m. The limit a is the number of keyed assemblies in the chain. The j summation runs over all the translation contacts within the ith keyed assembly, where ti is the number of translation contacts in the ith assembly. The terms r. and s.. can be 1 13 identified and related to the term 4k defined in Section 3 as follows. Vector s.i is the slip displacement of the jth translation contact counted within the ith keyed assembly. It is therefore the same as a vector qk spanning a translation

-29part m r4 keyed assembly I Uji u12 r2 s 21 Figure 9. Abstract-Representation of a Characteristic Situation in the Computation of Kinetic Energy Terms in Lagrange's Equation. The chain of parts and keyed assemblies shown runs from ground to the center of mass of part m. The chain is only a portion of a larger total mechanical system. Terms r. and tij are noted. LJ

-30contact. Vector ri is the vector sum of all the vectors spanning parts, in a given chain of vectors passing through the ith keyed assembly. It is therefore the same as the summation of all the qk vectors which span parts in the ith keyed assembly. Development of the kinetic energy terms in equation (4J1) requires squaring and performing differentiation operations on equation (4.4), in terms of the specific coordinate nR corresponding to the specific Lagrange equation under construction. When this is worked out a surprisingly simple relation is obtained: For angular coordinates, rQ = 8: d.Tm.Ta 1.. Tm _b dt XL k iTm 6,m m m m j iOl (4.5) i (b. x k) + si. (s.. x k) + 2 s.. ei(s ) ]~b For slip coordinates, qr = sQ: I-I t. dd~ W(a ti T mm 1 i(kxi- 1 1 L 2p S p i=l j=1 * ^A 0. ^ ^ ^A + s.i sij + 2 si.. i (k x Sij)AOS J (4.6) Some of the terminology of equations (4.5) and (4.6) is explained in terms of Figure 10. Following a chain of

-31part m q6 keyed assembly I / 4 Sip q2 Figure 10. Same Illustration as Figure 9, Notated to Define Terms 0,Q s 9 and 4k. bQ PQM

-32parts from ground to the part containing mass m, b is the vector spanning the keyed assembly having angle e,; s p is the unit vector of translation for the pth translation contact within keyed assembly Q. An individual qk vector is one vector segment in the chain of part and slip vectors running from ground to mass m. Equations (4.5) and (4.6) have a striking similarity to the classical four-term expression for the second derivative of a single position vector with respect to time, containing tangial, centripetal, radial and coreolis terms. In this sense they reduce the Lagrange equation of equation (4.1) to the somewhat more tractable form: ma = F (4.7) The potential force and constraint terms introduce little programming difficulty provided a limited variety of such effects are included. For example, the potential term might represent forces due to gravity and linear torsional and translational springs. The constraint term would account only for part-contact network constraint, unless higher pair contact is involved. The Qi term poses a more detailed computational problem. Assume that a total of p actual non-conservative forces Fk act on the system. The generalized force Qi is computed as a force acting directly o-n the coordinate ni, equivalent in effect to that of all the actual forces Fk:

-33pi k pi (4.8) k=l 1 Usually ik can be conceived as the vector between the two points of attachment of k.' Given a series of tree branches connecting these two points, the partial derivative terms are readily evaluated. If the Fk terms include only forces dependent on displacement and velocity, they are computed from the motion data obtained in the prior integration step. Dependence of Fk terms on reaction force (Coulomb friction) does not fit this pattern and may require introduction of a new constraint, for which the associated Lagrange multiplier corresponds to the required reaction force. As in equation (2.2) of the example problem, the constraint equations accompanying the Lagrange equations are stated in second order form so that the orders of all equations will be the same. The second derivative of equation (3.3) is: m E q. = 0 (4.9) i=l The expanded expression for ~i is: 2 - qi =[qi 8- 6i qi]qi +qi 2 qi i](k x qi) (4.10) The second order constraint conditions are formed as the orthogonal components of equation (4.9). If the partcontact loop includes a higher pair constraint, an additional accelerational term is required in terms of the

-34radii of curvature of the contacting surfaces. The above analysis method extends directly to three dimensional problems. Use would then be made of expressions similar in form to equations (4.5), (4.6) and (4.10).

S. REMARKS ON PROGRAM ORGANIZATION A preliminary generalized two-dimensional program based on the foregoing methods has been written and is presently under test. The program is written entirely in Fortran IV for compatibility, and therefore has a somewhat primitive data structure compared with what might have been utilized in PL/1 for example. Basically, all fixed and floating point data are stored in tables identified under one of the following categories, termed entities: 1. Fields. Conservative or non-conservative applied forces exerted on one part by another. 2. Generatorso Motion sources, inputing either a rotation to a keyed assembly or a slip to a translation. 3. Keyed Assemblies. Assemblies of one or more parts related entirely by translation contacts. (Internal use only.) 4, Loops. Sequences of parts and contacts forming an independent closed chain. (Internal use only.) 5. Markers. Points, unit vectors or points and unit vectors, fixed in translational or orientational position relative to a specified part. 6. Parts. Rigid, inextensible objects of zero or finite size, mass and moment of inertia. -35

-367. Rotations. Contacts allowing single-freedom rotation between parts. 8. Requests. Indications by the user of the pro_ gram of desired output (displacements, motions, forces versus time). 9. System. A category established under which information pertinent only to the overall system is stored. 10. Translations. Contacts allowing single-freedom straight line slip between parts. For each entity there are two tables, one for fixed point data; the other for floating point data. Variable length lists are referenced by pointers in the tables and stored end-to-end in long, one-dimensional data arrays. A set of simple access subroutines is employed to store and retrieve data in the tables.

REFERENCES 1. H20-0205-1, System/360 Scientific Subroutine Package (360A-CM-03XO), Version II, Programmer's Manual, IBM Corporation, Technical Publications Department, White Plains, N.Y. 2. Calahan, D. A., Computer-aided Network Design, McGrawHill, New York, 1968. -37

-38UNCLASSIFIED Securit Cfa s sification DOCUMENT CONTROL DATA - R & D (Security classification of title, body of abstract and indexing annotation must be entered when the overall report is classified 1 ORIGINATING ACTIVITY (Corporate author) tZ. REPORT SECURITY CLASSIFICATION THE UNIVERSITY OF MICHIGAN Unclassified CONCOMP PROJECT 2b. GROUP 3. REPORT TITLE A NETWORK-VARIATIONAL BASIS FOR GENERALIZED COMPUTER REPRESENTATION OF MULTIFREEDOM, CONSTRAINED, MECHANICAL SYSTEMS 4. DESCRIPTIVE NOTES (Type of report and inclusive dates) Technical Report 18 5. AU rOR(S) (First name, middle initial, last name) MILTON A. CHACE 6f REPORT DATE 78. TOTAL NO. OF PAGES 7b. NO. OF REFS May 1969 37 2 8a. CONTRACT OR GRANT NO. 9a. ORIGINATOR'S REPORT NUMBER(S) DA-49-083 OSA-3050 b. PROJECT NO. Technical Report 18 9b. OTHER REPORT NO(S) (Any other numbers that rrmuv bh assrpr'eit this report) DISTRIBUTION S T ATEMENT Qualified requesters may obtain copies of this report from DDC:1. -FdPPLEMENTARY NOTES 12. SPONSORING MILITARY ACT VITY Advanced Research Projects Agency 1 A B 5 i P.C T A vital component of computer-aided engineering design is the base program which computes the behavior of an arbitrary design, given a minimal input of both the structural identity and the design parameters. This paper considers the computer-aided design of multifreedom, constrained mechanical systems (realistic machinery). Characteristics of such systems and their computational representation and graphic display output are discussed in terms of an example machine system. An outline of mathematical methods useful for generalized representation is made. Finally an organizational scheme for implementation of a generalized program is briefly discussed. A preliminary program based on these methods is presently under test. DD, NV!NO 7 Unclassified Secitx AEa 1slI_'c _ _inn

-39Security Classification 14 r LINK A, LINK B L INI C KEY WORDS ROLE WT ROLE WT ROLE WT Sectit ClaI l A ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~I Security Classification

UNIVERSITY OF MICHIGAN 3 901 5 02654 031311111 111115 02654 0313