SIMPLIFICATIONS IN THE SIMULATION OF MECHANISMS CONTAINING FLEXIBLE MEMBERS by Thomas John Wielenga A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Mechanical Engineering) in The University of Michigan 1984-. Doctoral Committee: Professor Milton Chace, Co-Chairman Assistant Professor Joseph Whitesell, Co-Chairman Professor William Anderson Associate Professor Noboru Kikuchi Assistant Professor Mohammed Zarrugh

Thomas John Wielenga 1 All Rights Reserved

ACKNOWLEDGEMENTS I would like to express my gratitude to those who helped and encouraged me in this endeavor. First of all to my committee for their role in examining this dissertation. I would especially like to thank my co-chairmen, Milton A. Chace, and Joseph Whitesell. Professor Milton Chace consistently obtained financial support and research material through the years which were essential to my progress. I have also found him to be a loyal friend. Assistant Professor Joseph Whitesell provided encouragement and hope at a particularly difficult point in time. He also provided major help in discussions about the size of various terms in the equations. I would also like to thank Professor Jens Wittenberg for a particularly helpful discussion at a formative stage. His notation and perspective has been valuable in keeping the equations simple. I also appreciate the many discussions with Rajiiv Rampalli about the host of details involved in simulation which are important, but not at all obvious. I would also like to thank the corporations which supported me through various research grants. Deere and Company initiated my research through their support and guided it to some extent. The Robotics Division of CRIM at The University of Michigan provided help in supporting Brian Dent to code a two dimensional version of the equations. His help was greatly appreciated. CRIM also provided the means to type ii u

and publish this dissertation through the able help of Mary Ann Pruder. I would also like to thank Mechanical Dynamics Incorporated for allowing me to work on this dissertation while employed there. This support was pivotal in the final stages. Finally I would like to thank my wife Deborah and my parents for their unfailing support, encouragement, and steadfast faith in me. 111 11!

ABSTRACT SIMPLIFICATIONS IN THE SIMULATION OF MECHANISMS CONTAINING FLEXIBLE MEMBERS by Thomas John Wielenga Co-Chairmen: Milton Chace, Joseph Whitesell The equations of motion for linearly elastic bodies undergoing large displacement motion are derived. The equations of various types of joints and forces are also derived. This produces a set of equations which are efficient to numerically integrate and provides a general approach to mechanism simulation. The equations for the elastic members are formulated and simplified to provide as much efficiency as possible in their numerical solution. The techniques used are: 1) A local frame of reference is associated with each body. In this frame deformation retains its linearity, despite the non-linear overall motion of the body. 2) The equations of motion for the rotation of the frame are formulated in terms of angular velocity instead of angular coordinates. This allows much simpler vector equations. 3) By making the local reference frame "float" with respect to the body many of the terms in the equations are eliminated. This can be done by choosing free modes for deformation coordinates, or by constructing a set of coordinates with similar properties by using Gram-Schmidt orthogonalization.

4) Many of the remaining terms are shown to be of negligible size and are eliminated from the equations. A further efficiency is obtained through the use of Euler Parameters as angular coordinates. These coordinates have no singular positions and are easy to evaluate and integrate. The equations are presented in two forms for numerical integration. The first form is suitable for integrators requiring the equations to be explicitly solved for the first derivatives. This form is appropriate for non-stiff systems of equations. The second form is appropriate for those integration techniques which allow implicit equations and is suited to stiff sets of equations. The implicit numerical integration is extended to solve second order equations, further reducing the numerical effort required. Several examples are solved and the results are given. The formulation given is seen to be accurate and is expected to be efficient for many types of problems.

TABLE OF CONTENTS ACKNOWLEDGEMENTS...........................................................ii LIST OF FIGURES..*...*........... *............g........... 6 t.e............. Vii LIST OF TABLES...*******.........*...***************.. **......... ********ix LIST OF APPENDICES............................................................ x CHAPTER 1. PIrobleUC Definiio............................................................ 1 1.1 Problem Definition...................... **1 1.2 Previous Work.................................................................. 2 1.3 Dissertation Overview................................................ 5 1.4 Notational Preliminaries.................................................. 1.5 The Vector Application of the Fundamental Equation........................................................................ 11 1.8 Moving Frames of Reference............................ 13 1.7 Rotational Coordinates.................................................. 16 1.8 Numerical Integration................................................... 20 II. THE SINGLE RIGID BODY............................................... 28 2.1 Introduction...................................... 28 2.2 The Local Reference Frame and Vector Definitions....... 2 2.3 Linear Momentum............................................... 28 2.4 Angular Momentum............................................ 30 2.5 Kinetic Energy............................................................... 32 2.6 Application of the Fundamental Equation....................... 33 2.7 Matrix Form of the Equations........................................ 37 iv

II. THE SINGLE FLEXIBLE BODY ~~~~~~~~~~0~~,~~~~~~~~~~~~~~~~~0~ 42 3.1 Introduction.*.......................**...*.....................eeeee. ***,***42 3.2 Vector Definitions and Deformation Coordinates........... 42 3o3 Linear Momentum.4................5...... 16 3.3 Linear Momentum *oe................................................. 45 3.4 Angular Momentum........................................................... 47 3.5 Kinetic Energy.................................................. 50 3.8 Strain Energy.......................................*****. 52 3.7 Application of the Fundamental Equation...................... 53 IV. SIMPLIFICATIONS FOR DIFFERENT FRAMES OF REFERENCE................................................................. 80 4.1 Introduction oo...................................................................... 0 4.2 The Locally Attached Frame.............................61 4.3 The Floating Frames.................................................... 67 4.4 The BucTssensrand Frame....................................................... 75 4.5 The Tisserand Frame. ****................ e e.............. 77 4.8 The Free Modes as Deformation Coordinates................. 79 4.7 Simplifications Using Gram-Schmidt Orthogonalization......................................... 82 V. SIMPLIFICATIONS DUE TO DROPPING NEGLIGIBLE TERMS............................................. 84 5.1 Introduction...................................................................... 84 5.2 Sizing Assumptions...................................................... 84 VI. THE MECHANISM AS A WHOLE........................................ 94 8.1 Assembling Bodies, Joints and Forces............................ 94 6.2 Constraint Equations..........8.............................. 98 6.3 Joints Foe 03.................................................103 6.3 Applied Forces........................................ 105 V

VII. NUMERICAL RESULTS AND CONCLUSIONS................. 108 7.1 7.2 7.3 Introduction............... Testing and Validation Conclusions................. 108 109 120 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ APPENDICES. REFERENCES 122 143 vi

LIST OF FIGURES FIGUR: 1.1 1.2 1.3 2.1 3.1 4.1 4.2 4.3 5.1 5.2 6.1 0.2 6.3 8.4 6.5 E Vector definitions for a moving frame Planar rotation............................. Rotation around an axis................. Vector definitions for a rigid body..... Vector definitions for a flexible body Two mass problem............................. Locally attached frame.................. Floating frame.................................. Planar three mass problem................. Mode shapes for three mass problem Spherical joint................................ Universal joint..................................... Translation constraint....................... Joint types......................................... Line force........................................ * --- —-— e-eeeeeeeeeeeeeeeeeeeeeeeeee I.................................... 14 17 19 27 43 61 62 70 89 90 97 101 102 104 106 *VIl

7.1 Free vibration..........*..*..*****........*********.********* **** ***. 110 7.2 Boundary condition test cases *******............................. 112 7.3 Center point deflection of connecting rod for slider-crank (From Reference 7)............................ 114 7.4 Center point deflection of connecting rod for slider-crank (From DYMOS).............................. 115 7.5 Angular velocity of triangle for level sero.......................... 116 7.6 Angular velocity of triangle for level one-half....................... 117 7.7 First mode of triangle for level zero.................................... 118 7.8 First mode of triangle for level one-half................................ 119 B.1 Spherical Joint.................................................................... 134 viii V111

LIST OF TABLES TABLE 5.1 Size of terms for frame rotation............................................. 87 5.2 Size of terms for modal equations...................................... 87 5.3 Body constants for three mass problem................................. 91 5.4 Level one-half equations for flexible body.............................. 93 6.1 Constraint equations for joints............................................... 105 7.1 Boundary condition test results....................................... 111 ix

LIST OF APPENDICES APPENDIX A. The Second Order Implicit Integrator.................................... 123 B. The Spherical Joint Constraint Derivation......................... 133 x

CHAPTER I INTRODUCTION 1.1. Problem Definition The simulation of dynamical mechanisms undergoing large displacements involves the derivation and solution of highly non-linear equations of motion. In the past, the equations of motion for a specific problem were derived by hand, programmed and numerically integrated. Recently however, computer programs which automatically generate and numerically integrate the equations of motion have been developed. Similar progress has been achieved for problems undergoing small displacement motion, assumed to be in the elastic range of component materials. In this case, linear dynamical models are automatically formed through the finite element method and solved with modal analysis techniques. Typically, mechanisms composed of rigid bodies are modeled with a small set of equations which are highly non-linear. Dynamical elastic structures, on the other hand, are often represented with linear models having several thousand degrees of freedom. Recently these separate fields have begun to coalesce. The assumption that all of the members in the mechanism are rigid is not always a good one. The members are actually flexible structures and sometimes the interaction between the member's flexibility and the mechanism's overall performance is important. This is especially true in I

2 high speed mechanisms and problems involving impact. On the other hand one of the biggest problems in finite element analysis is determining the loads on the structure. These loads can be calculated through mechanism simulation. If the structure is a moving part of a mechanism, its own inertia also introduces significant loads. Those considerations cause us to view the structure as an integral part of the solution process. The dynamics of the member are considered simultaneously with the dynamics of the mechanism. This can cause significant problems. Problems arise because the large linear problem of structural dynamics becomes enmeshed in the highly non-linear dynamics of overall motion. If handled improperly, this can result in huge non-linear sets of equations. It should not be surprising that techniques from both fields can be used to simplify the equations of motion. This is the subject of this dissertation: Simplifications in the simulation of mechanisms containing flexible members. 1.2. Previous Work There have been many investigators pursuing this field since the early 1970's. The early researchers tried to solve the problem by considering the dynamics of mechanism separately from the response of the individual members. The dynamics of the mechanism were considered first with the assumption that the members could be considered rigid. Once the motion and forces were determined from rigid body dynamics, these were applied to the individual members and each member's response was determined. During this phase of the analysis a finite element model usually represented the dynamics of each individual member. This type of analysis was called Kineto-Elastodynamics [1-4].

3 Since these approaches separated the overall system dynamics from the dynamics of the flexible members, the accuracy of the results was limited. The largest disparity between theory and measurements was the magnitude of reaction forces at joints in the system. Because the flexible members absorbed energy during impact, for example, the actual forces were lower than computed forces. One approach that followed formulated the problem in coupled form analytically. This was a very tedious process and was done only for simple mechanisms [36]. Another approach to formulating the coupled problem was to formulate finite elements with large displacement coordinates [7]. This provided the missing coupling, but presented a larger problem. In order to model a member in detail many elements would be required. This made the number of coordinates to integrate astronomical. An approach which eliminated many of these coordinates was presented by Dubowsky, et al. [8,9]. A method called Component Mode Synthesis was used to reduce the number of coordinates by doing a modal reduction on the internal degrees of freedom. Modal coordinates for the lower frequencies were used to represent the deformation of the body. The interface nodes were constrained to move with interface nodes of other bodies to tie the system of bodies into a mechanism. The method eliminated the constraints and constraint forces and produced a minimal set of coordinates. This produced a highly non-linear and coupled set of equations to integrate. Similar approaches were proposed using Component Mode Synthesis by Shabana/Wehage, [10] and also Benson [11]. These used different approaches in the formulation of the overall equations and solution procedures. Another set of investigators in the spacecraft industry were solving the problem of flexible spacecraft by attaching a local reference frame to an assumed rigid body in

4 the satellite [14-16]. All of the deformations were measured with respect to this rigid body. Various improvements followed usually with this fundamental assumption. The equations of motion generated were fairly complex. One set of investigators [13] produced a program (called DISCOS) capable of simulating satellites composed of coupled flexible members. This is a fairly ponderous program using NASTRAN and DMAP. Here the investigators again used a locally attached frame in each flexible body and deformation coordinates which could be either finite element mode coordinates or modal coordinates. Only hinges are allowed as joints. Also in the spacecraft industry, DeVeubeke [17] showed that the single flexible body could be represented using a floating frame of reference. He examined two types of floating frames, the Tisserand frame and the Buckens frame. He showed that the free modes satisfy the Bucken's constraint equations and dramatically simplified the equations of motion. Canavin and Likins [18] in a later work claimed that for a linearly elastic member the Tisserand frame was approximately the same as the Buckens frame. This approximation allowed them to drop a non-linear coupling term from the constraint equations. This further simplified the equations for a single elastic body. This particular simplification does not seem to be justified in all situations and will not be used in the development in this dissertation. The development here will use the Buckens reference frame for each flexible body and free modes (or modes with similar properties) to represent the deformation. This gives substantial simplification for each body's equations. In addition, the remaining terms will be checked for their relative contributions, and terms with negligible contri

5 butions will be dropped. This gives more consistent simplifications than those presented by Canavin. Techniques from rigid body dynamics for the solution of constrained mechanical systems will be used to solve the equations. Several developments are unique to this dissertation. Determining the relative size of terms in the equations and dropping those that are negligible is one contribution. Another unique contribution is the presentation of the equations in two efficient forms: one form suitable for non-stiff systems, and one form suitable for stiff systems. An additional contribution is the presentation of an integrator which will integrate a mixed set of implicit equations consisting of second' order equations, first order equations and algebraic equations. This integrator is especially appropriate if the system is stiff. 1.3. Dissertation Overview The dissertation is organized as follows: First, preliminary material to the understanding of the overall problem is given in Chapter I. This includes mathematical notation, general concepts of moving reference frames, techniques for applying Lagrange's equations in vector form, and integration techniques. The second chapter derives some of the inertial properties of rigid bodies and their equations of motion. It also shows how they can be simplified by an appropriate choice of local reference frame. The third chapter derives the properties of flexible bodies and examines the newly added terms. The equations of motion are presented without any simplifications. In the fourth chapter these equations are simplified by using deformating on coordinates that are appropriate for a floating frame.

6 The fifth chapter presents further simplifications resulting from eliminating terms which are negligible from the equations. The sixth chapter shows how the individual bodies are assembled to make up a mechanism. Various constraint equations for different types of joints are given as well as a basic forcing function. Examples are presented in chapter seven and the strengths and weaknesses of the various approaches are discussed. 1.4. Notational Preliminaries The notation used in this dissertation will be a combination of several different kinds. First of all, vectors and tensors will be used in describing the overall motion and orientation of each body. The equations of motion for a body's' overall motion will be derived in vector form. This is the standard form for traditional rigid body dynamics. To program the equations, however, matrix form is preferred, and the vector equations can be converted fairly readily to matrix form using the techniques of this section. Most of these techniques are described more fully in Reference [22]. Along with the overall motion of a flexible body, the deformation of the body must be described. For the deformation coordinates of the body, index notation will be used. This seems to be the least confusing way to handle the deformation if the number of subscripts is kept low. A vector will be denoted by an arrow above the vector quantity. It can be represented as a linear combination of three unit basis vectors: v = ve + v2e2 + 3e (.4.1)

7 A matrix will be denoted by an underbar. For example v = [vI v2 v3lr. Equation 1.4.1 can be represented in matrix notation as -+ T V= v e -- (T (1.4.2) or v=e t where e is a column matrix with unit vectors as elements: = ['l 2 e AJ The column matrix v contains the coordinates of the vector v in base c. Unit vectors are orthogonal and of unit length. This can be stated in matrix notation: e I e * e2 [e e2 e3. e 3 This matrix operation is done in the normal way, except that the product between terms is the dot product. Expanding this expression into a matrix yields the following: e e'e= I l'e2= T - -+ -O -* -4 -- -J C ' 81 e3'2 Ce3'e Since the unit vectors are orthogonal, all the off-diagonal terms are zero, and the diagonal is unit valued. The expression becomes 1 0 e e 01 0 =I (1.4.4) 001, where I is the 3 x 3 identity matrix. More than one set of base vectors will be needed. In order to distinguish between sets of base vectors, a superscript in parenthesis identifies the particular base. For example e(S) are the base unit vectors of base 8. The set of base vectors itself is referred to as a reference frame. The unit vectors of one reference frame are linear combinations of the unit vectors in another reference frame. This can be written in

8 the following way: =() A 's (S) (1.4.5) where A ' is the direction cosine matrix. The direction cosine matrix has an important property. The inverse of A ' is equal to its transpose: (A rs)-1 (As)T =-A (1.4.6) If the vector v is resolved in two different reference frames, it will have different coordinates in each frame. The vector resolved in two different reference frames is: vs e )r () ()r (r) Using Equation 1.4.5 for e'), this equation becomes: -s)r v(s) ( e (s)r A sV(r) (1.4.7) or v(s) - Asr V( The direction cosine matrix then, also transforms the coordinates of one frame into the coordinates of another frame. For this reason it is also called the transformation matrix. The product of two different sets of unit vectors will not result in the identity matrix (as in Equation 1.4.4), but results in the transformation matrix. (r). (S)r A Ar(S). -()r = A'sI (1.4.8) = A's. Dyadics will be used in addition to vectors. A dyadic is also called a "second order tensor" (this is shortened to "tensor" in most cases). A dyadic is a sum of indeterminate products of two vectors each, and is denoted by a double arrow. D = alb + a2b + 2 (1.4.9) It is also a linear vector operator. It produces a vector when a vector is dotted with it.

9 D * v _s (-)+ 2(b2 =)+-a or U * D= (v )b + f (v*)b2 + - -D The unit tensor (I) has the property that it returns the same vector back again: i * u= 7v (1.4.10) and * - v= (1.4.11) It is possible to resolve a tensor in some coordinate frame. D can be written D =-D Dele + Dl2ele2 + D13e1e3 + D2e2e1 + D22e2e2 + D23e2e3 + D lee + D 32e 2 + D 3ese Rewriting in terms of matrix notation this becomes D le De (1.4.12) where D11 D12 D13 D= D2l D22 D23 D31 D32 D33 and D is called the coordinate matrix of D in frame e. The coordinates of D can be obtained in another frame by a transformation. Using Equation 1.4.12:.(s)r s O (S) = +(,)rD' e(') and Equation 1.4.5, -(s)r Ds (s) = ()As r D As e(s). e"sf D S e(s)e(s)t A rsD " A rs,s So: D' — A D A T (1.4.13) or D' = AsD'A. (1.4.14) The double vector cross product can be written in terms of a tensor:

10 a X (i X b) = (a b) - b(a * ) = [(;* b)I- b j-i (1.4.15) = D v where D = ( * b)I - j (1.416) or D- arb I-b aT The transpose of a tensor is defined D v. = -. D In terms of the matrix of coordinates, D becomes D. A tensor is symmetric if D T = D or D r D. A symmetric tensor's coordinate matrix can be reduced to diagonal form, by choosing an appropriate frame in which to resolve the tensor. The axes of the frame are called the principal axes of the tensor. A tensor is anti-symmetric if D r -D or D = -D. It has only three independent components. An example is defined by the cross product operator. n x v =n ~ v (1.4.17) f has the following coordinate matrix: F 0 -fs nf2 = ns 0 -ni (1.4.18) — n2 l 0 Equation 1.4.17 written in matrix notion is nv (1.4.19) The double vector cross product a X (b X;) could be written a b v. Comparing this with equation 1.4.16 ab b a aTb I (1.4.20) Calculating a b from this identity takes 12 multiplies as opposed to 27.

11 The cross product (a X ) X = (- a) v written [a blv. The quantity [a b\ can be written [abj ^ aTabT -(1.4.21) [a bl b a -a b An identity used fairly regularly to simplify equations is (1.4.22) b a + [a b] a b This can be seen fairly readily by expanding both sides using Equations 1.4.20 and 1.4.21. As the body moves and translates in an overall way it can also deform (in the flexible analysis). For flexible bodies in general, an infinite number of coordinates are possible. Usually, however, approximations are made which reduce the number of coordinates to a finite number. For these deformation coordinates, indicial notation will be used. The range of the subscript will usually run over the number of deformation coordinates. Repeated subscripts in an expression indicates summation. For example _~ N.. U- mini liA i U =,., -- S, A subscript in parenthesis indicates no summation on this index. U = w(s,2i The Kronecker delta is defined 8 1 -i ] s -= (0 otherwise and it also acts to interchange indices

12 i 6j- = u, (1.4.23) 1.5. The Vector Application of the Fundamental Equation The fundamental equation for dynamics can take many forms. The form used here is amenable to getting both vector equations and scalar equations. This make the resulting equations more compact and individual terms have more recognizable meaning. The fundamental equation can be written (Reference 24, Equation 8.6.1): N d OT T OV -JX k OJ q t s= O0 K =1...NC The first line of the equation has the terms in the system relating to the kinetic energy (T), and potential energy (V). The second line contains the generalized forces of constraints which were left out of the original variations. For every constraint force XB included here, there must be a corresponding constraint 4k included with the equations of motion. These constraints are shown included at the end. In the third line are the generalized forces. These terms are for applied forces left out of the potential energy function. Fe is the applied force and P, is the corresponding point of application or spanning vector. In this equation the quantities 6qj, Xk,,I, F,, P, can be either scalars or vectors. However they must be consistent. When Xt is a vector, then 4A must also be a vector. Likewise when F, is a vector PC must also be.

13 The products involving the dots must also be consistent. If Xt and <k are vectors the dot means a vector dot product in this term. If they are both scalar valued, then the dot means normal scalar multiplication. This is also true for the product involving F,, and also for the final trailing product with the variations. If the terms are vectors, a dot product is meant; otherwise scalar multiplication is meant. When a variation 6qj is a vector, a vector equation of motion will result. This necessitates taking partials with respect to a vector. The partials of a scalar with respect to a vector variation will give a vector quantity. The simplest way to obtain the vector result, is to take the variation of the scalar, collect like variations, and put it into the form of a coefficient dotted with each vector variation. The coefficient of each vector variation is its partial. 6(T)(-l T * 6R+I d 6 OR i | J For the case when a vector partial is being taken of a vector quantity a similar procedure can be used. In this case the result would be a dyadic. Obtaining this dyadic explicitly, can usually be avoided. In each case a vector will be dotted with it. To prevent the formation of a dyadic, the dot product should be done first, with the entire variation. After this the coefficients of the vector variations can be collected. For example: X k — R + - *R + d 1.8. Moving Frames of Reference In rigid body dynamics and in flexible body dynamics, moving reference frames are used to follow the motion of each body. The properties of the bodies are more easily defined in terms of these frames. However, the equations of motion require

14 derivatives measured in an inertial frame. For this reason, the inertial derivatives must be expressed in terms of the motion of the frame, and motion relative to the frame. In order to express derivatives measured in different frames a superscript notation will be used. The differential of vector v for example, measured in frame f is expressed as odv. The differential measured in frame b is expressed: ()dv. In the following development f will stand for the fixed inertia frame, and b stands for the body fixed frame. Consider the vectors defined in Figure 1.1. Since P = R + q an infinitesimal change in P measured in the fixed frame (dP is equal to the change in R plus the R ~ Inertial Frame z / Inertial Frame Figure 1.1 Vector definitions for a moving frame

15 change in q, both measured in the fixed frame. ()np ndR + ('d7 ( 1..1) (dR measures the change in position of the body frame. The change in 7 is made up of a change due to the change in orientation of the body frame, and the change of 7 relative to the frame. From Euler's theorem, any change in orientation can be represented as a rotation about a certain axis through a certain angle. An infinitesimal change in orientation will be denoted by dr. It is a vector in the direction of the axis of rotation with a magnitude of the angle of rotation. The change in q due to this rotation is therefore: dr X 7. The change in 7 measured relative to the body frame is ()dq. The total change in q is: (nd0 = <d X g+ (b)d7 (1.6.2) The total change in P becomes: dP = (dR + d' X q+ (b)dq (1.6.3) If these equations are divided by an infinitesimal increment in time they become dt dt X dt and OdP _(dR di r (&)dq dt dti d X q+ "di or q-= 2 X q +q (1.6.4) and P==R+ xq+ q (1.6.5) A dot over a vector means the time derivative measured relative to the fixed frame (dd |. An open circle over a vector means the derivative measured relative to the

18 body frame I )d. Equation 1.6.4 defines i as ~dti ~ ~ - J f i- -_ (1.6.6) dt and is 0 called the angular velocity. One of the properties that the angular velocity has is that its derivative in the fixed frame is the same as its derivative in the body frame. This can be seen by substituting f in Equation 1.6.4. n = n xn +n n = fi (1.6.7) One useful identity is that the derivative of a scalar product is the same no matter what frame it is measured in. (Ed (-< b) =~- ~ a e + a Od a- ' a b- + a -b = I(n x 3+ + — t '* + ( x b)+ -- * + a *-b + + a -n x b -a b + a b 1.7. Rotational Coordinates So far, only the change in orientation of a local frame has been investigated. An infinitesimal change in orientation can be represented as a vector; the orientation itself, in general, cannot be. To specify the orientation of a frame, a set of rotational coordinates are used. In planar problems all that is needed is one angle. For three dimensional problems, three angular coordinates are usually used. In deriving the equations of motion, it is easiest to formulate the equations in terms of the angular velocity vector fl,.which is the rate of change of angular orienta

17 tion. What is needed then, is a link between the angular velocity n and the angular coordinates. The equations for this link are called the kinematic differential equations. In addition to this, the transformation matrix A bf needs to be defined in terms of the angular coordinates. For planar problems the situation is simple. The bodies always rotate around the z axis (perpendicular to the plane). The angular velocity vector has only one component: f,. That is, f = [0 0 of n,'. The angular coordinate normally chosen is the angle the local frame has rotated around the z axis (positive according to the right hand rule). See Figure 1.2. The kinematic differential equation is simple 1n, = (1.7.1) The transformation matrix from inertial coordinates to local coordinates is: cosG sine 0 A f -sine cos 0O (1.7.2) 0 0 1 9Y 9 / Figure 1.2 Planar rotation

18 In three dimensions, the situation is more complicated. But, since the equations of motion will be derived in terms of the angular velocity, they are fairly independent of the particular choice of angular coordinates. As long as the kinematic differential equations and transformation matrix for a set of angular coordinates are known, any set can be used. There are various possibilities. The more well known include: 1) Euler Angles 2) Bryant Angles 3) Euler Parameters Both Euler Angles and Bryant angles use sets of three successive rotations. Two difficulties that both these choices share are: 1) The kinematic differential equations are unsolvable at particular orientations. 2) Both the kinematic differential equations and the transformation matrix are complicated trigonometric functions of these angles. Evaluating these functions can be a significant part of the cost of computation during a simulation. The Euler parameter approach does not share these difficulties. The kinematic differential equations for them are always well defined and are simple linear functions. The transformation matrix is a simple quadratic function of the parameters. Its disadvantage is that it has four parameters with an added constraint function. This is not much of a problem for the implicit formulation, because it is suited to including algebraic constraint equations. For the explicit formulation, extra bookkeeping is required to make sure the constraint is always satisfied. But the numerical advantages far outweigh this extra bookkeeping. In keeping with these arguments, the Euler parameters were adopted for three dimensional simulations. The physical significance of the Euler parameters is fairly simple. From Euler's Theorem any orientation of the local frame can be obtained from the fixed frame by a single rotation around a particular axis (see Figure 1.3). Let this axis be defined by the unit vector i. Let the angle of rotation around this axis ( measured positive

19 Figure 1.3 Rotation around an axis according to the right hand rule) be 6. Then the Euler parameters are related to 0 and u by Po = cos-, and p = u sin-. Since u is the axis of rotation it has the same coordinates in each frame. Therefore p has the same coordinates in each frame also. The constraint equation relating the parameters comes from 2 - - -~ 2. ~ 1 p +-C2 u Ps - i COS2 + sin- 8 p P P 2 2 so po2 + 1 The relationships for Euler Parameters are well known and are simply repeated here (see Reference 22). The kinematic differential equation can be formulated in two forms. The following form is suitable for an implicit integrator:

20 fIl -Pi Po Ps -P2 fS2 - 2 -P2 -Ps Po P1 P s -P3s P2 -PI Po with o 2+p2 +P22 + p 2 = 1 or in matrix notation Equation 1.7.3 is (1.7.3) _ = p (1.7.4) where -P 1 Po P s -P 2 B = 2 -P2 -Ps P. P1 (1.7.5) -Ps P2 -Pi Pa For an explicit integrator, the equations are solved explicitly for the highest derivatives by including the derivative of the constraint equation and inverting the resulting equations. This results in * \'0 p. P I 1., P2 2 02 PS 123 -nl 0 -23 n2 -n2 0 -1n, -nf In P 2 0 pJ IPi1 (1.7.6) or The transformation matrix from inertial coordinates to local coordinates is: inertial coordinates to local coordinates is: 2(po2+p I)-1 Abf = 2(p 1PPoPs) 2(p Ips+Po P2) 2(p lP 2+P oP s) 2(po2+p 2-P1 2(p2p S-Po P 1) 2(p p s-po p 2) 2(p 2p +p, p 1) 2(po2+p 2)-1 (1.7.7) 1.8. Numerical Integration The equations of motion arising from both rigid body dynamics and flexible body dynamics are much too complicated to be solved analytically. Therefore a numerical integration procedure must be used to solve the differential equations. There are many

21 different types of integrators. Each type has different characteristics as well as strengths and weaknesses. One of the most obvious differences between integration procedures is the form that it requires the equations to be in. Most integrators require that the equations be in first order form, explicitly solved for the highest derivatives. Y= [(y,t) Less common are those procedures intended for second order equations of the form My + Cy+_K y= t) This form of the equations arises naturally in structural dynamics for finite element equations. Another form which some integrators require is implicit. Here the equations are written as implicit functions. _, I, t) = 0 In this case the implicit equations do not necessarily have to be differential equations, some can be algebraic equations. The entire set is a mixed set of differential and algebraic equations. This form allows the analyst more freedom in formulating the equations. The second characteristic which is important in the context of flexible mechanisms is stability. Because the numerical integrator is approximately solving the equations, the nature of the approximations may cause instabilities in certain situations. When this happens, a constant step size integrator may produce wildly erratic results. A variable step size integrator, on the other hand, will cut its step size and become very inefficient. The situation which brings on this instability is commonly called "stiffness". Stiffness occurs when the solution to the problem has small derivatives, but there is a solution nearby which has large derivatives. In this case, error generated

22 while following the low derivative solution may cause a jump to the high derivative solution. For a constant stepsize integrator, the high derivatives cause wildly erratic solutions. A variable step size integrator, will cut its step size to accurately follow this nearby solution. A stable integrator on the other hand does not produce errors in the direction of the offensive nearby solution. This is not to say that it does not produce errors, only that the errors it produces do not excite solutions with higher derivatives. Consequently it can take large step sizes and. still accurately follow the desired solution. In the case of flexible mechanisms, stiffness often occurs. The flexible members usually introduce high frequency components. At times these components are excited and important and they should be followed meticulously. At other times the damping in the system causes these components to die out. At these times the system is stiff. The solution has small derivatives, but a nearby one has large derivatives. For this case a stiff integrator is desired. One final situation that arises in flexible mechanisms also gives trouble. This is the case where some high frequency components are excited, but are so small as to be negligible. In fact, the magnitude of the component may be smaller than the error tolerance. Ideally, these components should be ignored and the step size of the algorithms suited to follow components whose magnitude is larger than the error tolerance. At the present time, unfortunately, even a stiff integrator will not ignore these components until they are virtually zero. The current way to ensure that they become zero is to introduce damping. This drives these solutions eventually to zero where a stiff integrator will ignore them.

23 One requirement that a stiff integrator has, that a non-stiff integrator does not have, is that it requires a Jacobian matrix. That is, it needs a matrix of partial derivatives (or some approximate form) as part of the numerical solution process. This is not always easily obtained. In fact for explicit methods, it is very difficult to obtain it. The reason is that the equations of motion for mechanisms are highly non-linear. If they are put in explicit form, obtaining an explicit Jacobian is a formidable task. The alternative to this, numerical differencing, is numerically intensive and dominates the cost of simulation in larger problems. In order to understand how all these factors relate to the integration of mechanism equations a few remarks must be made about the characteristics of these equations. First of all, no one has found a way to put the equations for general mechanisms into explicit form without forming and decomposing a matrix. Secondly, for the implicit methods, a matrix must be formed and decomposed as well. The matrix which must be decomposed in this case is larger and more highly populated. Generally, it has been found that the non-stiff explicit methods are faster than the implicit stiff methods for non-stiff problems. The reverse is also true. For stiff problems, stiff methods are faster. Generally, for mechanisms containing only rigid bodies and where impact or other exotic effects are not present, the equations are not stiff. For these, use of the explicit methods will give more efficient solution of the equations. For mechanisms containing flexible members, or exotic forces the stiff methods will outperform the non-stiff methods (as long as there is sufficient damping present). If there is insufficient damping to drive the stiff components to zero, the explicit methods will outperform the stiff methods. For these reasons, the equations are formulated in forms suitable to both of these methods of integration.

24 For explicit methods the equations are put into first order form and explicitly solved for the first derivatives. Since most of the equations of dynamics are naturally second order this involves introducing new velocity variables and simple relationships such as: RI -= V and 1 - N1 The rotation equations, however, are best formed directly in terms of the angular velocity and will already be first order in angular velocity. The relationship between angular velocity and angular coordinates is more complicated, but can still be formulated easily in explicit form using Equation 1.7.6. P - QP _Q_ where Q =Q()M The equations which make up these definitions can be solved independently of the equations for the dynamics. The equations of dynamics will be presented in mixed second and first order form with the understanding that the equations just previously mentioned must be applied to actually convert them to first order form (except, of course, the angular velocity). For stiff problems the implicit form of the equations are presented. The reason for using the implicit form for a stiff method is that a Jacobian matrix is needed. It is much easier to obtain the matrix of partial derivatives for equations in implicit form. It is also easier to produce a Jacobian matrix which is sparse. One of the most important aspects of this formulation is the sparsity of the Jacobian matrix. If the Jacobian matrix has a high percentage of zero locations to non-zero locations (sparsity), the decomposition and back substitution can be performed efficiently. Since these two steps are the most time consuming, this dramatically reduces the overall effort required

25 to solve the equations. The sparsity of the Jacobian matrix is take advantage of by using subroutines especially adapted to sparse matrices. One way to keep the matrix sparse and at the same time, reduce the number of equations, would be to allow the equations to be in second order form, as well as first order and algebraic form. Examination of the derivation for the Gear integrator [27] showed that this indeed was possible. The original derivation done by Gear does not restrict the algorithm to first order form. It was done for arbitrary order. When it was applied, it was always applied to first order equations. Upon recognizing this, new integration subroutines were written which included second order equations - based on this same derivation. It is included as Appendix A. These subroutines will integrate mixed sets of second order, first order, and algebraic equations. The only major operational difference between these subroutines and the former is that the Jacobian has an extra component. Previously, the Jacobian was the sum of two Jacobians. [[OF 1r O[F1 II [L T j- + where 71 is a scale factor computed by the integration procedure. The new Jacobian matrix is the sum of three individual Jacobian matrices ~2 = a-' 7 + [ I + _ G2=[ r 9OTr J 1 9Fo 1 [ Lor ]r where both 7i and -72 are scale factors calculated by the integration procedure. For second order equations this preserves the sparsity of the Jacobian matrix. Non zero terms tend to land on top of each other when the Jacobians are added together.

CHAPTER 1 THIE SINGLE RIGID BODY 2.1. Introduction The field of rigid body dynamics is not new. It has been studied and expanded for centuries, thus many of the equations that will be derived in this chapter are very old. One of the reasons that this chapter is included is to provide a review of the equations and principles of rigid body dynamics. From this review many of the terms in the flexible body equations (which are more complicated) will be more easily recognized. It will also provide an opportunity to become familiar with the vector notation used throughout this dissertation. In addition some of the strategies used in simplifying the equations of motion for rigid bodies can be carried over to flexible bodies. 2.2. The Local Reference Frame and Vector Definitions In order to locate a rigid body in space, a reference frame is attached to the body. A vector from the origin of an inertial reference frame to the origin of this local reference frame is all that is needed to specify the location of the local reference frame. This vector will be R throughout the dissertation (see Figure 2.1). Only the position, not the orientation is specified by R. A set of angles (usually three) is needed to specify the orientation of of the local frame. These angles cannot 28

27 usually be represented as a vector (planar problems are an exception). The inability for the orientation of a body to be represented as a vector is a source of many complications in dynamics. Many of these problems can be circumvented, however. Although a finite rotation cannot be represented as a vector, an infinitesimal rotation can be. The infinitesimal rotation will be denoted dr (notice the arrow is over the entire quantity). The rate of change of orientation can also be represented as a vector. dz This is the angular velocity and will be denoted 1 throughout (s -= d) dt The.position of every-mass particle in the rigid body can be described relative to the local frame. The vector from the origin of the local frame to an arbitrary mass particle is denoted with a lower case 7 (see Figure 2.1). With respect to the local frame, the position of every mass particle is invariant. riz I / b,~~~~~~~ Figure 2.1 Vector definitions for a rigid body

28 The position P of every mass particle can be represented as the sum of two vectors: R - which is common to all mass particles, and r - which is the location of each particle relative to the local frame and depends on the orientation of the frame. P = R+r (2.2.1) An infinitesimal change in position P with respect to the inertial frame is (Equation 1.8.1): d;p = cdR + 0)d Since 7 is invariant in the local frame ( (b)dr 0 ), Equation 1.6.2 becomes: O(d7 = dr X r (2.1.2) and ')dP (dR + dr X r (2.1.3) where drf X r is the change in position r due to an infinitesimal rotation of the frame. In terms of time rates of change, these two equations become the velocity equations: - n x 7 (2.1.4) P -R+n xr (2.1.5) where P is the absolute velocity of a mass particle R is the translational velocity of local frame n is the angular velocity of local frame r is the position of a particle relative to the local frame. 2.3. Linear Momentum The linear momentum for a rigid body is the sum of all the individual linear momenta of the mass particles. This can be written as

29 P i P dm (2.3.1) where P is the total linear momentum of the rigid body, Pdm is the contribution of an infinitesimal mass particle, and the integral is over the entire mass of the body. The velocity P in terms of the local reference frame was given by Equation 2.2.5. Substituting this for P gives the following expression: P = f(R + X dm or P = Rdm + If x rdm Since R is common to all of the mass particles, it is constant with respect to the integral and can be moved outside the integral. The same is also true of fl. This yields P-= Rfdm +Ix FJdm These integrals can now be evaluated in the local reference frame. The sum of all the mass particles is just the total mass. The integral over the mass weighted position vectors gives the center of the mass vector. Let M = fdm (Total mass) (2.3.2) MrC = f7dm (Vector to C.M.) (2.3.3) The linear momentum becomes simply: P = MR + f X Mr (2.3.4) From equation 2.2.4. we can see that f x rcm is the absolute velocity of the center of mass (rcm) relative to the frame. So this equation states that the linear momentum of the rigid body is the sum of the translational momentum of the local frame, plus the

30 linear momentum of the center of mass relative to the frame. This is a statement of the well known fact that the translational momentum of a collection of particles is the translational momentum of the center of mass. If the frame origin was moved to coincide with the center of mass, rc,, would be zero and Equation 2.3.4 would become P - MR (2.3.5) 2.4. Angular Momentum The angular momentum of a rigid body is the sum of all the individual angular momenta of its mass particles. The angular momentum of a mass particle is always defined relative to some reference point. It is defined to be the cross product of the vector S from this reference point and the linear momentum of the mass particle -'SHy =~ S- X'mP "(2.4.1) H, =Ms x mP its own frame of reference. In this case the vector r is used for S. H = | 7 X Pdm (2.4.2) Again using Equation 2.2.5 for P H = Jfx (R+n x)dm Rearranging, H = f(-R x )dm + FJ X (f X 7)dm The R can be moved outside the integral in the first term and simplified

31 f(-R X i)dm = -R X fJdm — R x Mrcm The second term can be simplified by expressing the double cross product as a dyadic product (refer to Equation 1.4.16), fr x (fl x7)dm = f[(- )I- 7-] * ' dm Omega can be moved outside the integral, and it becomes [I(.)7r -rldm -n Now define J as the term in the integral - fJ(7 -r 7)- rldm (2.4.3) Jf x (f x 7)dm = j.n (2.4.4) The second order tensor J, is called the inertia tensor. It is a symmetric tensor and when its coordinates are found in some reference frame it can be represented as a 3 x 3 symmetric matrix. (fry2 + r2)dm -jfr ry dm -Jr, r, dm _-== f(r,2+r,2)dm -r r, dm (2.4.5) SYM f(r,+ry2)dm The angular momentum becomes: H -— R X Mrm + J (2.4.8) -rc x MR+J 'n Here it can be seen that the angular momentum of a rigid body around its reference frame's origin is the sum of two parts. The first term is the angular momentum of the center of mass around the frame origin. The second term is the angular momentum of the frame due to its angular velocity. Again it can be seen that if the origin of the frame coincides with the center of mass, the first term will be zero and a simpler equation results:

32 H =J- (2.4.7) 2.5. Kinetic Energy The kinetic energy of a mass particle is defined to be 1 x XTP m(P - P) where Tp is the kinetic energy of the particle, m is the mass of the particle, and P is the velocity of the particle in an inertial frame. The kinetic energy of the rigid body is therefore: T -=i PP dm (2.5.1) Expanding P from Equation 2.2.5 2T ~ JIR + fn xr.(R + fI X rdm T = - (R -R)dm + f(R - n Xrdm + 2 J(n Xr) (f Xrdm (2.5.2) The first term is simplified readily by moving R -R outside the integral: 2 f(R R)dm 2 (R- R)fdm 1 J^4 -4 - M(R- R) The second term is simplified similarly f -4 L -~ J(R ( Xi)dm = R X dm2.5.3) (2.5.3) -R' x Mrc The third term is simplified by recognizing that the dot and cross can be interchanged in this mixed product (let - = n x r and it becomes a triple scalar product) 2On mo g ouside te ingral this. f x( xbeom On moving fl outside the integral this becomes

33 i E. 'F 7x( n x7)ldm This can be further simplified by expressing the remaining integral in terms of the inertia tensor (Equation 2.4.8) 12 ~ n- n (2.5.4) 2 The kinetic energy for the entire body becomes 4 ~ * - 4 -T=- M(R * R) + R XMcm + n.J n (2.5.5) 2 2 The kinetic energy of the rigid body has three parts. The first part is the kinetic energy due to translation of the local frame. The second part is kinetic energy of the center of mass relative to the frame. The third part is kinetic energy due to the rotation of the frame. Notice again that the equation can be simplified by choosing the frame to be situated at the center of mass, and oriented to diagonalize the inertia tensor. In this case, the kinetic energy becomes: T = M(R R) + n * - n (2.5.6) 2.8. Application of the Fundamental Equation The fundamental equation for dynamics given in Chapter I can now be applied to the single rigid body. This will give its equations of motion. For now, the terms involving forces and constraints will not be considered. They will be added subsequently in Chapter VI. The variations to be used will be the vector variations of the position of the local frame 5R, and its orientation dr. The corresponding velocity variations are SR and 6l. In order to get the partial derivatives of the kinetic energy with respect to generalized velocities, the variaion of the kinetic energy with respect to velocity will Oq

34 be performed (this is done while holding position constant). 6{lj = -2 M(R R)+Rn x Mho M + 2 nJ-.J (2.6.1) The first term in this equation is evaluated using the product rule 1 4 JL 1 + A 1 4 4 6( MR Rj 2 MER *R + MR * 6R The ordering of the terms in the dot product is immaterial so the first term becomes: 6j 2MR-R =6R-MR (2.6.2) The second term has both R and f in it so it will give variations for both of these. R-n X Mrm = -R f XMr m + R * Sf XMMr Rearranging yields OR a x Micm}6R *inxMic, +6ni Mrcm xR (2.8.3) The third term has only f in it. } 2 n + Since J is a symmetric tensor, the order of the dot products is immaterial and this term becomes 6f 2 n7 -n =s-7n n (2.6.4) The velocity variation on the kinetic energy yields: { T}= R ~(MR + n xMrc)+ -(J -n + Mc xR) (2.6.5) The terms in the parenthesis dotted with each variation is the partial for that velocity.

35 T = MR + n x Mrm (2.8.8) aOR OT 4 -I =J n + Mr*c xR (2.6.7) an These partials are called the generalized momenta. For this choice of coordinates they turn out to be the linear momentum of the rigid body and the angular momentum, respectively. This can be seen by comparing them with Equations 2.3.5 and 2.4.5. The time derivative of the linear momentum follows: (P) (MR + n XMr 4 (f )d=-MR + dt (n xM,) dt The time derivative of the second term is best expressed in the local frame. This time derivative becomes < (n x MC,,) = n x (n x, )+ ( xM ) dt dt The product rule can be used again dd (n xMr, ) n xMrc,, + n x MrC dt The second term is zero because rm is invariant in the local frame. The time derivative of the linear momentum becomes ~~~~~d = =0 A 4 -+~~~(2.6.8) d (P) = MR + x (n XMr, )+ n Mr (2.6.8) The time derivative of the angular momentum is next: (/)d -* 0()d " -' L' d (H) d (J - n + Mc,, XR) The first term is most naturally measured in the local frame Tedrvtvoftes ond ter i e-wi 4h p-o* u rule: i - - - A d (I *n) = n x J n + J n The derivative of the second term is done with the product rule:

38 0I d '* Z*. L* td (Mr-c X R) (Rem ) X R + Mrm XR dt From Equation 2.2.4 r= 0 x Od (r x - dt (AMrm X R) = ( X rm) X R + M,, X R dt Finally the time derivative of the angular momentum is (Ond - (-2 ) 6 () (H)= J-n + x i n + (n x M7r, ) xR + M, xR (2.6.9) dt Since the kinetic energy in this case does not depend on the position or orienta81T a8T tion of the body, the partials with respect to them d, — are zero. Lumping all of the remaining terms into f (force) and t (torque) gives the following expression for the fundamental equation: 8R 'MR + n x(n XMrn ) + n X,= -f ) + 6. -. n + n xJ.n + (n xMrcm ) R+ M +. xR-t (2.-.10) Where f and t represent the additional terms which will be arrived at in later chapters. If the frame is at the center of mass the fundamental equation becomes R MR- f +6.?- n+n xJa-t] =- o (2.B.l) Since the variations are arbitrary the terms in the brackets must be zero. This gives the equations: 14 MR- f-= 0.n + n x J.n-t= (2.6.12)

37 2.7. Matrix Form of the Equations The purpose of this section is to transform the vector equations derived in the preceding section to matrix form. This will be done with the aid of the matrix representation of vectors and tensors given in Section 1.4. The vector R has derivatives measured in the inertial frame. If it is resolved in the inertial frame, the derivative of R is simply expressed in terms of the derivative of T..T its components. So R = R e( ). Likewise, 6R -= 6R r (, and R R = R ). The vector f will also be resolved in the inertial frame. The coordinate matrices for these vectors are: - R - R ' -2 -/2 R = RX, R = IRJ R =2 |R2) 1 J = IP Rj R3 R3 J If we choose the principal axis of inertia for the body as the body fixed frame, the coordinate matrix for J is diagonal and becomes - J2 (2.7.1) J3 Also if we locate the frame at the center of mass of the body the equations of motion themselves are much simpler. Since these equations are the ones actually programmed, they are done explicitly here. Equation 2.6.11 is repeated here with the f x replaced with its tensor equivalent fl. R. MR-t J4-r -n 4- n 4- -} + r. ' f1+fl. n. n-t -= 0

38 Replacing each vector with the appropriate matrix equivalent yields: SR T^^n M^n~R-tVI'Ir 6 1+(J). (M)R - e /~7}f + 6.I~ -(b),. + ev(r n Io)( ~ F(O)r J e(rnl -e' l} =0 Now the dot products of the unit vector matrices can be carried out. In this case each product is between unit vectors of the same base. An identity matrix is the result of each dot product from Equation 1.4. The equations reduce to: - MR - The variations multiplying the braces are arbitrary, so the terms in the braces must equal zero. MR- - 0O (2.7.2) IJ+nJn-t = o For an explicit integrator, the equations must be solved for the highest derivatives in each variable. The highest derivative in R is R. The highest derivative in n is f. These will be grouped on the left and the remaining lower derivative terms on the right. MI f - I - nJ+ (2.7.3) This can be written out explicitly.

39 M f - -I - - -2 (iJZ-J)fl2 + tI J2 03s (1-J%2)nn02 + ts I Js I=nn,+ j Notice that the matrix multiplying the highest order terms is diagonal and constant. To solve this set of equations, each equation is multiplied by the reciprocal of its diagonal element. Together with these equations, the kinematics differential equations must also be integrated. These equations were given as Equation 1.7.6. When there is more than one body and there are constraints involved, the dynamic equations will involve Lagrange multipliers. These must be determined simultaneously with the highest derivatives. The total system of equations is more complicated, but the simple equations for each body reduces the effort and complexity in solving the system. This will be treated in more detail in Chapter VI. For a stiff integrator, the equations are written in implicit form. In addition to the dynamics equations, kinematic differential equations must be included and integrated simultaneously. For the implicit integrator the appropriate form (Equation 1.7.3) is n B p with P r P 1 Together, the equations for a rigid body become MR -f = O ATk'- - 0 Tr E-1=0 Bp - n 0 IJ n + n J n - t o _ m m _ _ _

40 For the implicit integrator the partials of these equations must be calculated. The Jacobian matrix is the sum of three individual Jacobian matrices (Equation 1.8.2). FG2 - 72 + aF1 0- Oz xzm In the partial derivatives that follow, the partials corresponding to L and t will be ignored - they will be calculated separately in Chapter VI and added there. Most of the partials are very simple, only two are not obvious: 1) For 1f J n the partials could be calculated by matrix manipulation, but it is simpler to write out this expression and take its partials manually. (Js3- 2) f2 fSl n( nl (JI-3) - s J Il (J2 Jl) fix fl We define C to be the matrix of partials: 0 (3 - 12)nA (J3 - 2)n21 C (11- I3)2 0 (i - 3Is)1l (2.7.4) (12 - J)r2 (J2 - J1)n 0o 2) B p The partials with respect to p are simple: - B. The partials with ap respect to p are also simple but not at all obvious at first glance: aB = B. p - This can be seem by writing out the expression in Equation 1.7.3 and taking the partial derivatives. For the system of equations: MR-f O rT -2=0 Bp- =o0 (2.7.5) f J fl- t= _ _ _

41 The total Jacobian for the rigid body is: 2p_ B71- B -I 7i + C (2.7.6) where -PI Pe Ps -P2 B=2 -P2 -Ps P, P -Ps P2 -Pi P. o (J - J2)aS (3- J2)2 C- (JI- J3)n (0- ) (J2- Jl)f2 (J2- J1)l1 0~ J1 i 12 J3

CHAPTER III THE SINGLE FLEXIBLE BODY 3.1. Introduction In this section the equations of motion for a linearly elastic body will be derived. They will be derived in general first, so that all of the possible terms will appear in the equations. By using different frames of reference, different deformation shapes, and neglecting small terms, various simplifications of the equations can be made. The pros and cons of the various types of equations will be discussed and examples done to illustrate each one. The inertia and energy properties of the flexible body will be derived first, so that the results can be compared to those for a rigid body. This will cause new bodydependent constants to be defined and the meaning of these will be examined. The fundamental equation will then be applied and the equations of motion will result. 3.2. Vector Definitions and Deformation Coordinates As in the rigid body case, a local reference frame is associated with the body. Because this reference frame follows the body in some way, the deformations will be measured in this frame. This will preserve their linear character in spite of large overall motions of the body. 42

43 The definitions for the vectors used in the derivation can be seen in Figure 3.1. The vector R locates the body reference frame relative to the inertial origin. The undeformed shape of the body is represented by the dashed outline. The vectors r, are the positions in the local reference frame of all the mass particles before the body deforms. The vectors ui represent the change in position of each mass particle due to deformation. The vectors q are the instantaneous position of each mass particle relative to the local frame. The instantaneous position of any mass particle is P, and is the sum of the R and q. In general, there are an infinite number of i vectors. The deformation u depends on its location in the body (r), and it depends on time (i). It is both a spatial and a time dependent quantity [u - '(rt)j. Therefore an assumption is made that a separation of variables can be performed: / / I I \ o I I i I! R Figure 3.1 Vector definitions for a flexible body

44 N - u( r, t) -E (-r7)p () or in index notation urt) -,) (7) (t) Here it is assumed that the deformation of the body can be represented as the sum of shape functions, (7), multiplied by corresponding time dependent deformation coordinates. At this point we have assumed that the deformation is a linear function of deformation coordinates. All of these relationships are summarized in the following equations P R + (3.2.1) r + (3.2.2) (7,-t) =, (7)~, (t) (3.2.3) The differentials and velocity equations follow directly from the identities of Section 1.6: (dP (dR + (dq (3.2.4) (d q = d i Xq+ (b)dq (3.2.5) This time the body fixed change in q ( (b)d is not zero as it was for rigid bodies. From Equation 3.2.2: (b)dq (bd + ()du. Since 7 is invariant in the body frame its body fixed differential is zero. This equation becomes (b)d - (b)du (3.2.6) From Equation 3.2.3, the body fixed change in deformation becomes (b)du =,i () (b)d, (t) + (b)d^1 (7)t (t). Since the shape functions ti (7) are invariant in the body fixed frame the second term is zero, and this equation simplifies to: ()du; = (b)d r (3.2.7) Altogether, the differential in P is

45 (ndP (dR + dir X + i (b)dp, (3.2.8) The velocity equations corresponding to these differential equations are: P R+q (3.2.9) _~ n x~q + (3.2.10) e e e - q,, (3.2.11) = - '.i (3.2.12) P- R + nx + l (3.2.13) 3.3. Linear Momentum The linear momentum of a flexible body is again the sum over the entire body of the individual momenta. P - fPdm Expanding P from Equation 3.13 yields: P J(R + ni Xq +l jA )dm (3.3.1) The first term becomes again: fR dm - MR The second term needs further expansion. First the angular velocity can be moved outside the integral. Jn x q dm = f X fqdm. Secondly, q can be expressed as r+ u, where = iJA: f 7dm - fJ(; +,, t )dm The first of these (frdm) is just Mr,m as before. The second term introduces a new constant. Since /,j is purely a time dependent quantity it can be moved outside the integral over the body: fJ j, dm = Jf t dm s, The quantity remaining inside the

46 integral is purely spatial, so that integral can be performed for each mode shape i, once, before the simulation begins. Let us define M~ f Ij dm (3.3.2) then fI'dm becomes f7dm M(r7c + u8sj) (3.3.3) and fn x qdm becomes fl X M(rcm + A. ) -. 0 For the third term in the momentum equation (fi, i, dm) the A, can be moved outside the integral and the entire quantity becomes MsA i. The linear momentum of a flexible body becomes P = MR+ n XM(r + sj p) + Mis (33.4) The expression for the linear momentum of the flexible body has two new terms in comparison with the expression for rigid body. Both of these involved the deformation coordinates. The term i Aj, is simply the distance the center of mass has moved from its original undeformed location. Each i is the amount the center of mass is moved by the ja deformation coordinate. The term MZ, A, is the linear momentum due to the velocity of the center of mass within the body frame (a sloshing motion relative to the body frame). It may be noticed again at this point, that if the body frame of reference was located at the instantaneous center of mass (rom = 0, and ' = 0 for all i) the linear momentum would again be the simple expression:

47 P = MR 3.4. Angular Momentum The angular momentum of the flexible body around the origin of the body frame is: H =Ja x Pdm (3.4.1) If P is expanded using Equation 3.2.13, this becomes H =- fJ X (R + n Xq + s, M,, dm The first term can be expanded using previously defined quantities: Iq x Rdm =J(- + f(+ J ij)dm X R -= (,. +ip,) X MR The second term has a double cross product in it. It can be expressed in terms of a tensor (cf. Equations 1.4.15, 1.4.16) Jf x ( x q)dm = f[()q — q-qgdm* n where the flexible inertia tensor is: Jf[( I - lqqdm (3.4.2) This inertia tensor is not a constant as it was in the case of a rigid body. 7 can be expanded in terms of 7 and u.

48 J = f(i(+ [( u (7+i)I- (r+u-)(7r+ui) dm = Jf[(r 7) - 7idm + f[2(. * )7- (ru + -)ldm + Jl(( *)t?)/ - dm The first term is a constant and is the same as the rigid body inertia tensor. It will be defined to be O: JD= If(7r )7- drm (3.4.3) The second term is a function of the deformation. It can be reduced by replacing u with /Ai,.e f12(r * (;)I- (ru + ur dm = p)2(7r i )7- (t i + )] dm = (2(r i )I- (ti +-ii7)ldm j The quantity inside the integral is now completely spatial, so that it can be carried out 4-, once for each i. Let us define it as Jl, 7l = [12(7 r* i )I- ( +, 7])ldm (3.4.4) This part of the inertia tensor is J1, i. The coordinate matrix for this tensor is again symmetric; 2f (ry y, + r, *a )dm -f (rt, 'yi + ry i )dm - f(rs, i, + r,,, )dm] J1 = [| 2J(r. y, + r, qn )dm - J(ry, + r, ly )dm (3.4.5) SYM (2f(r, ti' + ry,i )dm) The third part of the inertia tensor is second order in displacement. Expanding u in terms of deformation coordinates gives: 4-u -dm * f[( U. )-UI -_]dm =- [ '(/ ' jA)I-, Jdm

49,= f JI( k)I - j, 1k Jdm A This term becomes 4," IJ2,t it where J2,t [(itj.* )I -,j Pk 1dm (3.4.6) Finally, the entire inertia tensor J becomes: I =J+ JO 1 + l i, + j.2,k (3.4.7) The angular momentum has become so far: H = (rc + ai i) X MR+J f + f(ii X )dm This last term can be expressed: f X tudm = (q x is, Idm = f(+ + mSmS) x A.rjdm = X r6dm A, + ^, Jf XI,dm t, Defining the first integral to be a, a, f Jr X,, dm (3.4.8) and the second to be G,, Gm - fn^ Xc dm, (3.4.9) this term becomes JF X;dm = (-, +, Gd"9)I' (3.4.10) Then the angular momentum becomes: H -= (r,, + 8l) XMR+J * + (a, + n Gme)IL,

50 Examination of the -various terms in the angular momentum will give insight into what each one means. The vectors 78 appear again in the term which represents the angular momentum of the center of mass around the frame origin. The term rm + 1 gi is simply the position of the center of mass relative to the frame origin. The angular momentum of body due to the frame's angular velocity is the second term. In this case the inertia tenser is not constant, but depends on the deformation of the body. As the shape of the body deforms the inertia tensor changes. The third and fourth term represent the angular momentum that the deformation has relative to the local frame. The first part of this is the contribution of each deformation mode independently. The a, constant measures the component of rotation for each deformation shape. This, times the velocity of the deformation, gives an angular momentum. The second part of the deformation angular momentum is due to coupling between deformation coordinates. The G., constant measures this "cross coupling", or "cross momentum". 3.5. Kinetic Energy The kinetic energy of the flexible body is derived from the same statement as the rigid body T =- JP i P.dm If P is expanded using Equation 3.2.9, this becomes T = — - 1. 4 *1 2 (R - R)M + R * jqdm + - jqqdm The second term can be expressed as

51 R- qdm = Rjf[ Xq + isi Jdm 44 _~. o R- R xJqdm + R f* J i dm si From Equation 3.3.3. fqdm = M(r, + 7,). Likewise from Equation 3.3.2. f J dm =- MZi. So this term becomes R fJqdm - R' * n M(rm Xi ji) + R * i (3.5.1) The third term is expanded in the following way. From Equation 3.2.10 q fl x q +'. With this substituting, this term becomes 2f Iq dm = Jn q + ul [n xq' + uldm 2 J(n x ) * (0 x -)dm oa0~~~~~ (~(3.5.2) + fi xq uidm + 2 J * dm The first term in this expression can be written 2 J(n x) (n) x ) dm = 2 n. |q x (n x q-)!dm by interchanging dot and cross. The part in the bracket is recognized as J f so this 1 4-. - 4 — term becomes f J fl. It must be remembered that J is the flexible inertia tensor and is no longer a constant. The term Jf x qa u"idm is also done readily. n x q0m n 0- x dm f -. f From Equation 3.4.10, fJ X udm = (a, +,,, Gme ).e. This term becomes n (e + A, G,, )/,e. The last term in this series introduces a new quantity. Expanding it in terms of deformation coordinates:

52 2 * udm = -- A ' *pr fp, dm 1 < -~ -* ~ 2 I- 7k & fp dm p, The spatial integral in this case is called the mass matrix (when in matrix notation). Mk, J fk p, dm (3.5.3) This term becomes 2 -fu -dm 2 Ak MkP Ip (3.5.4) The total kinetic energy of the flexible body then is + n * n. + n ( + A G) I 9 + J-f.+Q-(,, +,Gm,,e)p, I ~ + 2A Mkt, The expression for kinetic energy is much more complicated than that for the rigid body. The new terms which were arrived at in the momentum derivations are present here also. In addition the mass matrix is involved in the kinetic energy. 3.8. Strain Energy For the rigid body case, the body is assumed inextensible so there is no energy exchange with the internal energy of the body. In the case of a flexible body, the body can store energy as strain energy and exchange it with kinetic energy as the system moves. For the flexible body then, characterization of the strain energy is necessary. The characterization for strain energy to be used here is a quadratic function of the deformation. This is the case for a linearly elastic structure. U = 2rw(:E:7(u) (3.6.1) where U is the strain energy, () is the strai n at a point, E is the elastici constant where U is the strain energy, 7(u) is the strain at a point, E is the elasticity constant,

53 and dV is an infinitesimal volume. The exact nature of the elasticity constant E will not be discussed here. It is discussed in equivalent forms in finite element, and strength of materials texts. The more important, and recognizable quantity needed is the stiffness matrix for a particular choice of deformation coordinates. From Equation 3.2.3: u = Hi, A,. The strain energy becomes U = r )iE:-{^) dV 2 ~ The term within the integral is purely spatial and is called the stiffness matrix. U -= -, K;,, Cj (3.8.2) where K,, = f(i): E: (1j) dV (3.6.3) The matrix [K,, ]is called the stiffness matrix. 3.7. Application of the Fundamental Equation We are now in a position to derive the equations of motion for a single flexible body by applying the fundamental equation from Chapter I. The variations will again include the vector variations of the position of the local frame dR and its orientation r6. In addition to these, variations in the deformation coordinates 6I,, are included. First, the partials of the kinetic energy will be taken with respect to the cT aT aT corresponding velocities (, T ). This is done by taking the variation of dR an.

54 the kinetic energy with respect to velocities, holding position constant. I -i = 6{ 2j MR R + R * XM., + si M) + MR 8j ij 1 ' - -. O + 2n *. * n + n * e)A, + 2 Mp A)p) The variation of the first term is 6( M R1= 2 M R + MR R 6R- MR The variation of the second term is R - n xM( + s8 = 6R -S x M (ra + 8 it ) - 4 -+.+ - =R'flx M(r-n + s~ll) =6R fn M(c, + st) + 6f (Ir +ipi) X MR The variation of the third term is: -f ^ - 0! A* -^ 0 0.L {MR *sj 6j -R.Msj /j + 6 j MR * sj The variation of the fourth term is: 1 - - - 1 — # 4 - - 12 Jf 2+ 2 -- 4-f -_ 1n j *n The variation of the fifth term is:

55 { (a + GffG )P } (ac + Am sG )t + sAC, * (a, + A, G, ) The variation of the last term is: 1 10 0 6( 2. IMkLP 1p = -b 2 t Mp 'p + 1 Mk p6tp Since the M,. is symmetric this becomes: a 0 atk Mk, CPp The velocity variation of the kinetic energy has become: 6{T)} = R [MR + n x M(rn + Ml pi) + Sj iJ + sn l + (r^ + 8,i ) x MR+( + ( mG+ a.) + Ak [Mkp A + MRS * k + n '(ak + Pm Gm )1 The terms in the square brackets are the corresponding partial derivatives. Examination of the partial derivative with respect to frame velocity shows that this is just the linear momentum of the flexible body, P (cf. Equation 3.3.4). OT -' --- =P (3.7.1) OR Likewise, the partial with respect to the angular velocity is just the angular momentum of the body around the frame origin, H (cf. Equation 3.4.11). 0T -' -~. H (3.7.2) Writing the fundamental equation in these terms yields the following: -" ( d 0 T c9 U - R' * d (P) -- + - a + *dt OR -OR + dir ( (H) -- } dt air r3

58 (t d T OT + t {"8 ( a- ) - ' } = 0 where f, t, ak represent generalized forces not explicitly included as yet. Examination of the kinetic energy expression reveals that it does not explicitly depend on the body's position or its orientation. This is also true for the expression for the strain energy - there is no dependence on frame position or orientation. ThereaOT aT O U fore the following partials are zero in the fundamental equation: -R.,, - oR 8' aR a8r a-U. Both the kinetic energy and strain energy do explicitly depend on the deformation coordinates. Therefore -T and U are not zero in general. da,,a. Now the fundamental equation has become: R {P - f} + rX - (H- t) { eod aT 86T au + +,k........ ffl;~k 0 dt isk bsk + o a*k where P MR + S X M(r + 8 1#) + MsJ Aj H J n + ( + 8i A.) x MR + (a. + AmGM )ce< 8T 7' -4 - - - * Me pi, + MR * sk + n *(akc + /i Gn ) ant In order to obtain the partials of the kinetic energy with respect to deformation, the variation of T is taken with velocities held constant. 5{T})= MR R + R - fl XM +sl)l)-+ MR Sjoj +2 n J n + n * (a, + Am G+ c )ic

57 10 2 At Mtp f, j The only terms containing deformations are the second, fourth, and fifth. The second term is: 6{R- n xM(7 + S, t )} This becomes 4 - _ MR. 0 X sl il The fourth term has the inertia tensor in it which is a function of deformation. 1 ' 1 -- *4 -14 6{- -f2l (2} 2 flS-J —f 2 ~2 # ++4 -J-= li, 6s, + J2,, 6o,p, + 1J20,,, p (J1 + 2J2,, i, )6 A, This term becomes: 0 * (Al, + J2zp Ap ) -ln bA The fifth term is S{f i (a, + A, Gm. )p}e f * Gm, Pe 6s OT Altogether -- becomes: Oti 4 1 4 _ MR fl xsi + fnl Gl,* + nl- (Jli + 2J2lp, p) n To obtain the partials of the strain energy with respect to deformations, its variation is taken 6U - =6{~1,K,, + 1 = 1 6,^sKCis +,- pK,, Sw,, U-2 2 Since K,, is symmetric this becomes: 6U = K, j, 6 SiThe partials ~u are then

58 Kk j l. To finish off the deformation equations and put them in their final form, the time derivative of their velocity partials must be taken. The time derivative of the partial with respect to deformation velocities is: (td.. _. - _ -.,di [M, Ap + MR * 8k + n (k + amGmk )1 The first term is a scalar so it becomes 00 The second term becomes MR sk + MR (a Xs) The third term can be expressed using Equation 1.6.8: n f(ak +, i, Gk) + fl Am, Ct Together this becomes: Mkp ip + MR sk + f (ak + Am Gmk) 1 -._ = - _ _ 0 + MR n x si + nl Gii m Finally all of the terms corresponding to the variation of a deformation coordinate have been derived. Collecting them all inside the braces yields 00 - - _ _ 6t Mkp t +MRK- sk + f (ak + m Gmk):.+.=+ + -- -_ 0 + MR (f X sk) + n ~ Gmk Am 4 -4 0 * 1 t -[MR- (f X k)+ fl Gk, u, + n (Jlk + 2J2,kpp) ) + Ka )

59 This can be simplified slightly by noticing that the terms involving MR * (0 Xs) -. -. 0 -.- ~. 0 cancel. One further simplification results from rewriting - (0 - Gk,,) as 0- G,, m and adding to the existing identical one. With these simplifications, the fundamental equation becomes ~R' P-f| 6R P- f + 6IH-t + 61 M, p + Kkj I + + MR sk 0 G 0 + f2 *(k + A Gk ) + 2fl 'G st - 2 fl (J + 2nk2 p) ) J O J. 0 -. 4.. -. ~ where P =- MR + fl XM('ca + si#1 ) + M sJ (j H J Q n + (rA, + si, ) x MR + (e, +. m G, ) At this point two approaches can be taken. The first is to leave the equations in this form. Here we have defined variables for the linear and angular momenta, and their derivatives appear in the fundamental equations. The second approach is to explicitly take the derivative of each momentum in the fundamental equation. One of the disadvantages for doing this is that the equations will get more complicated. Taking a derivative explicitly usually generates more terms. A disadvantage of the first method is that there are six extra equations to carry along.

CHAPTER IV SIMPLIFICATIONS FOR DIFFERENT FRAMES OF REFERENCE 4.1. Introduction The equations of motion for a flexible body derived in the preceding chapter are fairly complicated. They were derived for the general case, so that the simplifications which arise for different choices of reference frames and deformations coordinates can be compared. Just as the rigid body equations became much simpler from a suitable choice of a reference frame, the flexible body equations will also. Along with the reference frame choice, different types of deformation coordinates will also simplify the equations. In fact the two choices are interrelated. The two major types of reference frames to be discussed in this chapter are the locally attached reference frame and the floating reference frames. The floating reference frames will be seen to give the greatest simplifications to the equations of motion. Simplifications to the equations beyond those presented in this chapter are possible. These further specifications are due to the fact that some of the terms in the equations are insignificant in most cases. The simplifications due to this fact will be discussed in the following chapter. 60

61 4.2. The Locally Attached Frame The body frame of reference in this case, is attached to a mass particle in the body and its orientation is the orientation of this local area. All of the deformations are measured relative to this local frame of reference. In general, no simplifications occur in the equations. A simple example will show what is involved. Example 1: Two mass problem with locally attached frame The problem to be analyzed is shown in Figure 4.1. It consists of two point masses connected by a spring. The dynamics equations for this system are written below it as if a finite element model had produced them. These equations, as they stand, are good only for small displacements and no rotations. The equations include one rigid body mode and one straining mode. N SU2 m/ /f Figur12. — kk x fim ] (i f: -jt =i ro Figure 4.1 Two mass problem

62 To describe the large displacement dynamics, we attach a local frame of reference to mass. one, and orient the z axis so that it always lies along the spring and passes through mass two. (see Figure 4.2). Now, because the frame is attached to mass one, both 7r and iu are always zero. Letting the distance between the masses in the undeformed state be r, then 72 is simply r5. The vector u2 is simply po. r = 0 r2 -- r Ul 0 iO 2 = -a Now the constants associated with this body can be calculated Total mass (M) M fdm - (ml + m2) Underformed position of C.M. (rm ) U2 a y H 1-r y A x Figure 4.2 Locally attached frame

63. - I rj.; rcm = f J drm = (Oml + rm2) rm2 _ z ml + m2 There is only one deformation coordinate which we are calling s. The deformation vectors u (rt) have only two spatial positions r and r2.,,(,,t) == -. -. = (A u2 r =r2) = Identifying this with -i',,, gives s-{~l Now we can determine 7 and a Ms= Jl dm 2 -= z,, i m,, = (OmI + m2) m 2 -= m22.... 2 2 2 = r t,; X t,, -0 This shows that the deformation does change the position of the center of mass,

64 but does not have an angular component. The Inertia Tensors JO, Jl,, J22, Because the chosen axis are lines of symmetry, the inertia tensor will be diagonal. In fact the moment of inertia around the Z axis is zero because both masses are point masses. The only remaining terms in the undeformed inertia tensor are JOyy and JOz,. JO Y= JOyy JOY z JOyy = J(r2 + rj2) dm ZZ = J(r2 + ry2) dm Since r, and ry are zero for all masses, so these terms reduce to JO = f r2 dm 2 E ri2m, = (02m + r2m) 2 2 mr2 so JO [o m2r2 m r2 which could have been expected. The Jl, tensor has the same structure. There is only one deformation so J1- JIY1y Jyvy = 2f(rz, t + rt',) dm J1zZ = 2J (r2, + ryty) dm Again both r, and ry are always zero so these terms reduce to

65 J = 2f(r O ) dm 2 r,i Or, mi = [(O)(O)m + (r)(l)m21 2m2r so Jl s 2m 2r 2m2 The J2co tensor has again the same structure as ia and and its non zero term is J2 = J2dm - Eb2m2 M (02M, + 12M^) 1M2 The coordinate matrix for the inertia tensor becomes: J = A + A +J2 t2 ==[ m2r2 ]+ 2m2r2 MOr2 4 + m2 2m2r2 2r M2 _ m2(r2+ ru + u2) m2(r2+ ru + u2) The only remaining body constant which needs to be defined is symmetric in i and j, GI, must be zero. GCI. Since G,Y is ant The mass and stiffness matrices were previously defined without the attached local reference frame. With a local frame attached, the mass and stiffness matrices must reflect the fact that the frame defines the rigid body motion of the body. Therefore there cannot be any rigid body modes in the stiffness matrix. These are eliminated very simply in this case. Attaching the frame on mass one ensures that 1i, is always zero. Therefore E61, 1 and 1i are zero also. These facts can be reflected in

66 this mass and stiffness matrices by eliminating the rows and columns associated with 0s1. The mass matrix becomes just m2, and the stiffness matrix is just k. ) [m m2] ( o) - m - -_ The equations of motion are P= -4 -_ H =t where P = MR + f XM(ia+ + sI) + Ms-'l H= * fl + (r,, + i) X MR oo ~ *-. 1 * 1 and m2s + kIM + MR * s - 2 f *(J 1 + 2J2A) ~ nl M = + m2 -~ __ rm2 rcm 2 _ m2 A MI + m2 M 1M A+ iA J =- J0 + J1i + sJ2s 0 JO -- m2r2 m 2r A - 2m2r2 2m2r2 J2 m2 m2

87 4.3. Floating Frames The equations of motion can be simplified by using a floating frame of reference. Determining suitable deformation coordinates appropriate for a floating frame requires more preliminary numerical work than in the locally attached case. But since the equations of motion are solved repeatedly, the floating frames simplicity is a great advantage. There are several types of floating frames, but every one has one thing in common. Each one puts the local reference frame at the instantaneous center of mass. In terms of the constants which we have derived for the local reference frame, rc, = 0, and 8i = 0 for all i. These can be done by first moving the reference frame to the underformed center of mass, and then redefining all of the deformations so that they do not disturb the center of mass. With this, the equations are already much simpler. This completely decouples the position, velocity, and acceleration of the reference frame from angular and deformation coordinates. The equation for linear momentum becomes P =MR This can be directly substituted in the fundamental equation and still maintain a simple result. 6R. MR - f+... The angular momentum also becomes simpler: -_ 44 _- o H- - 'f + (a +d e Gc, )t, The deformation equation is also simpler.

68 00 4 -._ - 0 S1Ik[Mkp p + Kkj, + 0 -(a + m Gt ) + 21 m * Go+ G -2 (nJ' + 2,s2,p1) --- 0 = 0 2 A look at the kinetic energy will give a clue as to what can be done to further simplify the equations of motion. Repeating Equation 3.5.5 here + 4 A 4 _4, 4 T- = MR R + MR fl X M(m + slr)+ Rsj j 1 -* _ G + n J * n + n -(a + u 1 Gm )i,6 1 0 + 2- Mkp,p Putting the reference frame at the instantaneous center of mass has already reduced this to: I o 3 T= 7 MR.R + 2Tk Mkp 7p If it could be done, eliminating the entire quantity in the parenthesis would give the greatest simplification. In general this is not possible by simply redefining the deformation coordinates. What can be done is to eliminate a, for all e. This can be done by a simple redefinition of the deformation coordinates. The deformations are made to be orthogonal to rigid body rotation. The frame of reference associated with this is called the Bucken's Frame. In some cases the term G.m is entirely zero. But in general this is not true. In the following chapter it will be shown that some of the terms coming from GC, are insignificant and can be neglected. Another way to get rid of this term is to introduce a non-holonomic constraint equation to be solved with the equations of motion. It is

69.- e -. t1m Gme r = 0 The problems that this creates, far outweigh any simplification of the equations of motion. For one thing, this constraint is meaningless whenever all l, are zero or all tl, are zero. This causes singularity problems in the system of equations. The other problem is that it doesn't really simplify the equations. The three extra Lagrange multipliers which must be added, reintroduce most of the terms that were eliminated. The frame of reference associated with this constraint is called the Tisserand Frame. In the case where G^, is identically zero the Bucken's and Tisserand Frame are identical. Another simplification which can be seen from kinetic energy is related to the mass matrix. When the deformation coordinates were redefined to eliminate i and a,, a process of orthogonalization was performed. To eliminate 8, and ge the deformations were made orthogonal to the rigid body modes. If they are made orthogonal to each other, the mass matrix becomes diagonal. The kinetic energy becomes 1 -* -* T 2- MR -R 1 - - -# _ + J n flJ n + n * I,, GM' l 2 + -L a + 2 M(k),; tt The final simplification that is possible concerns the potential energy. So far we have assumed a full stiffness matrix. If the deformations chosen would diagonalize the stiffness matrix, the equations would be even simpler. The free modes do just that. They have all of the preceding properties of being orthogonal to the rigid body modes and to each other, plus the stiffness matrix is diagonal. The strain energy becomes U = T K(=))k 7k 2 —

70 To illustrate the simplifications obtained from a floating frame, Example 1 is repeated with a floating local reference frame. Example 2: Two mass problem with a floating frame The problem again is illustrated in Figure 4.1. This time the frame must be located at the instantaneous center of mass. For this to be true the frame must be located at the undeformed center of mass when the deformations are zero, i.e. rcm = 0. See Figure 4.3. This is easily done by subtracting r, from every position vector r. Denoting the new vectors with a prime rt r-,.-- = _ Fcm review: r -0, 2= O = r, r ----- z. The new position vectors are calculated ml + m2 as, I -r' r2 RR., XY Figure 4.3 Floating frame

71 rl r - rce = (O- rm2 -rm2 ml+ m2 r2 r2- rcm =- (r- ) rm I z inl + i2 ml+ m2 A check to see if the new rcm is zero confirms that it is. rcm m - rc- 0 Now the underformed position vectors r have been redefined in terms of the new frame. Next the deformation coordinates can be redefined so that they do not disturb the center of mass. This is done similarly to what was done with the underformed position vectors. There was only one deformation coordinate in example one. The shape function for this coordinate was (2 r rj M The amount that this disturbed the C.M. was 8 = M z. The new shape hi is simply 4 — M -42 1 m2z~ M mX A check to see if the new 7 is zero confirms that it is, 8 = fdm -O —m- O,- +M =0

72 In Example 1 the 'a term was zero already. A check of it with the new deformation shape shows that it is still zero. al= f X ' dm 2 s- a Et x tOil2)mi =0 The cross momentum term G,p is also still zero, because there is only one mode. The inertia tensor can be recalculated. It still has the same nonzero locations. The undeformed inertia tensor is =D [ jY with J =yy JO -=- f rs 1 2dm 2 -rm2 rm TaL 2.M M I mM 1 2 1 m2- M (m1 + m2)m2m1 =mi r MIM1132 2 M The first order inertia tensor is l = Jlyy with

73 Jlyy = Jl 2f(r,', )dm =j2 I,', 1I,m, 1 1 M mM 1 1 = 2 = | ml + j.Mj2] M 'M M -m mim2 = —2 M r The second order inertia tensor is J2 J2, with J2yy = J2 = f( )2dm 2 = ~(d)2m J = = + 'M' m)2 ] m m2 M mIm2 If we let M - the inertia tensor can be written. M i J= + Jlu + uJ2u where JO=[ Mr2 Mr2 iJO = 2Mr 2M1 J2 M M. Now the equivalent mass and stiffness matrices can be calculated. In terms of the original model of the problem, the deformation shape is

74 0 = ^ —m2 M M I The new mass matrix is calculated as O MO M - m[ r 2 1 {M -m2 I _ M m2 ~mlM Mimi [ml -mml 1 M2 "ml "'mlm"2 m2 miOm2j 1 2 2 = -mim2 + m2m1 m 1m2 M M The new stiffness matrix is calculated similarly OT K - 1 L2 } -- M mI -m2 1 1 -m -k(ml + m2)1 MT2 Imil *(ml+ m2)1 -- k In summary, the equations of motion are: MR= f H - J fl -L -- H =t Mi + kt - 1 (J l + 2J 27) nf = a where M = M

75 1 -M2 l ^l -m } [ M(r2+ 2rrn + ) J M(r2 + 2rri + 02 4.4. The Bucken's Frame In this frame of reference the frame will follow the body in such a way as to minimize the square of the local deformations. Because of this it will lie at the center of mass, and no deformation coordinate by itself can introduce a rotation relative to this frame. This can be shown by the following development. Minimize: C = - fI * udm (4.4.1) Taking the variation of C: 6C = - f6u * + u -6udm - f. ii From Equation 3.2.13 6P - SR + dr X q + 6u Since the position of the mass particles are fixed and only the frame and deformations are varying, 6P = 0. Therefore 6R + dr x q + 6u -= 0, or solving for Su gives 6u =- (6R + dir X ) Substituting this into the variations of C and setting the equation equal to zero will give stationarity conditions for the minimization of C. J - (bR + da x j) i udm = 0 Dropping the minus sign, and distributing gives: f6R uI + Jf Xq 'udm O

78 or 6R * Jfdm + Sn X idm =O Since the variation 6R and dir are arbitrary, two constraints result: Jfdm = 0 and X iudm = 0 If the frame is initially located at the center of mass then fdm - 0, and the first equation becomes fqdm = 0 (4.4.3) The second equation (Jfq x udm = 0 can be simplified by expanding q. f7 x -dm = f(r + ) X idm. Since i X u is always zero this becomes.7 Xudm - 0 (4.4.4) The first of these constraints puts the frame at the center of mass. The second is an angular constraint. If the separation of variables is used to express u (u-i, vi ), these become fI(T + i i )dm = 0 and f ( Xqj )dm - 0 In terms of the already defined body constants these are MR,, + Z i,. = 0 (4.4.5) and a, ra 0 (4.4.6) There are a total of six constraint equations here. - The first three fix the frame at the instantaneous center of mass. The second three prevent the deformations from producing a rotation relative to the frame. The simplest way to ensure these are always satisfied is to define the deformation shapes such that - = 0, and a, - 0 for all i. Of course rm can be made zero by moving the frame of reference to the center of mass. If all these terms are zero, the reference frame is a Bucken's frame.

77 4.5. The Tisserand Frame The Tisserand frame minimizes the local kinetic energy. It too, will lie at the center of mass. It's angular constraint however will be non-holonomic, and is closely related to the body's internal angular momentum. The minimization principle is 1 e Minimize: C -. uidm (4.5.1) Taking the velocity variation of C gives 6C = - (6 * u +u * 6 u)dm --- u -udm Taking now the velocity variation of Equation 3.2.13 4 4 ^ -~ _ -L0 0 6P =- R + 6n Xq + iu When set to zero this gives an expression for 6u Su=-(6R + Sn Xq) This can be substituted in the variation of C which in turn is set to zero "Lo 0 J - (6R + 6 Xq) ~ u'dm -0 This results in o o 6R * J idm + 5 *J q Xu — = 0 Since the velocity variations are arbitrary the following constraints result. fudm = 0 Jq x udm 0 The first constraint can be integrated to fJ( + +i)dm = 0 or Jfdm = 0

78 The second term cannot be integrated to holonomic form and is therefore nonholonomic. In terms of deformation coordinates (u = i,, ), these become (from Equations 3.3.3 and 3.4.10): MAc, + i-i, = 0 (4.5.2) - O0 and (a, + q7 G^,, ), 0 The first term again establishes the frame at the center of mass. The second term can be recognized from the angular momentum equations, as the angular momentum of the deformation relative to the frame of reference. This constraint forces this to be zero. The angular velocity of the frame of reference, in this case, will completely specify the angular momentum. The deformation can be redefined to zero si and a,. A simple redefinition however cannot zero Gi,. In some cases it will be zero because of the nature of the body, but this is not generally true. If it is zero the Bucken's and Tisserand frames are identical. In order to use this frame when it is not zero, the constraint,,, G,, r, = 0 must be included with the equations of motion. This causes a problem with the current methods of solution. When either all a, are zero or all ia, are zero the constraint is automatically satisfied. In this case the constraint "turns off" and produces three zero columns in the matrix being inverted. For this reason the Tisserand constraint was abandoned. Another reason for abandoning this constraint, was the it did not really simplify the equations overall. The dynamics terms were simpler, but including the constraint reintroduced most of the eliminated terms anyway.

79 4.8. The Free Modes as Deformation Coordinates As was stated earlier, the free modes of deformation give the greatest simplification to the equations of motion. The reasons for this are the "nice" properties of the free modes. First of all, the free modes are the natural modes of vibration for a body freely floating in space. This means it is not attached in anyway to anything else. For general mechanisms, of course, this is not the case at all. Each body is attached in some way to outside influences, and it does not vibrate with its free modes. This does not stop us from using these modes as deformation coordinates, however. These additional boundary conditions are enforced by adding constraint equations corresponding to each boundary condition. This introduces a Lagrange multiplier for each constraint equation. If the deformation shapes of the body with the new boundary conditions can be spanned by the free mode shapes, the equations will be correct and accurate. If, however, the new deformation shapes are not easily spanned by the free mode shapes, a lot of free modes may need to be carried. In this case, the techniques of the following section can be used to generate shapes with most of the properties of free modes, and at the same time span the space with fewer modes. The first important property of the free modes is that all of the modes are orthogonal to each other with respect to the mass of the body. We will presume that the modes are scaled such that fi, * dm - (4.6.1) In matrix notation this states that the mass matrix for the free modes is the identity matrix. Since the body is freely floating in space, six of the modes are rigid body modes which introduce no strain in the body. The first three are translations. For each of

80 these translation modes every particle translates equally in one direction. These modes can be written as =a -eb) = 1,2,3 where e b) is a unit vector for some frame on the body. Substitution of this into the orthogonality property shows that every other mode is orthogonal to the rigid body translation modes. J 'b), dm = 0 a 1,2,3 j > 3 The unit vectors are constant over the body and this becomes,a *J dm =0 a = 1,2,3, j>3 Since the unit vectors are arbitrary the following is true f dm - 0 j>3 (4.6.2) In terms of previously defined body constants this is -= 0 for j > 3. This indicates that no free mode (other than the rigid body translation modes) disturbs the position of the center of mass. The same can be done for rigid body rotations. In this case, a rigid body rotation needs to be described. In a rigid rotation, a point moves perpendicular to the vector representing the rotation. The motion is also perpendicular to its position vector from the origin of the rotation. The magnitude is proportional to the length of its position vector. This can be stated concisely for rotations around unit vectors C(b) as A = w- ed(X r a 1,2,3 Substituting this in the orthogonality condition yields feb) Xr.* dm = 0 a 1,2,3 j >6 Exchanging dot and cross and moving the unit vector outside the integral gives eb^ f Xq dm = 0 j>6

81 Because the unit vectors are arbitrary, the integral must be zero and we have J X j dm = 0 j > 6 (4.6.3) This is just the body constant aj. This shows that for all the modes which cause straining in the body i = 0, and aj = 7. This gives the desired simplifications which had been hoped for. Since these were also the conditions for a Bucken's frame, the local frame of reference is a Bucken's frame. One further property of free modes gives another simplification. The modes are also orthogonal with respect to the body stiffness. 4-+ f(q7 ):E: ( )dV= W(2 6, The stiffness matrix is diagonal and the potential energy becomes: 1 U I- w. n( o, (4.6.4) With these properties of the free modes, the large displacement equations for the linearly flexible body become 6R MR- f? lH-t ) (4.6.5) + 6tr (k + W(k)lk - 2 a (J- k + 2J2kp i)'f A4 _ - e 1 + f A', Gmk + 2fl G #ma -a -rj 0 where H - n + 4.. G. t/0 H = J - fl + 17m G., t17

82 4.7. Simplifications Using Gram-Schmidt Orthogonalisation As was stated in Section 4.6 the free modes do not always easily represent the deformations actually occurring in the body. In this case, very many modes would have to be carried to get an accurate representation. The major advantage of the free modes was their orthogonality to the rigid body modes and to each other. But they are not the only shapes which have this property. What is desired then, is a set of shapes which easily represent the deformations of the body with its new boundary conditions, and are orthogonal to both rigid body modes and each other. This can be done by starting with a set of shapes which more closely satisfy the boundary conditions. These shapes can be static deflection shapes, mode shapes or other assumed mode shapes. They can be analytically derived, derived from finite element models, derived from engineering judgement, or derived from modal test data. They will not, in general, be orthogonal to the rigid body motion of the body. If they are the natural modes of the body with the boundary conditions, they will be orthogonal to each other. This is not necessary though. The deformation shapes are made orthogonal to the rigid body modes by subtracting from each shape any rigid body component. This orthogonalizes each mode to the rigid body modes. Next each mode can be orthogonalized to every other. When this process is done, every deformation shape is orthogonal to both the rigid body modes and every other mode. The resulting shapes must still span the same space as the original shapes. This process is called Gram-Schmidt Orthogonalization. After the new shapes have been defined, the new body constants are calculated Since the shapes are orthogonal to the rigid body modes both bs, = 0 and i, = 0 for

83 all i. The mass matrix will be diagonal, but the stiffness matrix will not be (unless they turn out to be free modes).

CHAPTER V SIMPLIFICATIONS DUE TO DROPPING NEGLIGIBLE TERMS 5.1. Introduction In the previous chapter we saw that the equations of motion for a flexible body could be simplified considerably by a suitable choice of reference frame and deformation coordinates. In this chapter further simplifications are explored by examining the relative contributions of each term in the equations. The terms will separate into levels according to their relative contribution for each equation. 5.2. Sizing Assumptions In order to attach relative sizes to the terms in the equations some assumptions need to be made. These assumptions may not always be accurate, but they do give an indication of the relative importance of various terms. The equations to be used in this analysis are the equations for the free modes because they are the simplest. The first assumption that is made is that all the terms from the rigid body equations are important and will be kept. For example MR and JO - n are both found in the rigid body equations and are therefore kept. The second assumption is that the solution to the it modal deformation equation is sinusoidal in character. We will let the magnitude of the solution (It, 1) be e, and the 84

85 frequency of oscillation be w0. Taking the time derivative of the solution gives the magnitude for the rate of deformation (Y, ) as w. e. Another derivative gives the magni00 tude of u, as w2e. If the members are sufficiently stiff we can assume that the solution frequency is less than the mode's natural frequency (w, < wi ). Thus we make estimates of the magnitudes of u, and r, by using w, instead of wo: 1 Ii, 1 - 7 { - wa, e (over estimated) 1l7, w = e,2 (over estimated) In the case of transient vibration, these magnitudes are good estimates instead of over estimates. 99 9 The modal equation in a non-rotating frame is: v, + 2,)w(,)v, + w(2li = ai. Writing the magnitude estimates below each term results in: it + 2jp(,), + W(vqs. - as we 23)w(iE^ W(2)E le I Since all these terms may be important they are assigned a value of approximately one 1 for convenience. This gives the relationship: w,2e = 1 or w, = c 2 From this, the e oe magnitudes of u7, and 7r, can be expressed in terms of e. k l-,1 0e Now these expressions can be used in the equations for a flexible body to see the relative size of the terms. The flexible body equations will be expressed with the derivative of the angular momentum explicitly carried out. This is needed because the

86 derivative of a quantity is of different magnitude than the original. For accurate comparisons then,the derivatives are completely carried out. _R( MR- 1 +.(7 n +n x7 4n+ n (5.2.2) 1 1 1 -4 — + O 0 + n x~ lm Gm# c, + ilm Gm<7 l, - t 2 1 +&l, W(. + w(. + fI - i GO, + 2 * Gm *~ 1 1 + 62 2n (ni. + J2 + 2, ) n-, -= o 1 6 1 J ' J+ t rit + J2e,,ris At 1 e e2 J ti Ak + 2J2,,,, e 1 3 2 2 Written below each term in each equation is its corresponding size according to Equations 5.21. The center of mass translation equation is entirely unaffected by the modes, so it is entirely on the order of one. The frame rotation equation does have different orders of magnitude, as well as the modal equations. These different orders of magnitude will be assigned levels. Level zero are those terms on the order of rigid

87 body equations and straight modal equations. These are the terms with one under them (e -- 1). The other levels again correspond to the power of epsilon. These lev1 3 els are i, 1,, and 2. Table 5.1 contains the terms added at each level for the 2 2' frame rotation equation. Table 5.2 contains the terms added at each level for the modal equations. LEVEL SIZE 0 1 1 1 TERMS JO n + n x Jo n - t 1 e1 3 Jlt l, ~fl AJ, 'i * n + nf x As,, ' f + tmi Gme, te —.. 0 -. m.. r,., ~4 - + -* N J2m. 7,, y. + f+l X IN J2, -f I 3 2 2 Table 6.1 Size of terms for frame rotation LEVEL 0 SIZE 1 TERMS h71 + 23,w(9(,l + W(Iq 1 - - 1 2 1 1 f2?-+ -4 2 n G.,,,,,, f1 A4 -4 1 -* n * CGh 1 - 2 — n 2J2,,?, 2n 2 Table 5.2 Size of terms for modal equations

88 Level zero in the angular rotation equation corresponds exactly with the rigid body equations for angular rotation. At level one-half a term involving the velocities of the modal deformation is added. At this level the deformations can directly influence the rate of change of angular momentum. The next level down (level one), part of the change in the inertia tensor from the deformations is included along with a cross momentum term. Another level down includes another cross momentum and second order velocity change in the inertia tensor. The last level brings in the second order change in the inertia tensor. In some techniques for analyzing flexible body dynamics, the rigid body problems is solved first, independent of any deformations. The forces and motion derived from this analysis are then applied to the flexible bodies individually to determine the deformation. Level zero corresponds, in part, with this analysis. The expression for the angular rotation equation is entirely independent of the deformations, as it is in the "separated" technique. However, this is misleading because even at level zero the angular rotation equation does depend on deformations. The dependence comes through the torque term. The torque term includes torques due to reaction forces at the joints. These same reaction forces are also affecting the deformation coordinates which are in turn affecting them. All of these effects are happening simultaneously through the action of the Lagrange multipliers (explained in more detail in Chapter VI). Although the rotation equation is not explicitly dependent on the deformation, the reaction torques are, and this is already a better approximation to the actual dynamics than the "separated" technique. For the modal equation, level zero corresponds to a straight modal analysis problem, with the inclusion of the "centrifugal forces" in the analysis. This covers most of

89 k k z Figure 5.1 Planar three mass problem the ground needed for an accurate representation of the deformation. Level one half brings in the Coriolis effects due to the rotating frame. Level one includes both an Euler acceleration term, and a second order centrifugal force term. As it is expected, the Coriolis term is as large as the modal terms and their centrifugal term. Again as expected, the Euler acceleration term is even smaller. To examine the validity of this sizing procedure a simple example was run, in which these various levels were implemented. The problem was run at each level and the results compared. The example was a planar problem, but each term in every level existed so that going to a new level always changed the equations. The example was a single body made up of three point masses connected by three linear springs. See Figure 5.1. The mode shapes for this body can be seen in Figure 5.2. The body constants such as total mass, inertia tensor, cross momentum, etc.. are given in Table 5.3 The

90 X v \ / / /- -A Figure 5.2 Modal shapes for three mass problem three mass body was given an initial angular velocity of 10 radians/sec and rl — r2 -2.=111, 3s - 0. Because this body is not connected to any others, only terms in the equations effect the solution, and there are no applied forces and torques to complicate the interpretation. Examination of the plots for omega showed a big change in going from level zero equations to level one-half equations. (Some of these plots are given in Chapter VII: Figures 7.5 - 7.8). For level zero equations, omega is constant. For level one-half, the

91 angular velocity starts to oscillate because of expanding and contracting of the body from oscillations in vs. Going to level one changes this oscillation slightly. Higher levels show barely any change. _1st mode 2nd mode 3rd mode I-Ix1 0-l y 0.5i - 0.866y 0.866B + 0.5y 0.5i+0.886y -0.866z+0.5g ox + iV o0.866 - 0.5y -0.866-0.5 J 3.0 2Mr I w A.Z m 1.5 0 1.5 0 0 0 0 O} MJ 0 -M 0 M 0 0 O O 0 O} 0 M = 3m Jo = Mr2 Table 5.3 Body constants for three mass problem

92 In percentages these changes are: level 0 to level 2 25% level -to level 1 3% 2 level 1 to level 3 not noticeable Examination of the output for the modal deformations shows similar results. The shape of the solution changes dramatically from level zero to level one-half for both vi and t2. In going from level one-half to level one there was no perceptible change in any of the deformations, or for any further levels. From this it seems that a reasonable choice is to include level one-half in the dynamics equations and leave out the remaining terms (shown in Table 5.4). This gives a simple set of equations to program with little numerical work and an accurate approximation to the equations of motion. In most cases the deformation equations are approximate to some extent anyway. With this choice it is unnecessary to calculate the second order inertia tensor 12, which had 9 N2 terms in it. The crossmomentum term still appears in the modal equation, so it still needs to be calculated. Canavin and Likins in the equations that they propose completely eliminate the cross-momentum term [18j. This was done with an assumption at the very beginning -. 9 that tlm Gm,,s was second order in deformation and could be eliminated from the Tisserand constraint equation. This corresponds to eliminating this term from the kinetic energy. If all second order terms in deformation had been eliminated from the kinetic energy, there would be no deformation equations at all. Here we see that this term is usually small, but should be included in the deformation equations because it does affect the solution. Notice that with this choice of levels, the highest derivative

93 terms are again completely uncoupled. This will again give a diagonal and constant matrix for the explicit integrator. 6R.MR - f + si ri i + 2q,)ow(,), + w(,9 -n -1ni 1f + 2n * GM., - a- 0 Table 5.4 Level one-half equations for a flexible body

CHAPTER VI THE MECHANISM AS A WHOLE 6.1. Assembling Bodies, Joints and Forces Up to this point we have only examined the equations of motion for a single body. This of course does not cover the field of mechanisms, where there are many bodies in each mechanism interacting with each other. In order to cover this territory, the equations for the various kinds of interactions will also be derived. The interactions between bodies can be separated into two major categories. The first category contains all of the different types of constraints. A constraint is a connection between two bodies which prescribe a particular relative motion between the bodies. The first kind of constraints are joints. For example a hinge (called a revolute joint) is a connection between two bodies which only allows relative rotation around the hinge axis, and no relative translation of the two bodies at the hinge point. A ball joint (called a spherical joint) also prevents two bodies from translating relative to each other at the joint, but allows all three relative rotations around the joint. If the relative motion prescribed between two bodies is time dependent, this is no longer called a joint, but instead is called a generator. It is called a generator because it is generating some relative motion as a function of time. An example of this is a 94

95 four bar mechanism, in which the crank angle is a prescribed function of time. The constraint that prescribes this motion as a function of time is called a rotational generator. The second category of interactions between bodies contains the applied forces. The term applied forces means forces that are applied to, or between, bodies which are an explicit function of displacements, and/or velocities, and/or time. These are not to be confused with reaction forces. Reaction forces are forces on bodies due to the action of some kind of constraint. They are not explicit functions of displacements, velocities, and/or time, but are the forces at the constraints which are necessary to impose the constraint. An example of an applied force is a spring. The force that a spring applies between two bodies is a function of its length. A shock absorber is another example of an applied force. In this case, the force applied is a function of the velocity of extension or compression. Once the equations for each particular interaction have been derived, they can be used over and over again. The equations become a kind of library from which a particular interaction can be called whenever needed. In the case of constraints, the underlying variations which were previously assumed independent, have become dependent. There are several techniques for solving for a set of independent variations. The technique used here is suited to automatically generating the equations of motion. Lagrange multipliers are used to enforce the constraints. In this case, as many multipliers as there are constraints are added to the set of coordinates. Each multiplier is scaled by the amount that its constraint affects each equation and then is added to that equation. This is shown in the second line of the fundamental equation. To make the whole equation set solvable, the constraint

98 equation is appended to the equations and solved simultaneously. This procedure is easy to automate. But it results in a large set of equations. Happily, this set of equations is fairly easy to solve, even though it is larger. This stems from the fact that the equations are loosely coupled. This causes the matrices involved to be sparse (have mostly zero terms). Subroutines exist which can decompose a sparse matrix in far less time than full matrices. In this way, both ease of automation and speed of solution are simultaneously obtained. Applied forces, however, do not make the variations dependent. All that is needed is to add to each equation the component of the force in the direction of the equation's underlying variation. These components are called the generalized forces and are shown on the third line of the fundamental equation. Again this is an easily automated procedure, and computationally efficient. It is not necessary to include extra equations in this case, but it is usually done for the implicit formulation to improve sparsity. For the implicit formulation, the partials of the generalized forces would need to be taken. Since the expression for the force is usually complicated, this can create a great many complicated terms. New variables are defined to remedy this situation. In terms of these new variables the equations become much simpler and consequently the partials become simpler also. The equations defining these new variables are appended to the equations of motion and solved simultaneously with the rest of the system. This is not needed in the case of the explicit integrator, because the applied forces are not functions of accelerations or constraint forces. Therefore the generalized forces are evaluated on the right hand side of the equation set and do not enter into the matrix multiplying the highest derivatives.

97 F., AR Figure 6.1 Spherical joint In both the explicit and implicit forms for the equations, a matrix decomposition must be performed for general mechanism problems. In the explicit case, a matrix decomposition is needed to solve for the highest derivatives. In the implicit case, a matrix decomposition of the Jacobian matrix is needed. Particular attention has been paid to put the equations in a form which simplifies both of these matrices.

98 8.2. Constraint Equations The constraint equations for various types of joints are given in this section. There are three fundamental constraint equations from which many different joints can be formulated. The Spherical Joint Constraint The first constraint corresponds directly with the spherical joint. This vector constraint is illustrated in Figure 6.1. The constraint equation states that the sum of the vectors around the loop is zero. $- AR+ A _(BR BQ) = 0 (6.2.1) This constraint will be imposed by a Lagrange multiplier vector X. The contribution of X to each equation is - X. This is as shown in the second line of the fundamen9q tal equation. The constraint must also be included in the equations of motion as shown at the bottom of the fundamental equation. The form that these equations take is different for the explicit (non-stiff) integrator, and the implicit (stiff) integrator. The derivation for both these forms is included as Appendix B. For the explicit integrator the equations for body A, and body B, and the constraint equations for the spherical joint take the following form:

99 M(a)l Jo (') 1 M(b) _1) 0 A qf af R(a) n(a) fl(a) R(b) (b) ) (6.2.2) JO (b) I (q()A a' )T ((a) A a ) #10 -. The constraint adds a column corresponding to the associated Lagrange multiplier. It also adds the bottom row, which corresponds to the constraint equation in second order form. This matrix has a very nice form. It can be partitioned in the following way: m r I = i where M is a diagonal matrix of constant terms and r is a matrix of time varying terms. This structure can be taken advantage of in solving for Q and X.

100 The implicit equations are similar. The Jacobian matrix in this case is: (6.2.3) q(-,(B x( ()A -1M(& X)B +X()+ l) 7 )A bf (X_)_'))T _ -B1)T A_,' A ~f |I ($()AI)Tbi (.)-T A)!rq 1 (-q(b)A)TB ')r s)A r t )r The Jacobian matrix has more terms in it than the matrix used for the explicit case. The terms populating the center come from taking the partials of the terms in the equations which were added by the constraint. They involve second partials of the constraint equation. The matrix is no longer symmetric as it was for the explicit case. The partials with respect to angular coordinates produce terms which are nearly symmetric, except the terms are post multiplied by the B matrix. The Rotational Constraint The universal joint is similar to the spherical joint, but constrains one more degree of freedom. The degree of freedom constrained is rotational in nature. The equation of constraint is best understood from Figure 6.2. The universal joint contains a cross bar with two perpendicular bars. No matter how body A and body B move, the cross pieces will always be perpendicular. The unit vector a moves with body A

101 b BODYA / BODY B Figure 6.2 Universal joint and represents its cross bar. Unit vector b moves with body B and represents the other cross bar. The fact that the two remain perpendicular can be stated simply as a * b 0 (8.2.4) In the case of rigid bodies both a and b are fixed relative to the local frame. These unit vectors can move relative to their local frames if the bodies are deformable. In the case of deformable bodies this small change in orientation will be represented by the vector 9. Since the vector 9 is small the unit vector can be represented as the undeformed unit vector a' plus the change due to 6. Concisely stated this is a = a0 + 8 xa0. The variation of a measured in the fixed frame is () = (a) X a + (a). Because a, is fixed in the local frame (Sa, =0. The variation of a becomes =a = a) Xa+ x a and similarly for b.

102 The vector () is the sum of all angular changes of the modes. Therefore = (a)t, where 0, is the participation of mode i in an angular orientation change of this point. The variation of a in terms of frame rotation and deformation is 6a ^=(a) Xa + 5 a () X ao (6.2.2) These properties can be used to derive the partials of the constraint equations., -('P( - i '*) - * - o Figure 6.3 Translation constraint

103 This will give the contribution to each equation in each body's equation set. The Translational Constraint Another remaining type of constraint used in joint formulation is translational in nature. This constraint constrains the translation of points on different bodies in a certain direction. Figure 6.3 illustrates this situation. The point b on body B is constrained to lie on the surface defined by the unit normal a on body A. The constraint equation is =- (p()- p()). 0 (6.25) If body A deforms, a it will not correctly specify the plane of contact. Therefore this constraint should be used when body A is assumed rigid, but body B may be flexible. 8.3. Joints The joints are generated by combining the one or more of the constraint equations until the characteristics of the joint have been fully described. Many different joints can be described unsing these fundamental constraint equations. Figure 6.4 shows some of these possibilities. Table 8.1 gives the characteristics of these joints and the constraint equations used to formulate them.

104 SPHERICAL UNIVERSAL REVOLUTE I I I I I 1 I I I I I SLOT CYLINDRICAL TRANSLATIONAL PLANAR POINT CONTACT PLANAR Figure 8.4 Joint types

105 Name Degrees of Constrained Degrees Uses Freedom of Freedom Equations Spherical 3 3 1 Universal 2 4 1,2 Revolute 1 5 1,2,2 Slot 4 2 3,3, Cylindrical 2 4 3,3,2,2 Translational 1 5 3,3,2,2,2 Planar Point Contact 5 1 3 Planar 3 3 3,2,2 1 - Spherical Joint Constraint (6.2.1) 2 - Rotational Constraint (6.2.4) 3- Translational Constraint (6.2.5) Table 8.1 Constraint equations for joints 6.4. Applied Forces In addition to joints influencing the motion of bodies, there are also applied forces. These are forces applied on, or between, bodies which are a function of displacement, and/or velocity, and/or time. They are not explicit functions of acceleration or reaction forces. The simplest and most common forcing effect is gravity. Gravity It is easiest to include gravity in the formulation by including it in the calculation for potential energy. Since the frame is located at the center of mass the potential energy due to gravity can be written simply as U, = - M * R where M is the mass of the body, g is the acceleration of gravity vector, and R is the position of the center of mass. The variation is done to obtain the partials giving: 6U, - Mg 6R The contribution to the equations occurs solely in theequations for the position of the center of mass, and is - Mg. The center of mass equations become

108 14 MR Mg=0 O or MR =- M In matrix notation this is MR - Mg The Line Force This applied force acts between bodies along a straight line connecting them. This is illustrated in Figure 6.5. The forcing element itself is assumed to be massless. The force by convention, is positive it if acts to separate the bodies. P- -pB Body B Pb Body A pPe Figure 8.5 Line force

107 In order to calculate the direction and magnitude of the force new quantities will be defined. AP-P' - Pb l vaP -AP P is the spanning vector between the two bodies along which the foe acts. The disAP is the spanning vector between the two bodies along which the force acts. The distance between bodies (positive) is I, and v is the rate of change of the length. The magnitude of the force itself is a function of I, v, t and possibly other quantities. f = f(l,v,t) Common types of forcing functions are linear springs, dampers, and constant forces. All of these can be combined in the simple forcing definition f = K(I - l ) + CV + f. The force vector is just the magnitude times the unit vector along the spanning vector C4a f l f The contribution to the equations of motion for this force is: F. These are the 9q generalized forces as seen in the fundamental equation in Section 1.6 These equations can be used to formulate the contribution of the force in either explicit or implicit form. In the explicit case the equations are used to evaluate the force and its contribution to each equation in the equation set. For the implicit case the equations are actually included in the equation set so that the partial derivatives are simpler to formulate.

CHAPTER VII NUMERICAL RESULTS AND CONCLUSIONS 7.1. Introduction The goal of this chapter is to confirm the validity of the formulation for the flexible mechanism problem presented in the previous chapters. This will be done by reviewing some of the results of testing the formulation at various stages. Most of the examples which were used as test cases were simple enough to either: (1) predict the solution analytically, (2) the effects which were being tested could be anticipated, (3) the results could be compared to the output of other existing programs. Simplicity was preferred in testing so that effects being validated could be isolated, tested and understood without loosing them in other effects. The following areas were considered in validation of the formulation. 1. Rigid body problems 2. Free vibration 3. Vibration with boundary conditions 4. Coupling of flexibility to mechanism behavior 5. Cross-coupling effects 6. The non-stiff explicit formulation 7. The stiff second order implicit formulation 8. Planar Problems 9. Spatial Problems 108

109 7.2. Testing and Validation Most of the initial testing was done for planar problems. This was because it was simpler and easier to program. The formulation that was used was the non-stiff explicit type of formulation. As the testing proceeded, the program that evolved (DYMOS1) became fairly generalized. It was an experimental program and was used for testing correctness, not efficiency. It did not use sparse matrix techniques, but instead used full matrix decomposition routines available on the operating system. The first area tested for correctness was how well the formulation gave the rigid body response of mechanisms. To validate this area, the response was compare to the output of DRAM2 At least six different rigid body problems were solved: 1. Single pendulum under the action of gravity 2. Single pendulum with an impulse applied at the tip 3. Single pendulum with a rotational generator imposing a rotational motion about the pin 4. Double pendulum with gravity 5. Three body robot mechanism 6. Four bar mechanism. In all of these runs the results were indistinguishable from those produced by DRAM. This gave a high level of confidence for the validity of the terms in the formulation which represented the rigid body motion. The next area of testing validated the free response of a single flexible body with no interactions with other bodies. Several modes of vibration were given a flexible body at rest. Since in this case there is no interaction between the rigid body motion and the deformation, the response should be simple oscillation at the frequencies of the 'The program's name DYMOS stands for DYnamic MOdal Simulation. Inquires about it should be directed to the person responsible for coding it: Brian Dent, ADAPCO, 60 Broadhollow, Melville, New York, 11747. Phone: (516)549-2300. 2DRAM is an acronym for Dynamic Response of Articulated Machinery. This planar rigid body program is marketed by Mechanical Dynamics Incorporated, 555 S. Forest, Ann Arbor, MI 48104.

110 free-modes. This indeed was the case. with the frequencies and damping ratios pie output. The period and rate of decay corresponded which were given. See Figure 7.1 for an exam The next step in validation was to test the structural vibration response of a member with boundary conditions. The test consisted of applying a constant load of one pound to beams with different boundary conditions. Three free-modes of a 36 inch long beam, with a cross section that was 0.25 inch square were used for each case. The first three free-mode natural frequencies for this beam were 40.3, 111.2 and 217.9 Hz. Two things were considered in each case. First, the resulting static deflection was r ME 0.0 I.OOOOOE-02 2.OOOOOE-02 3.OOOOOE-02 4.OOOOE -02 5. 0000E-02 6.OOOE-02 7.OOOOOE-02 6.OOOOOE-02 9.OOOOOE-02.OOOOOE -01 1.10000E-01. 20000E-01. 3000OE-01. 40000E-01 1.50000E-01. 60000E-01 1.70000E-01 t. 80000E-0. 90000E-01 2.OOOOOE-01 2. 1OOOOE-1 2.200000E-1 2.30000E-01 2.40000O-01 2.50000E-01 2.60000E-01 2.70000E-01 2.80000E-Ot 2.90000E-01 3.OOOOOE -0 3.0000E-01 3.20000E-01 3.30000E-01 3.40000E-01 3.50000E-01 3.60000E-01 3.70000E-01 3.80000E-01 3.90000E-01 4.OOOOOE-0 4.10000E-01 4. 20000E-01 4.30000E-01 4.40000E-O0 4.500OE -0 4.60000E-O1 4.70000 -01 4.80000E-01 4.90000E-01 5.000OOE-01 ETAI EFA2 1.00000E00 6.0000-0 9.91950E-01 5.80790E-01 9.68872E-01 5.26595E-01 9.31271E-01 4.410786-01 8.80089E-01 3.30285~-01 8. 16459E-01 2.01700E-01 7.41686E-01 6 37968E-02 6.57214E-01 -7.46234E-02 5.64597E-01 -2.04978E-01 4 65461E-01 -3.9382E-O0 3 61572E-01 -4.11122E-01 2.54608E-01 -4.75073E-01.46333E-01 -5.07951E-01 3.84653E-02 -5.08490E-01 -6.73237E-02 -4.77457E-01 -.69437E-01 -4. 1756tE-0 -2.66375E-01 -3.33236E-01 -3. 56755E-01 -2.30323E-01 -4.39331E-01 -.15667E-01 -5.13011E-01 3.31938E-03 -5.76668E-01 1.t9152E-O0 -6.30149E-01 2.24712E-01 -6.72288E-01 3.13701E-01 -7.02904E~01 3.80993E-01 -7.21806E-01 4.22976E-01 -7.28992E-01 4.37688E-01 -7.246396-01 4 24885E-01 -7.09102E-01 3.86034E-01 -6.82901E-01 3.24192E-01 -6.46709E-01 2.43776E-01 -6.01337E-01 1.50239E-01 -5.47722E-01 4.97112E-02 -4.86902E-01 -5. 14003E-02 -4.20006E-01 -1.46808E-01 -3.48225E-01 -2.30733E-01 -2.72800E-01 -2.98255E-01 -1 94996E-01 -3.45592E-01 -1 16086E-0 -3.70311E-01 -3.13255E-02 -3.71440E-01 4.00610E-02 -3.49497E-01 1.14902E-01 -3.06418E-01 1.86094E-01 -2.45400E-01 2.52619E-01 -1.70676E-01 3.13556E-01 -6.72231E-02 3.68093E-01 -4.357211E-04 4.15541E-01 8.42188E-02 4.55335E-01 1.61542E-01 4.87046E-01 2.26918E-01 5.10380E-01 2 76590E-01 5.25816E-01 3.07865E-01 5.31431E-01 3.19263E-01 ETA3 3.OOOOOE-01 2.62383E-01 1.63152E-01 2.69495E-02 - 1.12305E-01 -2.20413E-01 -2.71372E-01 -2.54103E-01 -. 74690E-01 -5 49420E-02 7.51242E-02 1.83189E-01 2.43222E-01 2.41898E-01 1.81199E-01 7.76841E-02 -4.20452E-02 - 148095E-01 -2. 14876-01 -2.27006E-01 - 183069E-01 -9.54012E-02 1 31749E-02 t 15137E-01 1.85627E-01 2.08535E-01 1.79822E-01 1.08053E-01 1. 19533E-02 -.41340E-02 - 156467E-01 - 187909E-01 - I.72040E-01 - 114359E-01 -3.04664E-02 5.81424E-02 1.29545E-01 1.66749E-01 1.61727E-01 1.17093E-01 4.51241E-02 -3.55799E-02 - 1.048781-01 -1.46100E-01 -1 50031E-01 -.16932E-01 -5.61619E-02 1.642386-02 8.25863E-02 I.26267E-01 1. 37463-01....4.... _........4.-.......4....4........4+...3 1 3 21 1 3 2 1 1 3 2 2 I S 3 2 1 t 3 2 1 3 2 a t 1 3 2 1 I 1 23 2 t 1 2 3 1 4 12 1 3 1 2 1 3 t 2 21 3 2 32 1 12 3 3 + 2 1 3 + I 2 1 3 1 1 3 2 I 1 31 2 I 1 3 2 1 3 2 t 1 3 2! t 3 2 3 2 lI 3 32 I I 3 2 l fI 3 2 1 1 3 2 t 1 1 3 2 1 13 3 2 + 1 3 2 2 1 3 2 t 2 1 3 1 2 2 3 1 I 2 1 3 + 1 2 3 1 2 2 1 3 ~ I 3 1 1 312 3 32 2 1 3 2 3 2 1 t 3 0 1 3 2 1 I 3 2 1 1 + 3 2 1 I 3 21 2 3 2 t 1 312 2 *.1, C -* -—.** - - -* ----+ - 32 —4 —4.-+.4 ~= -. a = o.5, =.03, = 12.57, w2 = 25.13, a = 50.27 Figure 7.1 Free vibration

111 compared to the known theoretical value. Secondly, the fundamental frequency of oscillation in the response was compared to the theoretically correct value. The different test cases are shown in Figure 7.2. The results from DYMOS are compared to the theoretically correct value in Table 7.1. Deflection (in.) Frequency (Hz) Theoretical Calculated Error Theoretical Calculated Error Case 1.0994 0.0992 0.2% 17.46 17.77 0.5% Simply Supported Case 2 1.59 1.1.45 28.0% 7.8125 6.34 19.0% Cantilever Case 3 0.0435 0.0277 36.3 % 35.71 27.73 22.3% Fixed-Hinged Case 4 0.0248 0.0136 45.2% 62.5 40.3 35.5% Fixed-Fixed Table 7.1 Boundary condition test results First of all, for every case, the boundary conditions imposed by the Lagrange multipliers have pushed both the deflection and the fundamental frequency in the right direction. For the simply supported case, both the deflection and frequency are correct. In fact, the first of the free-modes accounts for almost all of the deflection and frequency determination. Removal of the second and third modes changed the calculated values less then one percent. An examination of the mode shapes for a simply supported beam shows that they are very similar to the free modes in shape. The cantilever beam, however, does not show such good results. In this case the deflection is nearly 30%o too low and the fundamental frequency is nearly 20% too low. Calculations by hand showed that the first mode contributed 53% of the total deflection, the second mode contributed an additional 11%, and the third mode contributed

112 Case) ---- _ _ Case 2) Cue 3) _Case 4) Figure 7.2 Boundary condition test cases an additional 10%. Carrying another mode would have contributed another 4.6%, a fifth another 4.2%, a sixth an additional 2.5%. Altogether, six modes would still account for only 86.4% of the total deflection. These poor results are related to how well the free modes approximate the shape that the deformation needs to take. An examination of the cantilever modes shows that they are quite different in shape than the free modes. In order to get accurate results for this case the Gram-Schmidt process should be used. Starting with either the cantilever mode shapes, or possibly the static deflection shape, an appropriate set of modes orthogonal to the rigid body modes could be constructed.

113 Even poorer correlation is seen for cases three and four. This is again related to the disparity of shapes actually used and those natural to the boundary conditions. Again this could be rectified by using the mode shapes natural to the boundary condition and orthogonalizing them to the rigid body modes. In order to demonstrate that the formulation was indeed coupling the flexibility to the overall mechanism behavior, DYMOS was altered to behave as superposition approach does. In this uncoupled mode, the program solved the rigid body problem independently of the deformations and then applied the resulting forces to the flexible members. Three different problems were run and compared with the results from DYMOS in coupled mode. One observation from this experiment was that uncoupled bodies vibrated at higher frequencies than the bodies in the coupled run. A second observation was that the uncoupled runs allowed the bodies to separate at the joints because of deformation. A third observation was that the reaction forces were independent of deformation. And finally, rigid body motion was unaffected by the flexibility of the members. All of these undesirable characteristics were not present in the runs in which the formulation was expressed in coupled form. One further experiment was done in the area of coupling flexibility to mechanism behavior. This involved the simulation of a slider crank mechanism derived analytically using one transverse mode and one longitudinal mode by Chu and Pan [36]. It has served as the basis of several comparisons of flexible formulations [7]. The output given by Reference [7] is shown in Figure 7.3. The output produced by DYMOS is shown in Figure 7.4. The flexible beam used in DYMOS was represented by the first three free modes of vibration. From this we see very good agreement in both magnitude and shape of the resulting deflection.

114 0 1 2 3 4 5 6 7 8 at Figure 7.3 Center point deflection of connecting rod for slider-crank (From Reference 7) So far we have examined rigid body motion, free vibration, structural vibration, and the coupling of flexibility to mechanism behavior. All of these must perform correctly to have high confidence in the formulation. In all of these problems, however, we have not examined the coupling among angular velocity and deformations due to the cross-coupling" terms in the equations. Up to this point, these terms were always zero. As was stated in Chapter V, one formulation [181 dismisses these terms entirely in every case. An experiment already referenced in Chapter V was run to examine the effect of these terms and other higher order terms.

115 SL I DER-CRANK 0.010 %. P-2 0 **ee.M Z - u 0 J -e.ees Z a. 8. -0.010 -e.*is TIMECSEC/SEC) Figure 7.4 Center point deflection of connecting rod for slider-crank (From DYMOS) Since DYMOS was formulated without the inclusion of these terms,3 a separate experimental program was used to test this area. In this case, the problem was formulated using the implicit second order form of the equations. The terms in the equations were assigned levels as shown in Tables 5.1 and 5.2. Each level was run and the results compared. It was found that significant changes occurred going from level zero to level one-half. But from level one-half to higher levels, no changes of significance 3All problems run on DYMOS used linear beams as flexible members. For transverse vibration of planar beams all of the cross coupling terms are zero, i.e. the Tisserand and Buckens frame are identical.

116 were observed. The change from level zero to level one-half can be seen in Figures 7.5 -7.8. The angular velocity maintains a constant value of 10 radians/sec for level zero equations (Figure 7.5). This is caused by the rotational equation being independent of deformation at this level. For level one-half, an oscillation appears in the solution for angular velocity (Figure 7.6). This is due to the expanding and contracting of the body represented by the third deformation mode. At further levels the angular velocity is only slightly different from level one-half. Ille 1 o.4 a I< 0 18.2 1.e X I I - I I I I I I I IIL= IIIIIIIII I I. I I I II I J II.... I i 0.8 0.6 9.4 I 0. 9.9 8.2 0.4 9.6 TIME (sec) 9.8 1.9 1.2 1.4 Figure 7.5. Angular velocity of triangle for level zero

117 The first deformation oscillates at its own natural frequency at level zero (Figure 7.7). This is because there is no coupling to other deformations at this level. For level one-half a fundamental change in the nature of the solution occurs (Figure 7.8). This change is due to interaction with mode two through a Coriolis type of effect. Further increases in level show virtually no further changes in the solution of the deformation. This experiment validated that there is coupling between angular velocity and deformations. The only terms of importance in this coupling are defined in Tables 5.1 and 5.2 as level one-half. II a a 4 Io= - - a I - I I I I I I I I I I I I I I 7..a.2 0.4 6.,,.6 TIME (sec) 1.8 1.2 1.4 Figure 7.8. Angular velocity of triangle for level one-half

118 0.1i 81o8 8 0s. a5 8.o0 -0. aS -4. 10 e. is.5a S.2 8.4 o.a 9. 1. 1.2 1.4 TIME (sec) Figure 7.7. First mode of triangle for level sero

119 @.s 1 ' ~\ I _, -. I 6.1 E / \.. *..02 *.4 *.* *.* 1. 1.2 1.4 TIME (sec) Figure 7.8 First mode of triangle for level one-half At this point all of the major effects that the equations represent have been accounted for and validated. Both formulations - explicit and implicit - have been successfully used for planar problems. The one remaining test was to see if the formulation was correct for three dimensional problems. This was done by running at least three specific problems: 1. De-stabilization of a free beam, spinning around the length axis, by transverse beam vibration. 2. A three dimensional pendulum, consisting of a rigid body and a spherical joint. 3. A three dimensional pendulum, consisting of a flexible body and a spherical joint.

120 Although these were not compared with the results of other programs, nothing was observed that indicated difficulty in the formulation. From the test cases cited in this section, all of the areas for validation were tested. In all of the tests, the formulation gave good results. In the case of the boundary condition test, it was seen that many free modes would have to be carried to approximate the cantilever and other fixed end condition beams. This can be remedied by using modes derived from these end conditions and then made orthogonal to the rigid body modes. 7.3. Conclusions The equations of motion for both linearly elastic members and rigid members were formulated in the preceding chapters. They were simplified as much as possible to provide efficiency in the solution process. The techniques used to simplify the equations were: 1) A local frame of reference was associated with each body. In this frame the deformation retained its linearity, despite the non-linear overall motion of the body. 2) The equations of motion for the rotation of the frame were formulated in terms of angular velocity instead of angular coordinates. This gave simpler vector equations. 3) By making the local reference frame "float" with respect to the body many of the terms in the equations were eliminated. This was done by choosing free modes for deformation coordinates, or could be done by constructing a set of coordinates with similar properties by using Gram-Schmidt orthogonalization. 4) Many of the remaining terms were shown to be of negligible size and were eliminated from the equations. The equations for various joints and forces were also formulated. Together with the equations for rigid and flexible bodies, these can be used by a computer program to automatically generate and solve the equations of motion for general mechanisms.

121 The equations were presented in two forms for numerical integration. One form was suitable for integrators requiring the equations to be explicitly solved for the first derivatives. This form was appropriate for non-stiff systems of equations. The equations were also presented in implicit form. This implicit form was more appropriate for stiff systems of equations. The implicit numerical integrator was extended to solve second order equations which further reduced the numerical effort required in the solution process. The formulation was validated using the test cases outlined in the previous section. The simplicity of the equations over other formulations holds great promise for high speed computation of flexible mechanism problems. This was not tested, but is anticipated in light of the fact that both the implicit method and the explicit method are more sparse and more compact than other methods proposed until this time. Another advantage of this formulation is that the terms in the equations have readily identifiable meanings. This gives the analyst a better understanding of the dynamics for any particular problem. A further advantage is that finite elements are not specifically associated with the formulation. They are just one method for obtaining the properties of the flexible body. Other techniques such as modal testing, or assumed modes can also be used to generate the desired properties needed for this equation formulation.

APPENDICES 122

123 APPENDIX A THE SECOND ORDER IMPLICIT INTEGRATOR SUBROUTINE DRIVE(NE NNZ) IPLICIT REAL8(A- HO - Z) C C THIS SUBROUTINE WILL INTEGRATE IMPLICIT DIFFERENIAL EQUATIONS C OF IORDER (2 OR 1) C SYSINT CONTAINS THE INTEGRATION PARAETERS C STATE HOLDS THE LOCATION OF THE STATE ARRAYS C COMMON /DDATA/ MD DIMENSION DD(1), MD(1) EQUVALENCE (MD(1),DD(l)) C COMMON /SPACE/ MSIZE MAVAIL, MAXNUMK MEAS, ERFLG, MBACK LOGICAL ERRFLG C COMMON /DEVICE/ NIN, NOUT, NERR, NTP, NDALSr, NTYIN COMMON /DUMP/ INDMP, ISTDMP, ICDM, IANDU', IOTDNP C COMMON /YSINT/ TSTART, TEND, HNIT, ESI, NSTEPS, IORDER, F COMMON /GRUSED/ SE NQUSED, NINTS, NFE, NJE, NBSUB NOPS COMMON /INTVEC/ IYY, YS, IYMAX, IERROR C COMMON /EVALU/ IG, NG, IRS, IDY, NEQ, INP0, NP, NPERR NRR COMMON /STATE/ IY IYD C COMMON /OUTD/ ITIM, OUT, NOV, NOUT LENPR C DIMENSION HINV(2) C WRITE (NERR,0) C(~ GET SPACE FOR ARRAYS NEEDED CALL GSPACE('DOUBE, ITIM, NSEPS+1, 1) CALL GSPACE('DOUB1, IOUT, NOV, NSTEPS-+) CALL TEMP(I, IMRIN1, IMRK2) CALL GSPACE('DOUB', IYS, NEQ, 7) CALL GSPACE('DOUB', IYY, NEQ. 7) CALL GSPACE('DOUB', IYMAX, NEQ 1) CALL GSPACE('DOUB', IERROR, NEQ, 1) CALL TEMP(2, IMRKI, IMRK2) IF (ERRFLG) CALL TERMDRIV', ', 1) C C,* SET UP EQUATIONS FOR INTEGRATION HNV(1) =.ODO / HINT HINV(2) = 2.0DO / (HNIT * HINT) CALL SETEQ(-1,.FALS&E,NENZDD(IYI),DD(IYD),TSARTHINV,-1.D0) C C* CALCULATE OUTPUT STEP SIZE HOUT - (TEND- TSTART) / NSTEPS TOUT = TSTART KERR =-1 NQUSED = 0 NOPS - 1 CALL OUPUT(TOUT, DD(nI), DD(IYD) NEQ DD(TIM), DD(IOUT), 1 NOV, NOPS) C Cm MAIN INTEGRATION LOOP. NGRATE IS REPETIVELY CALLED UNTIL TEND. DO 10 1 = 1, NSTEPS TOUT = TSTART + HOUT * I CALL NGRATE(IORDE, HNIT, TSTART, TOUT, NEQ, DD(nIY, DD(IY), 1 DD(II), DD(YD), DD(IYMAXn DD(IERROR) MD(INPERR), NRR,

124 2 DC(IRHS), EPSI, KE, MF) IF (KERR LT.- 1) GO TO 20 NOPS NOPS + 1 10 CALL OUTPUT(TOUT, DD(m), DD(IYD), NEQ. DD(TIM), DD(IOUT, 1 NOV, NOPS) C 0c ALL DONE (OR QUIT) 20 CALL TEMP(S, IMRK1, IMRK2) WRITE (NERR40) MAXNUM NOUTIS = NOPS RETURN 30 FORMAT (X/' ENTERING INTEATION PHASE...') 40 FORMAT (' INTEGRATION PHASE COMPLETE/'/' MAXNUM, 18) END SUBROTINE NGRATE(IORDER, HINIT, TZ TOUT, N, YY, YS, YI, YD, 1 YMAX, ERROR, NPERR NRR, REP,,KERR MT) IMPLICIT REAL*8(A - HO - Z) DIMENSON YY(N,1), RO S), 1), 1YMAX), NPERR(1) Y(l), 1 YD(N,2), YS(N,1) DIMESION HINV(2) COMMON /GRUSED/ HSED, NQUSE, NINS, NFE, NJE, NBSUB NOPS C COMMON /DEVICE/ N, NOUT, NERR, NTM, NDALST, NTYIN COMMON /DUMP/ INDMP, IS'DMP, ICDNP, IANDP, ICrTDh LOGICAL HCUT C C THIS SUBROUIINE INTEGRATES FROM TIME TO C TOUT BY CALLING NISIEP TO INTEGRATE UNTIL C IT GETS TO TOUT. C C MHIT HOW TO OBTAIN OUTPUT VALUES C 1 - NTERPOLATE TO OUPUT C 2- HIT OUTPUT TIME EXACTLY C CO INTIALIZE IF FIRST CALL IF (KERR.NE - 1) G TO 40 TIME ='TZ HOUT = TOUT TZ HI =.D0 / HOUT IF (MHIT.EQ. 2) HMAX =HOUT IF (MHIT.NE. 2) HMAX - 10.D0 HOUT HMIN HOUT * 1.0D>4 NQNEW = IORDER HNEW = HINIT HCUT.FALSEF HINV(1).ODO / HINIT HINV(2) 2.0DO / (HINITHINT) C*Z CONVERT YI,YD TO YY DO 10 I 1, N 10 YY(I) YI(I) DO 20 J -1, IORDER DO 201 =1, N 20 YY(I,J + 1) = YD(I,J) / HINV(J) C* FILL YMAX WITH ONES FOR ERROR CHECKING DO 30 I-1, N 30 YMAX(I) 1.ODO C Cs FOR EXACT OUTPUT HIT, JUMP TO OTHER SECTION C 40 IF (MHIT.EQ; 2) GO TO 90 C C* * INTERPOATION SECTION, C0 DO WHILE INTEGRATION TIME LESS THAN OUTPUT TIM 60 IF ((TIME - TOUTI)HI.GE 0.ODO) GO TO O8 C CALL NTSTEP(NQNEW, HNEW, TIME, N, YY, YS,YD, YMAX ERROR, NPERR, 1 NRR, RHS, HINV, HMIN, HMAX, EPSI, KERR) CO* WRITE OUT VECTORS AT END OF STEP (DEBUG) IF (IANDMP.LT. 6) GO TO 70 WRITE (NOUT,230) TIME, HUSED, NQUSED WRITE (NOUT,250) (YY(J,1),J1,N) DO 00 1 1, IORDER 0 WRITE (NOUT,260) (YD(JI),J==,N) C*. ERROR CIHECKING 70 IF (KERR.LT. - 1) GO TO 180

125 GO TO 60 C C** INTERPOLATE BACK TO OUTPUT TIME 80 CALL INTERP(TOUT, TIME HUSED, Yn, YD, IORDER, YY, NQUSED + 1, N) GO TO 1e0 C C * * EXACT HIT SECION * * * go IF (HCUT) HNEW HWAS C C* DO WHILE INTEGRATION TIME LESS THAN OUTPUT 100 IF ((TIME - TOUTHI.GT. -.01O0) GO TO 140 HCUT -.FALSE C- CAN HNEW BE STIRECHED 10% TO HIT oPUTI? JUMP IF NOT. IF ((TIME + 1.1DO*HNEW- TOUr)*HI.LE O.D0) GO TO 110 C- WILL THE CURRENT VALUE HIT OUTPUT WITHIN 1%? IF SO JUMP. IF (DABS((TIME + HNEW- TOUT)HI).LT..01DO) GO TO 110 HCUT.TRUE. HWAS - HNEW HNEW - TOUT TIME C 110 CALL NTEP(NQNEW, HNEW, TIME, N, YY, YS, YD, YMAX, ERROR NPERR, 1 NRRHS, INV, HMN, HMAX EPSI, KERR) C*- WRITE OUT VECTORS AT END OF STEP (DEBUG) IF (IANDMP.LT. 5) GO TO 130 WRITE (NOUT,230) TIME HUSED, NQUSED WRITE (NOUT,260) (YY(J,1),J=1N) DO 1201 1, IORDER 120 WRITE (NOUT,20) (YD(JJ,),J,N) O* ERROR CIIOKING 130 IF (KERR LT. - 1) GO TO 180 GO TO 100 C 0C LOAD OUTPUT ARRAY 140 TOUT TIME DO 160 I 1, N 160 YI( I)= YY,1) IF (IANDMP.GE 6) RETURN C C * ** RETURN SECTION 180 IF (IANDMP.LT. 4) REURN WRITE (NOUT,240) TOUT, HUSED, NQUSED WRITE (NOUT,260) (YI(J),J=l,N) DO 1701 I 1, IORDER 170 WRITE (NOUT,260) (YD(JJ),J,N) RETURN C C * * * ERROR SECTION * * * * 180 WRITE (NERR,190) TIME, HUSED, NFE, NJE, NINS IF (KERR.EQ - 2) WRITE (NERR,2O) IF (KERR.E. - 3) WRITE (NERR,10) IF (KERR.EQ - 4) WRITE (NERR,220) RETURN 100 FORMAT (' ERROR AT T=, 1PD124, ' Hm, D12.4, ' NFE, I6, 1 ' NJE=, 1, ' NINTSI, IS) 200 FORMAT (' CORRECTOR FAILED FIVE TIMES) 210 FORMAT (' REQUED STEPSIZE SMALLER THAN HMN') 220 FORMAT (' REQUESED ERROR TOO SMALL FOR PROBLEM) 230 FORMAT ('AT STEP: T=', 1PD12.4, ' H, D12.4, ' NQ, 13) 240 FORMAT ('AT OUTPUTIl. T, 1PD12.4, ' H=, D12.4,' NQ, IS) 260 FORMAT (IX 1P10D12.4)END SUBROUTINE INTERP(TOUT, TIME HUSED, YI YD, IORDE, YY, KDER, N) IMPLICIT REAL(A - H,O- Z) C C INTERPOLATES THE NORDSIEK VECTORS TO THE OUTPUT TIME C FOR BOTH THE VALUES AND DERIVATIVES C DIMENSION YI(1) YD(N,1) YY(N,1), SC(132)X HINV(2) C DATA SC /0.D0, 1.D0, 2.DO, 3.D0, 4DO, 6.DO, 8.D0, 7.0, 8.D0, 1 9.DO, 1O.DO, 11.Do, 12.Dz, 0.D0, 0.D0, l.D, 3.nD, 8.DO, 2 10.DO, 16.DO, 21.DO, 28.D, 38.D0, 46.D5, 66.DO, 6.DO/ C

128 C0 INITIALIZE THE ARRAYS DO 10 I 1, N 10 YI(I) = YY(I,1) DO 20 KK = 1, IORDER D0201 =1, N 20 YD(I,KK) =.DO C CG CALCULATE SCALE AND H VALUES H = TOUT- IME HINV(1) = 1.D / H HINV(2) =2.DO / (H*H) S -H / HUSED C 0C CLCULATE NEW VALUES DO 40 1= 1, N S1 = 1.DO DO 40 J 2, KDER S1 = S * S YN = YY(IJ) * Sl YI() YI(I) + YN DO 30 KK - 1, IORDER 30 YD(II) = YD(I,KK) + YN* SC(JKK) 40 CONTINUE C C0 CONVERT TO USCALED DERIVATIES DO 60 KK 1, IORDER DO60 1 =1, N 60 YD(I,K) = Y(KK) * HINV(KK) C RETURN END SUBROUTINE NISTPNQN, NEW, TIME N YY, YS, YD YMAX ERROR, 1 NPERR, NRR, RHS, HINV, HMN, HMAX EPSI KRR) IMNLCIT REAL*(A - OH,O- Z) C C THE PARAM1ERS TO THE SUBUROTINE N1STEP HAVE C THE FOLLOWING MEANINGS.. C C NQNEW ORDER TO BE USED ON THE NEXT STEP C HNEW THE STEP SIZE TO BE ATTEMPD ON THE NEXT STE C EVAL NAME OF SUBROUTNE WHICH EVALUATES THE JACOBIAN AND RI3 C TIME LAST SUCCESSFUL TIME C N THE NUMBER OF DIFFERENTIAL EQUATIONS. C YY AN N BY 7 ARRAY CONTAINING THE DEPENDING VARIABLES1 C AND THEIR SCALED DERIVA7TVES Y(IJ+) CONTAINS C THE J-TH DERIVATIVE OF Y(I) SCALED BY C H**J/FACTORIAL(J) WHERE H IS THE CURRENT C STEP SIZE C YD ARRAY WITH THE UNSCALED FIRST AND SECOND DERVATIVES C YMAX AN ARRAY OF N LOCATIONS WHICH CONTAINS THE MAXUM C OF EACH Y SEEN SO FAR. IT SHOULD NORMALLY BE C SET TO 1 IN EACH COMPONENT BEFORE THE FIRST ERY C NPERR VECTOR CONTAINING INTEGERS INDICATING THE DiFFERENTIA C ORDER OF THE CORRESPONDING EQUATION C 0 ALGEBRAIC C 1 UPPER FIRST ORDER (NO DOUBLE DOT ICS) C 2 2ND ORDER AND LOWER 1ST ORDER C ERRSV AND RHS ARE INTERNAL WORK VECTORS C YMAX,ERROR,ERRSV,AND RHS MUST BE DIENSIONED N OR GREATER C ERROR AN ARRAY OF N ELEMENI S WHICH CONTAINS THE ESIMATED C ONE STEP ERROR IN EACH COMPONENT. C HMIN THE MINIMUM STEP SIZE THAT WILL BE USED FOR THE C INTEGRATION. NOTE THAT ON STARTING THIS MUST BE C MUCH SMALLER THAN THE AVERAGE H EXPECTED. C HMAX THE MAXIMUM SIZE TO WHICH THE STEP WILL BE INCREASED C EPSI THE ERROR TEST CONSTANT. SINGLE STEP ERROR ESTIMATES C DIVIDED BY YMA(I) MUST BE LESS THAN THIS. C THE STEP AND/OR ORDER IS C ADJUSTED TO ACHIEVE THIS C KERR A CODE WITH THE FOLLOWING MEANINGS. C 0 STEP SUCCEEDED BUT ENCOUNTERD TROUBLE C 0 THE STEP WAS SUCCESFUL. C -1 INITIALIZE SUBROUTINE C -2 CORRECTOR FAILED 6 TIMS C -3 REQUIRED STEPSIZE SMALLER THAN HMIN C -4 THE REQUESTED ERROR IS SMALLER THAN CAN

127 C BE HANDLED FOR THIS PROBLEM. C C INTERNAL VARIABLE MEANINGS ARE AS FOLLOWS: C MAXDER THE MAXIMUM DERIVATIVE THAT SHOULD BE USED IN THE C METHOD. SINCE THE ORDER IS EQUAL TO THE HIGHEST C DERIVATIVE USED, THIS RESICTS1 THE ORDER. C T INTEGRATION TIME C MINDER THE MINIMUM ACTIVE OOLUMN LENGTH OF YY ARRAY C H CURRENT STEPSZE C NQ CURRENT ORDER OF METHOD C NXTCHK PREVENTS.ORDER CHECKING FOR A NUUMER OF ST-EPS C ITFAIL KEEPS TRACK OF CORRECIOR FAILURES C INTFAL KEEPS TRACK OF INTEGRATION ERROR FAILURES C ZIP IF ABS(RS) IS ESS THAN ZIP, ITS CLOE ENOUGH TO ZERO C C DIMENSION YY(N,1), ERROR(1), RH(1), YMAX(1), NPERR(1), YD(N,2), 1 YS(N,1) C COMMON /DEVIC/ N, NOUT, NE NTM, NALST, NTDIN N N COMMON /DUMP/ INDMP, ISTDMP, ICDMP, IANDMP, IOTDMP C COMMON /GRUSED/ HUSED, NQISED, NINIS, NFE, NJE, NBSUB, NOPS C LOGICAL CONVRG, RZERO C C DIMENSION A(13), ER3S) DIMENSION KMODE(10), HINV(), GSC2) DATA KMODE /2, 1, 1, 1, 2, 1, 1, 1, 2, 1/ DATA ITMAX /10/ C C THIS SUBROUTINE INTEGRATES ONE STEP IN TIME C C Co IF FIRST CALL INITIALIZE SUBROUTINE IF (KERR.NE - 1) CO TO 10 IORDER NQNEW MINDER - IORDER + 1 MINDP1 = MINDER + 1 KCON =-1 H -HNEW NQ =0 MAXDER 7 FNRR = NRR CALL COSET(NQ, KDER, NQNEW EPSI, GSC, RSC, HINV, NXTCHI KERR, 1 TIME FNRR,, ER) IDBUG = IANDMP - ZIP =- EPSI * 1.OD-3 BND -EPSI * 1.D1 NINTS 0 NFE= 0 NJE =0O NBSUB 0 C C* * INITIALIZE VARIABLES FOR STEP l * C 10 KERR 0 ITFAIL = O INTFAL = 0 C* SAVE NORSIEK VECTORS AND ASSOCATED H IN CASE STEP FAILS DO 20 J 1, KDER DO20I =1, N 20 YS(IJ) YY(I,J) HOLD = H C C* * GET READY FOR SaTEP * C C* SET UP POSSIBLE NEW ORDER 30 IF (NQNEW.NE NQ) CALL COSET(NQ KDER, NQNEW, EPSI GSC, RSC, 1 HINV, NXTCHK, KERR, TIME FNRR, A, ER) IF (KERR.LT. 0) RETURN C C* RESCALE NORDSIEK VECTORS IF STEPSIZE CHANGED 40 H = HOLD IF (HNEW.NE H) CALL RESCAL(H. HNEW, YYYS, N, KDER, HINV, GSC, 1 A, TME HMIN, HMAX, NXTCHK KERR) IF (KERR.LT. 0) RETURN

128 C* * * PREDICT NEW VALUES * * * C T ==TIE+H C CO NORDSIEK VECTORS DO 60 J 2, KDER DO 60 J1 =J, KDER J2- KDER- J1 +J- 1 DO60 I =1, N YY(I,J2) YY(I,J2) + YY(,J2 + 1) 60 CONTINUE C C0 COMPUTE PREDICIED DERIVATIVES FOR EVALUATOR DO 80 J _ 1, IORDER DO0 I - 1, N 60 YD(lJ) = YY(J + 1) * HINV(J) C C* ZERO ACCUMULATED ERROR VECTOR DO 70 1 - 1, N 70 ERROR(I) = 0.ODO C C *~ ITERATE ON CORRECTOR * C DO 110 L = 1, ITMAX C 0C SOLVE EQUATIONS FOR DY (RETURNED IN RHS) IF (IDBUG.GE 2) WRITE (NOUT,240) CALL SOLVE(KMO E(L) KCON T, RH YY, YD, N, RSERR, GSC, RSC, 1 RZERO, ZIP, IDBUG) IF (RZERO) GO TO 150 IF (IANDMP.GE. 6) WRITE (NOUT,230) IF (IANDMP.GE 6) WRITE (NOUT,26) (RHS(K),KK-,N) C C*. CHECK FOR DIVERGENCE, IF SO JUMP IF (L.GE 4.AND. RHSERR.GT. 2*RHEOLD) GO TO 120 RHEOLD - RHSERR C C0 UPDATE STATE AND DERIVS, ACCMULATE ERROR, CHECK CONVERENCE CONVRG =.TRUE, DO 100 1 =1, N DO 0S J 1, MINDER 80 YY(I,) = YY(I,J) + A(J) RHS() DO 90 J 1, IORDER 90 YD(IJ) = YY(I,J + 1) * HINV(J) ERROR(I) = ERROR(I) + RS(I) YMAX(I) = DMAX1(YMA)I),DAB(YY(I,l))) IF (DAS(RHS(I))/YMAX(I).GT. BND) OONVRG -.FALSE 100 CONTINUE IF (OONVRG) GO TO 160 110 CONTINUE C O * CORRECTOR ITERATION FAILED * 120 KERR 2 IFAIL =1TFAIL +1 IF (IANDMP.GE. 1) WRITE (NERR200) T, H C Os IF CORRECTOR HAS FAILED 6 TIMES GIVE UP IF (ITFAIL.LT. 6) GO TO 10S KERR — 2 REITURN C 0* DETERMINE PROBABLE CAUSE OF FAILURE 130 IF (RISERR.GT. EPSI*FRR) GO TO 140 C C*0 ILL-CONDTIONED MATRIX HNEW. 10.ODO * H HMAX = DMAXI(HMAXJ-IEW) GO TO 40 C CO STEPSIZE TO BIG 140 HNEW - H / 4.0DO GO TO 40 C C * * * CO)RRECTOR CONV GES * * * * C C* IS THIS A BLIND STEP (NO ERROR CONTROL)L

129 160 CONTINUE C IF(.NOT. INTCHK) GO TO 210 C Co TEST ERROR D = O.ODO DO 180 I = 1, N IF (NPERR(I) NE 2) GO TO 100 D - D + (ERROR(I)/YMAX(I)) - 2 160 CONTINUE IF (D.L ER(2)) GO TO 180 C C * INTEGRAION ERROR TOO BIG. * *.* C INTFAL INWAL + 1 KERR =IF (IANDlPW.GE 1) WRTE (NERR,280) T, D, ER(2) C C. MANY INTEGRATION ERRORS REQUIRES DRASTIC ACTION IF (INTFAL.LT. 6) GO TO 170 NQNEW = IORDER HNEW H / 10.ODO GO TO 30 CO FIND BETTER ORDER AND STEPSIZE 170 CALL HORDER(NQNEW, HNEW, NQ, H D, KDER, MAXDER, MNDER, HMAX, YY, 1 YMAX, ERROR NPERR, N, T, NXTCK, KERR, A, ER) GO TO 30 C C*** INTEGRATION ERROR OK ** C C0* UPDATE THE REMAINDER OF THE NORDSIEK VECTORS 180 CONTINUE IF (KDER.LT. MINDP1) GO TO 200 DO 190 J = MINDP1, KDER DO 190 I -1, N 190 YY(I,J) = YY(I,J) + A(J) * ERROR(I) C C- DETERMINE IF ORDER, SEPSIZE CHECKING IS DESIRED 200 IF (NXTCHK.GT. 1.OR. (H.GE HMAX*.99 AND. L.LE 3)) 1 GOTO 210 C C- CHECK ORDERSTPSZE CAL HORDER(NQNEW, HNEW, NQ H D, KDER, MAX MINDE, HMAX, YY, 1 YMAX ERROR, NPERR, N, T, NXTCYK, KERR, A, ER) C C* PREVENT FURTHER CHEK(ING FOR 10 SEPS? IF (KERR.NE. 0.OR. NEW.GT. 1.1H.OR. L.T. 3) GO TO 210 NQNEW - NQ HNEW H NXTCHK =-10 C * * RETURN S ON ** * C 210 NXIICHK == NXTCHK -1 TIME - T HUSED = H NQUSED = NQ NINTS = NINTS + 1 IF (IANDMP.GE. 2) RTE (NERR,270) T, HUSED, NQUSED, L, NINTS, 1 NFE, NJE, NBSUB NOPS IF (NXTCHK.GT. 1.OR KDER.GE MAXDER) RETURN C C* SAVE ERROR FOR ESIIMATES IN HORDER DO 220 I 1, N 220 YY(I,MAXDER) = ERROR(I) RETURN C 230 FORMAT (' DELTA X S AT ITERATION:') 240 FORMAT (' DIFSUB JACOBIAN AND RHS:') 250 FORMAT(' DX: ', 1P6D1&./61i(T9,D16.6/)) 200 FORMAT (' CORRECTOR FAILED TO CONVERGE. TRYING NEW STE SIZE/ 1. ' TIME, 1PD12., ' H-, D12.6) 270 FORMAT (' INTEGRATION STATUS AT T', 1PD12.4, /' HUSED=-, D12.4, 1 ' NUSED==, 12,' NCORm', 12/ NINTSIm, 16,' NFE=, 16, 2 ' NJE,, 1 NBSUB, 1, ' NOPS'S, 16) 280 FORMAT (' INTEGRATION ERROR FAILURE'/' T=a, 1PD12.4,' ERRORm', 1 D12.4, ' ERRLIM, D12.4) END

130 SUBROUTINE COSET(NQ. KDR, NQNEW, EP, GSC, ISC, HINV, NXTCHK, 1 KERR TIME; FNRR, A ER) IMPLICIT REAL(A - HO - Z) C C SET THE COEFICIETS OF INEfGRATOR ACOORDING TO ORDER C DIMENSION PE'TST(12,2,3), GSC1), HINV1) COMMON /DEVICE/ NIN, NOUT, NERR, NTMP, NDALST, NT IN COMMON /DUMP/ INDMP, ISTDMP, ICDNP, IANDNV, IOTDMP C DIMENSION A(1), ER(1) C DATA PERTST / 1.,1.,2.,1.,.3168.07407,.01391,.002182, 1.0002946,.00003492,.000003892,.0000003624, 2 1.,1.,.,.1687,.04167,.00827,1.,1.,1.,1.,1.,l., 3 2.,12.,24.,3.89,,63.33,70.08,7.97,106.9, 4 125.7,147.4,188.8,191.0, 6 2.0,45,7.=333,10.42,13.7,,7.15,1,11,,1,1.,1 8 12.0,24.0,37.89,6333,70.08,87.97,10e.9, 7 128.7,147.4,188.8191.0,1., 8 3.0,0.0,a17,12.6&l,16,l,lll.,1l1.lll. / C IF (IANDNM.GT. 2) WRITE (NERt0) TIME, NNEW, NQ NQ == NQNEW KDER = NQ+ 1 GO TO (10, 20, 30, 40, 0, O0), NQ C 10 A(l) -1.0DO A(2)= -i.D GO TO 70 C 20 A(l) =-0. 6DO A(2) — 1.D A(3) — 0.334,343,343644333DO GO TO 70 C 40 A(l) = -0.4oooooo4444 A(3) =.A(l) A(4) = -0.09o90909090910o GO TO 70 C 40 A() =-0.48000000047 A(3) =- -0.70000000087 A(4) = -0.20000000028 A(6) 0.00000000744 4744 GO TO 70 C 60 NA) =-0.4370620437OI 2D0 A(3) = -0.82113832 160788D A(4) = -0.310218981021D89 A() = -0.064744626474462 DO A(N) -0.0036496360493604)D0 GO TO 70 C 80 A(l) -0.408183 M1, 220N A(3) = -0.92003492063492D O0 A(4) _ -.418888MOMMMe(7o A(6) -0.09920634920/3492DO A(i) = -0.0119047819047619DN A(7) =8-0.0006=9342400282DO C 70 NXDTCH KDER + 1 DO 80 I =-1, 3 80 ER(I) = FNRR * (PERTST(NQAI)EPSI) *0 2 RSC= -1l.ODO / A(l) GS(1) -FHINV(1) / A(1) GSC(2) A(3) HINV(2) / A(1) C IF (ER(2).EQ. 0.0DO) KERR -=-4 RELURN 90 FORMAT (' C06ET. TIMa, 1PD12.6, 'NQNEW- 12, 'NQ=i, 12) END SUBROUTINE HORDER(NQNEW, HNEW, NQ, H D, KDER, MAXDER, MNDER, 1 HMAX, YY, YMAX ERROR, NPERR, N, TIE, NXCrCHK, KERR, 2 A, ER) IMPICIT REAL*8(A- H, - Z)

131 C C0 CALCULATE FACTORS H COULD BE DVIDE BY FOR ADJACENT ORDERS DIMENSION YY(N,1), YMAX1), ERROR(1), NPERR(1) C DIMENSION A(I), ER(1) C COMMON /DEVICE/ NN, NOUT, NERR, NlTP, NDALST, NTYIN COMMON /DUMP/ INDNM, ISDMP, ICDMP, IANDMP, IOTDMP C C 0C* THIS ORDER ENQ2 -.DO / (NQ + 1) PR2 = ((D/ER(2))*ENQ2) * 1.2 C C* ORDER INCREASE PR3 = l.OD+20 IF (KDER-.GE. MAXDER.OR KERR NE 0.OR. NXTC.GT. 1) 1 GO TO 20 D 0.ODO DO 10 I =1, N IF (NPERR(I) EQ. 0) GO TO 10 D = D + ((ERROR(I) - YY(IMAER))/YMAXI)) 2 10 CONTINUE ENQ =.BDO / (NQ + 2) PR3 = ((D/ER(3))*EN3) * 1.4DO C C0 ORDER DECREASE 20 PR1 = 1l.OD+20 IF(KDER.LE MINDER) GO TO 40 D ==.ODO DO 30 1= 1, N IF (NPERR(I).EQ 0) GO TO 30 D = D + (YY(I,KDER)/YMAX(I)) 2 30 CONTINUE ENQI =. DO / NQ PR1 = ((D/ER(1))*.ENQ1) * 1.3D C 40 CONlINUE PRIT - l.ODO / DMAX1(i.OD10,PRI) PR2T = l.ODO / DMAX(1.OD.10,PR2) PRST l.ODO / DMAXl(I.0D-10,PR3) C COn DEERMINE BEST ORDER IF (PR2.LE PRS AND. PR2.LE PRI) GO TO 60 IF (PR3 LE PR2 AND. PR3.LE PR1) GO TO 60 C C* REDUCE ORDER NQNEW=NQ- 1 R = 1.ODO / DMAX1(PR1,I.OD-4) GOTO80 C 0C SAME ORDER 60 NNEW NQ R = l.ODO / DMAX1(PR2,1.0D.4) GOTO 80 C C- INCREASE ORDER 60 NQNEW = NQ + R =.ODO / DMAXI(PR3,1.0D,4) C! ESTIMATE NEXT DERIVATIVE DFLK = KDER DO 70 I _ 1, N 70 YY(I,NQNEW + 1) ERROR(I) * A(KDER) / DFLK C 80 CONTINUE HNEW = HMAX / IFIX(SNGLHMAX/DMN1(HMAXH*R) + 0.DO)) IF (IANDMP.GE 3) WRITE (NER,90) TIME, HNEW, H PRIT, PR2T, PRST RETURN 90 FORMAT (' ORDER: TIMEI=, 1PD12.6, ' HNEW, D12.6,' H=, 1 IPD12.6, /' RD=a, IPD.2, ' R=a, 1PD.2, ' RUm, IPDQ.2) END SUBROUTINE RESCA4(H, HNEW, YY, YS, N KDER, HINV, GSC, A, TIME 1 HMN, HMAX. NXTCMK, KERR) IMPLICIT REAL*(A- HO - Z) C C RESCALES THE NORDSIEK VECTORS WITH NEW ST SIZE

132 c COMMON /DEICE/ MN, NOUT, NERR, NTP, NDALST, NTTYIN COMMON /DUMP/ INDMP, ISTDMP, ICDMP, IANDPS, IOTDIP C DIMENSION YY(N,1), GSC(1), HINV(), A(l), YS(N1l) C 0C ERRCO REIURN IF REQUESTED TEPSIZE BELOW MNIMUM IF (HNEW.GE..99HMIN) GO TO 10 KERR =-3 RE7TURN C C0 RESITRE SAVED VEIORS IF ERROR, AND LIMIT S1r 10 IF (KERR.E. 0) GO TO 30 DO2 J 1, KDER DO201 1, N 20 YY(IJ) YS(IJ) IF(KRR.EQ 1) HNEW - DMIN1(HNEW, C 00 LIMTr STEPSIZE TO HMAX AND 100LD 30 HNEW - DMIN(HNEW,HMAX,10.0D*H) RH HNEW/ H C CI INCREASES MUST BE MORE THAN 10% IFRH.LT. 1.ODO.OR. RH.GE L1.DO) GO TO 40 RH = 1.ODO HNEW H IF (IANDMP.GE 3) WRITE (NERR,0) HNEW, H RH NCHK = KDER + 1 RETURN C C0 RESCALE NORD6IEK VECTORS 40 RI =' 1.0 DO 60 J 2, KDER Rl = RI * RH DO 60 I =1, N 60 YY(I,J) = YY(LJ) * RI IF (IANDMP.GE. 3) WRITE (NERR,) HNEW H, RH H - HNEW C C* PROHIST CHANGES TEMfORARILY NXTCHK - KDER + 1 C 0* UPDATE STEPSIZE DEPENDANT VARIABLES HINV1) =.ODo / H HINV(2) 2.0DO / (H*H) cGSq) =-HNV() / A(1) GS2) = A(3) * HINV(2) / A(1) RETURN 60 FORMAT (' SCALE: HNEWd, 1PD12.6,, LPD12., RATIo, 1 1PD126) END

133 APPENDIX B THE SPHERICAL JOINT CONSTRAINT DERIVATION The spherical joint has one vector constraint equation associated with it. The chain of vectors is shown in Figure B.1. The constraint equation is simply the sum of the vectors is always zero. ~ = A A BR + B) = 0 Only half of ths equation will be expanded because the second half is similar, but with opposite sign. A B.O <_. A+_ sB _ 0 In order to apply this constraint a Lagrange multipler X is introduced to enforce the constraint. To get the components of X for body A, the partials of the constraint are needed. To do this the variation of the constraint is performed. 6(A ) 6(A R) + 6 (A ) At this point the explicit reference to body A is dropped. The variation of q from Equation 3.8 is 8 =- 6r x q + 6, 6tl,. The variation of 4 has become 64 = 6R + 6r X q +, 6t1, In this case the vector represented by q does not range over the body, but is the vector from the origin of the reference frame to the location of the spherical joint on body A. Let this point be denoted by a. Therefore q = - (). The shape vectors 4, also, are

134 P.) F6 *R Figure B.1 Spherical Joint the shape vectors' value at point a. Therefore, =,i(3). The variation in these terms is S=. SR + Sn X q(a +, (a6l i The fundamental equation requires the dot product of X with the partials of the constraint equation. This is done by dotting X with the variation and collecting the partials.

135 -..-4 -4 4 -.* -. 4 -... x * 6 =X * bR + X * 6 X q() + X( * i,(7)&, = R * X + & 7 q(7 X X + ii X * () From this the terms are X = -- X OR X * _ q(;)X X The equations of motion become for body A: 6R *MR +X +6{J fl + n x. +,. fl + i() x + n, {e; + 2(,)(,), + (as, 2 2 nm Ai n + 2n G,,q + X, (a) -_ o Here we see that component affecting the position of the frame is the reaction force X. The component affecting frame rotation is the torque produced by this reaction force x(a) x X. The component affecting mode i, is the ith mode shape evaluated at a dotted with X. These could all have been expected. These same terms appear for body B, but with reversed sign. In addition ~ and i, are evaluated at a point on body B. In addition to the terms in the body equations the constraint equation itself must be included with the equation set. For the explicit integrator two derivatives are

138 taken of the constraint equation. 4= R + n xq(a) +. (1). -R+ x(n xfq(l)+ 0 xq(a +f xu XR f+n Xjq(4+a) + 2i(ais + n x (n x () + 2n x r, ()T, Now the various vectors can be resolved in frames of reference. The new vectors X and 4 will both be resolved in the inertial frame. The equations become: 6R T MR + ) 7 nJo + n Jo n + n + q(a)A X ai, (, + 2,)(,), + W(2),) 2 nr J, l + 2fnr GmI, +, (a)T A t The constraint equation can be written as * = R - q(a) xn +?, ();ti + x ( x * () + 2n x 4, ()l, = R -A"' g(a)n+A"' _i(a)Ai + Af n n 2(a) + 2Af' n le (a),

137 These equations can be written in a larger matrix form (including body B) as: JO(a) 1 M(b)I I -I 1 -1()A bf Abt)rT O R(a) now e00 n(a) R(b) n(b) (0 n1(b) X JO (b) e l-4")r (,(IA~) -1 (b)A bf )T -_ - _{A ' ( (b)T I * a - b * V The system mass matrix for the explicit form of the equations has the form: M rl \ L OJ where M is a strictly diagonal and constant matrix, and r has columns of non-linear terms, and z are the remaining terms in the constraint. In this case +' -(Al an. ) n(a) (a) - 2 AfA n2(,) (a)i,(a) +A/ nb) ^q(ba ) -2Afb nfi(b)kb)/b) For the implicit integrator a different approach is taken. In this case the constraint equation itself in algebraic form is included with the equations. The equations become MR+X - 0

138 4.. -, -,0 0 -. JO n + n XJo f n + J1kk l +q(a) x =- o.i + i )w(i + W(?15i 1 4 - - _ - n Jl, nfl+2G G,,,do,, +X, (a) 0 In order to integrate these equations the matrix of partial derivatives is needed. The partial derivatives for the terms associated with the body were done in Chapter V. All that remains is to find the partials of the terms corresponding to the newly added constraint terms. The easiest way to calculate these partials is to leave the equations in vector form and take the variation. The constraint equations variation was already performed (equation 6.5). Its partials are obtained by expressing it in terms of dyadics. 6= = I -8R - q 6 -x +Vi ((a)ti, Be resolving this vector in the global frame a matrix equation results. 64 = R - A1o q6 + Af" i 6q From this we obtain the partials in terms of matrices: 9 IOR - = -A' o q 0?, -- a -- One difficulty remains to be solved. The partials of the constraint are needed with respect to actual coordinates. Both R and,i are actual coordinates used in the equations. However, r is not a coordinate; its infinitesimal variation is only used to generate the equations of motion. The actual rotational coordinates are something different, in this case Euler parameters. What is needed then are partial derivatives of the constraint with respect to the Euler parameters.

139 The kinematic differential equations provide a simple solution to this problem. This equation set relates the time rate of change of angular orientation | d-t to time rate of change of angular coordinates I d-p. dt dr dp dt -dt If the equation's divisor of increment in time is eliminated this becomes dr = B dp or in terms of a partial derivative O-B 9p The partial derivative can be used to get the explicit dependence on the angular coordinates by using the chain rule: 9p Oi9 ap aG_ ap The partial of $ with respect to angular coordinates becomes dp - - Op The only remaining partials are those for the equations involving the lambdas. The extra term in the translation equation is just X. The partial of this is just the identity matrix again. The rotation equation has the term X(a) in it. The variation of this is 6(a) * X )= e() * X + + ) () q X6A The third term is the term already calculated for the explicit equations, occurring in r. The first and second term have not appeared before. Because these terms will result in the partials of the partials of the constraint, they are called the second partials of the constraint equations. The second term can be expanded as follows

140 e() * q X X = e() e (67r Xq + Oi (a)(i ) X X = -j(b) X X (r X q + ib 6, ) - e X (q* 7r - X X O, 6i, i X -- --. = (b) ~* nrn.e. a(b)r ) -() =A 4f XA f/. - A f XA fa, b, The first term in this expression can be further expressed in terms of angular coordinates using Equation 6.7: 6fr = B6p. This expression becomes A _ XAf1 _JB6 A a/ XA _ _ 7 _ Computing this expression can be simplified by recognizing that A 'XA4 is just X expressed in body coordinates. If we let X( A 1 X the triple matrix product becomes simply (a). The expression is now x(4)qBc 6 - x(W)O 61i, The first term has the variation of the unit vectors associated with the body. Since they are rotating with the body this variation is (a") - 6fr x t. This term becomes g(q x X ') (q X) = 6R *e() x (q XX) - (a) (q X X) X - qX(a) 67r = qX() 1 B6 This term can be combined with the other containing a variation in angular coordi

141 nates [IqX_( + X(a)q B6p From Identity 1.4.22 the term in brackets becomes qX. The variation in angular coordinates becomes qX(a)B 6E. The partials are: 0(qAax) r = =( qA qr X) p - - 071, ( qA X) = X(A ax - The partials for the term in the modal equation is done similarly. S(, X) 6= xr Xqi *X + i,6 =_ X X* 2 + _, *_ -. -9, -9 +, -= i X *6r +,i * aX The second term was again expressed previously in r. The first term is the term corresponding to the second partials. It must be expressed in the body frame. Using X(a) this becomes i X(a) a or in terms of angular coordinates 0 X(O)BSp The partials are (A (ax), A dX a(_, A ' _)

142 The matrix of partials for the spherical joint for the two bodies looks like: (4-)B -X(6) B -(')')B OF( )A '/ [)r A f' 1 (X(')x)Tf B (()r A If) -1 - (-I_ _A&1 -I (')A ' _,jts, sO I K' (." ()A f)r B (4,)r A 1)r

REFERENCES 143

144 REFERENCES [1] Erdman, A. G., Sandor, G. N., Bakberg, R. G. "A General Method for KinetoElastodynamic Analysis and Synthesis of Mechanisms," Journal of Engineering for Industry, ASME 71- WA/DE-6, New York, 1971. [2] Midha, A., Erdman, A. G., Frohrib, D. A., "A Computationally Efficient Numerical Algorithm for Transient Response of High-Speed Elastic Linkages," ASME 78-DET-54, New York, 1978. [3] Erdman, A. G., Imam Imadad, Sandor, G. N., "Applied Kineto-Elastodynamics," Paper presented at 2nd Applied Mechanism Conference, Oklahoma State University, Oklahoma. [4] Hossne, A. J., Soni, A. H., Harvey J. W., "Kineto-Elastodynamic Analysis of Mechanisms by Finite Element Formulations," ASME 80-DET-10S, New York, 1980. [5] Midha, A., Erdman, A. G., Frohrib, D. A., "A Closed-From Numerical Algorithm for the Periodic Response of High-Speed Elastic Linkages," Journal of Mechanical Design, ASME 78-DET-15, New York, 1978. [6] Naganathan, G., Wilhert, K. D., "New Finite Elements for Quasi-Static Deformations and Stresses in Mechanisms," ASME 80- WA/DSC-35, New York, 1980. [7] Song Ji Oh, Haug Edward, "Dynamic Analysis of Flexible Mechanisms," Materials Engineering Technical Report No. 55, The University of Iowa, May 1979. [8] Sunada, W., Dubowsky, S., "The Application of Finite Element Methods to the Dynamic Analysis of Flexible Spatial and Co-Planar Linkage Systems," Journal of Mechanical Design, ASME 80-DET-87, New York, 1980. [9] Dubowsky, S., Gardner, T. N., "Design and Analysis of Multilink Flexible Mechanisms with Multiple Clearance Connections," Journal of Engineering for Industry, Vol. 99, No. 1, February 1977.

145 [101 Shabana A., Wehage, R. A., "Variable Degree of Freedom Component Mode Analysis of Inertia Variant Flexible Mechanical Systems," Journal of Mechanical Design, ASME 82-DET-93, New York, 1982. [11] Benson, David J., The Simulation of Deformable Mechanical Systems Using Vector Processors, Ph.D. Thesis, University of Michigan, 1983. [12] Bodley, Carl, Devers, A. D., Park, A. Colton, Frisch, Harold, A Digital Computer Program for the Dynamic Interaction Simulation of Controls and Structure, Vol. 1, NASA Technical Paper 1219, NASA 1978. [13] Grote, P. B., McMunn, J. C., Bluck, R., "Equations of Motion of Flexible Spacecraft," Journal of Spacecraft and Rockets, Vol. 8, No. 6, June 1971. [141 Likins, P. W., Finite Element Appendate Equations for Hybrid Coordinate Dynamic Analysis, JPL Technical Report 32-1525, Pasadena, 1971. [15] Likins, P. W., Fleischer, G. E., Large - Deformation Modal Coordinate for Nonrigid Vehicle Dynamics, JPL Technical Report 32-1565, Pasadena, November 1972. [161 Likins, P. W., Dynamics and Control of Flexible Space Vehicles, JPL TR-32-1329, Pasadena, February 1969. [17] DeVeubeke, B. F., "The Dynamics of Flexible Bodies," Journal of Eng. Sci, Vol 14, pp 895-913, 1976. [18] Canavin, J. R., Likins, P. W., "Floating Reference Frames for Flexible Spacecraft," Journal of Spacecraft, Vol. 14, No. 12, December 1977. [19] Kane, T. R., Wang, C. F., "On the Derivation of Equations of Motion," J. Soc. Indust. Appl. Math., vol. 13, No. 2, June 1965. [20] Laskin, R. A., Likins, P. W., Longman, R. W., "Dynamical Equations of a FreeFree Beam Subject to Large Overall Motions," Journal of the Astronautical Sciences, October-December 1983. [21] Orlandea, N., Chace, M.A., Calahan, D. A., "A Sparsity-Oriented Approach to the Dynamic Analysis and Design of Mechanical Systems - Part 1, Part 2," Journal of Engineering for Industry, ASME 76-DET-19,20, New York, 1976. [22] Wittenburg, Jens, Dynamics of Systems of Rigid Bodies, B. G. Teubner, Stuttgart, 1977. [23] Marion, J. B., Classical Dynamics of Particles and Systems, 2nd Edition Academic Press, New York, 1970.

146 [24] Pars, L. A., A Treatise on Analytical Dynamics, Heinemann, London, 1965. [25] Greenwood, D. T., Classical Dynamics, Englewood Cliffs: Prentice-Hall 1977. [26] Sohoni, Ul, Haug, E. J., "A State Space Method for Kinematic Optimization of Mechanisms and Machines," Journal of Mechanical Design, ASME 8O-DET-43, New York, 1980. [27] Gear, C. W., "Simultaneous Numerical Solution of Differential-Algebraic Equations," IEEE Transactions on Circuit Theory, Vol. CT-18, No. 1, January 1971. [28] Dubowsky, S., Grant, J. L., "Application of Symbolic Manipulation to Time Domain Analysis of Non-linear Dynamic Systems," Journal of Dynamic Systems, Measurement, Control, ASME 75-Aut-J. [29] Witham, C. R., Dubowsky, S., "An Improved Symbolic Manipulation Technique for the Simulation of Nonlinear Dynamic Systems with Mixed Time-Varying and Constant Terms," Journal of Dynamic Systems, Afeasurement, and Control Transactions of the ASME, September 1977. [30] Cipra, R. J., Uicker J. J. Jr., "On the Dynamic Simulation of Large Nonlinear Mechanical Systems: Part I - An Overview of the Simulation Technique Substructuring and Frequency Domain Considerations," ASME 80, 1980. [31] Cipra, R. J., Uicker, J. J. Jr., "On the Dynamic Simulation of Large Nonlinear Mechanical Systems: Part II - The Time Integration Technique and Time Response Loop," ASME, 1980. [32] Case, P. L., Pratt, D. T., "A New Integration Algorithm for Chemical Kinetics," Presented to the 1977 Spring Meeting of the Western States Section/ The Combustion Institute, WSS/CI 77-23, 1977. [33] Gear, C. W., Numerical Integration of Stiff Ordinary Differential Equations, Report No. 221, Department of Computer Science, University of Illinois, January 1967. [34] Gear, C. W., "The Automatic Integration of Ordinary Differential Equations," Communications of the ACM, Vol. 14, No. 3, March 1971. [35] Calahan, D. A., Computer-Aided Network Design McGraw Hill, Inc., New York, 1977. [36] Chu, S. C., and Pan, K. C., "Dynamic Response of a High-Speed Slider-Crank Mechanism with an Elastic Connecting Rod," Journal of Engineering for Industry, ASME, Series B, Vol. 97, No. 2, May 1975, pp 542-550.