UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Dynamics of deformable multibody systems Huang, Xiaodan 1998

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_1999-389014.pdf [ 6.27MB ]
Metadata
JSON: 831-1.0089233.json
JSON-LD: 831-1.0089233-ld.json
RDF/XML (Pretty): 831-1.0089233-rdf.xml
RDF/JSON: 831-1.0089233-rdf.json
Turtle: 831-1.0089233-turtle.txt
N-Triples: 831-1.0089233-rdf-ntriples.txt
Original Record: 831-1.0089233-source.json
Full Text
831-1.0089233-fulltext.txt
Citation
831-1.0089233.ris

Full Text

DYNAMICS OF DEFORMABLE MULTIBODY SYSTEMS By Xiaodan Huang B.A.Sc. Tsinghua University, Beijing, China 1984; M.A.Sc. Tsinghua University, Beijing, China 1987  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY  in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF MECHANICAL ENGINEERING  We accept this thesis as conforming to the required standard  THE UNIVERSITY OF BRITISH COLUMBIA  October 1998 © Xiaodan Huang, 1998  In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis forfinancialgain shall not be allowed without my written permission.  Department of Mechanical Engineering The University of British Columbia 2324 Main Mall Vancouver, Canada V6T 1Z4  Date:  Abstract  This thesis presents an investigation of the nonlinear dynamics of arbitrary deformable multibody systems that undergo large translation and rotation movements and small elastic deformations. The objective of this study is to develop an accurate and efficient modeling method to meet the requirements in system design and control. A general implicit formulation based on the joint coordinate method for arbitrary tree or closed-loop deformable multibody systems (MBS) is developed by defining a new topological matrix. The newly-developed formulation and code have been verified numerically by investigating the total energies and strain energies of two different conservative rigid-flexible systems. The absolute error of the total energy should be several orders smaller than the strain energy to ensure the validity of the small elastic deformations. An experiment study on the dynamic responses of a 3-D test rig with both joint and link flexibility was conducted to verify the simulations. The results of the simulation and experiment show good agreement in both the time and frequency domains. A simulation comparison amongst the joint coordinate method, Order N method and absolute coordinate method was performed. It was shown that the Order N formulation method may induce chaotic behavior in nonchaotic systems due to the propagation and enlargement of numerical errors. The commercial software ADAMS was used as a representative of the absolute coordinate method. The results demonstrate that the newlydeveloped implicit formulation has some advantages compared with other methods. Two geometrical nonlinear effects are discussed in this thesis. The simulation and experimental results of the test rig show that the formulation including the foreshortening effect overestimates the nonlinearity. The thesis also presents simulations of modal expression, nonlinear coupling effects and chaos.  ii  Table of Contents  Abstract  ii  Table of Contents  iii  List of Tables  vi  List of Figures  vii  Nomenclature  xi  Acknowledgments  xvii  1 Introduction  1  1.1 Motivation  1  1.2 Preliminary Remarks  2  1.2.1 System Description  2  1.2.2 Equations of Motion  3  1.2.3 Numerical Solution and Error  5  1.2.4 Nonlinear Systems and Chaos  6  1.3 Literature Review  8  1.3.1 Formulation Methods  9  1.3.2 Software  18  1.3.3 Validation  20  1.3.4 Geometric Nonlinearities  21 iii  1.4 Objective and Scope  22  2 New Formulations for Deformable MBS  25  2.1 Topological Description  25  2.2 Background for Description of Deformable Bodies  27  2.2.1  The Position Description by the Assumed Vibration Modes  27  2.2.2  The Position Description by Finite Elements  28  2.2.3  The Orientation of the Body Coordinate System  32  2.3 Dynamic Equations for a Single Constrained Deformable Body  33  2.4 Velocity Transformation for DMBS  43  2.5 Dynamic Equations for DMBS  52  2.6 An Alternative Method for Deriving Velocity Transformation  53  2.7 Discussions  58  3 Numerical Validation of the General Purpose Software  61  3.1 Introduction  61  3.2 Software Implementation  61  3.3 Total Energy Validation for Two Different Examples  63  3.4 Modal Representation  79  3.5 Flexibility Coupling Effects  82  3.6 Summary  83  4 Experimental Validation 4.1  85  Introduction  85  4.2 Physical Description  85  4.3 Instrumentation  86 iv  4.4 Calibration Experiments  88  4.5 Comparison Between Experiment and Simulation Results  91  4.5.1 The Modelling of the Test Rig  92  4.5.2 Experiments Performed and Comparison with Simulations  93  4.6 Comparison with the Simulation Results of ADAMS  102  4.7  107  Summary  5 Simulation Comparison of Different Formulation Methods 5.1  Introduction  109  5.2 Recursive or Order N Formulations  109  5.3  Simulation Comparison  112  5.3.1  Numerical Validation of Derived Order N Formulations  113  5.3.2  Numerical Comparison  115  5.3.3  Comparison with ADAMS  122  5.4 The Chaotic Behavior in Simulation Results due to the Formulation Method 5.5  Summary  6 Geometrical Nonlinearities 6.1  Preliminary Remarks  6.2 The Two Geometrical Nonlinear Effects  7  109  ..126 129  130 130 130  6.2.1  Nonlinear Effect Induced by Axial Forces  134  6.2.2  Foreshortening Nonlinear Effect  135  6.3 Comparison of Simulation and Experiment Results  136  6.4  143  Summary  Summary and Conclusions  144  List of Tables  3.1  The Parameters of the two links  64  3.2  The parameters of the three links  71  3.3  The parameters of the eight-link manipulator  75  3.4  The initial and final positions of the joints  75  3.5  Eigenvalues and eigenvectors of a two-link manipulator  81  5.1  The parameters of a two-rigid-link mechanism  113  vi  List of Figures  1.1  A closed-loop deformable multibody system  2  1.2  Nonchaotic behavior when y = 0.6, pi = 10.0  1.3  Chaotic behavior when y = 0.06, /I = 10.0  2.1  Topology of a Spanning Tree  25  2.2  Coordinate systems for assumed mode method  28  2.3  Coordinate systems for finite element method  29  2.4  Relative motion among bodies  43  2.5  A branch of a tree system  54  2.6  Two adjoining bodies  56  3.1  Flow chart of the developed software  62  3.2  A two-link manipulator  63  3.3  Elastic displacement along the axis direction of joint two ( Z )  3.4  Elastic displacement perpendicular to the axis direction of joint two (Y )  66  3.5  Angular displacements of both joints (upper-joint 1; lower-joint 2)  66  3.6  Total energy  67  3.7  Strain energy  67  3.8  Kinetic and potential energies  68  3.9  Total energy comparison  69  7 7  65  2  2  3.10 Strain energy comparison  70  3.11 Kinetic and potential energy comparison  70  3.12 A three-link manipulator  71  3.13 Angular displacements of joint 1 and joint 3 (upper-joint 1; lower-joint 2)  72  3.14 Translation displacement of joint 2  72  3.15 Elastic displacement of link 3 (perpendicular to the axis direction of joint 2)  73  vii  3.16 Elastic displacement of link 3 tip (parallel to the axis direction of joint 2)  73  3.17 Total energy and strain energy  74  3.18 An eight-link manipulator  75  3.19 Displacements of joint 1 ~ 4 vs time (s)  76  3.20 Displacements of joint 5 ~ 8 vs time (s)  76  3.21 Tip elastic displacement of link 1 in the y direction of the body coordinate system  77  3.22 Tip elastic displacement of link 1 in the z direction of the body coordinate system  77  3.23 Tip elastic displacement of link 2 in the x direction of the body coordinate system  78  3.24 Tip elastic displacement of link 2 in the y direction of the body coordinate system  78  3.25 Comparison of the six elastic coordinate case and the two modal coordinate case  81  3.26 Comparison between rigid-rigid and rigid-flexible manipulator (joint 1)  82  3.27 Comparison between rigid-rigid and rigid-flexible manipulator (joint 2)  83  4.1  The test rig  86  4.2  Strain gauge wiring  87  4.3  Test devices  87  4.4  The calibration of potentiometer 1  88  4.5  The calibration of potentiometer 2  89  4.6  The calibration of strain gauges  89  4.7  Damping measurement method  90  4.8  Tip acceleration of the experiment  90  4.9  Tip acceleration of the simulation  91  4.10 Joint damping and flexibility  92  4.11 Response of joint 1 at initial condition 0, = 0° and 6 = 30°  93  4.12 Response of joint 2 at initial condition 0, = 0°and 0 = 30°  94  2  2  4.13 Jointed-end strain of the second link at initial condition 0, = 0°and 0 = 30° ( Z )  94  4.14 Spectrum of strain simulation response at initial condition 0, = 0°and 0 = 30°  96  2  2  2  viii  4.15 Spectrum of strain experiment response at initial condition ©, = 0° and 0 = 30°  96  4.16 Response of joint 1 at initial condition 0, = 0°and 0 = 45°  97  4.17 Response of joint 2 at initial condition©, = 0°and © = 45°  97  2  2  2  4.18 Jointed-end strain of the flexible link at initial condition©, = 0°and © = 45°( Z )  98  4.19 Spectrum of strain simulation response at initial condition©, = 0°and 0 = 45°  98  4.20 Spectrum of strain experiment response at initial condition©, = 0°and© = 45°  99  2  2  2  2  4.21 Spectrum of strain simulation response at initial condition ©, = 0° and © = 40°  100  4.22 Response of joint 1  100  4.23 Response of joint 2  101  2  4.24 Jointed-end strain of the flexible link ( Z )  101  2  4.25 Response of joint 1 at initial condition©, = 0°and © = 30°  102  4.26 Response of joint 2 at initial condition©, =0°and © =30°  103  4.27 Response of joint 1 at initial condition©, = 0°and © = 45°  103  4.28 Response of joint 2 at initial condition©, = 0°and © = 45°  104  4.29 Strain signal comparison at initial condition ©, = 0°and © = 30°  105  4.30 Strain signal comparison at initial condition ©, = 0°and © = 45°  105  4.31 Spectrum of the strain signal by ADAMS at ©, = 0°and © = 30°  106  4.32 Spectrum of the strain signal by experiment at ©, = 0° and © = 30°  106  4.33 Spectrum of the strain signal by ADAMS at ©, = 0°and © = 45°  107  4.34 Spectrum of the strain signal by experiment at ©, = 0°and 0 = 45°  107  5.1  Angular displacements of the two joints (upper-joint 1; lower-joint 2)  114  5.2  Total energy  114  5.3  Kinetic and potential energies  115  5.4  Angular displacements of the two joints (upper-joint 1; lower-joint 2)  116  5.5  Elastic displacements of second link tip along the axis direction of joint 2  117  5.6  Elastic displacements of second link tip perpendicular to the axis of joint 2  117  2  2  2  2  2  2  2  2  2  2  ix  5.7  Total energy  118  5.8  Strain energy  118  5.9  Kinetic and potential energies  119  5.10 Total energy  121  5.11 Angular displacements of the two joints (upper-joint 1; lower-joint 2)  121  5.12 Elastic displacement of second link tip along the axis direction of joint 2  122  5.13 Angular displacement of joint 1  124  5.14 Angular displacement of joint 2  124  5.15 Tip elastic displacement along the axis direction of joint 2  125  5.16 Tip elastic displacement perpendicular to the axis direction of joint 2  125  5.17 Angular displacements of the two joints (upper-joint 1; lower-joint 2)  127  5.18 Elastic displacement of the flexible link tip  127  5.19 Total energy  128  5.20 Strain energy  128  6.1  Beam vibration under axial forces  134  6.2  Foreshortening effect  135  6.3  Comparison between nonlinear and linear under initial condition 8 = 30°  138  6.4  Comparison between nonlinear and linear under initial condition 6 = 30°  139  6.5  Comparison between nonlinear and linear under initial condition 9 - 45°  139  6.6  Comparison between nonlinear and linear under initial condition 6 = 45°  140  6.7  Comparison of experiment and  2  2  2  2  simulation including only axial forces 6 = 30°  140  2  6.8  Comparison of experiment and simulation including both nonlinearities 6 = 30°  6.9  Comparison of experiment and  2  simulation including only axial forces 6 =45°  141  2  6.10 Comparison of experiment and simulation including both nonlinearities 9 = 45° 2  x  141  142  Nomenclature cij  acceleration vector of a point on body j  {a ,a ,a Y, (b , b , b ) x  y  z  x  y  z  J  node coordinates of element j referred to the ith body coordinate system  A'  rotation transformation matrix of ith body  B, B  global velocity and acceleration transformation matrices transformation matrix represents the connectivity of element j of body i  B'  transformation matrix of body i expressing the boundary  2  conditions B  the modal matrix of the ith body  l m  [C]  damping matrix  C, D, E,U, C, D, E, U C' , C J  ,J  matrices defined in equation (2.137) transformation matrices between the y'th intermediate element coordinate system and the ith body coordinate system  C  the viscous damping matrix of the system  d  D'  spatial differential operator  <J>  E'  elastic coefficient matrix  F'  conservative force vector of the body i  F,  the ith generalized active force  F*  the Ith generalized inertia force  <J>  e  G',G',G'  transformation matrices associated with generalized orientation vector <JO  Hy, H  component block matrices of B and B  I  identity matrix  tj  xi  moment of inertia matrices of the jth unit length element of body i with respect to the element coordinate system moment of inertia matrix (of element j) of body i referred to the body coordinate system inertia coefficient matrix due to elastic deformation component velocity transformation matrices of B, defined in equation (2.164), (2.152), (2.165) the stiffness matrix (of the j element) of the ith body with respect to the global inertia frame element stiffness matrix with respect to the element coordinate system the torsional stiffness matrix of the system the nonlinear stiffness matrices due to axial forces and foreshortening length of element j of body i Lagrangian of the ith body inertia matrices of the whole system and (of the jth element) of the ith body inertia mass of a point on body j component inertia matrices of the matrix M'  <J>  component inertia matrices of the matrix M'  <J>  generalized inertia matrices generalized inertia matrices total body number and joint definition point number of the system total element number and mode coordinate number of body i xii  N'  fl' ,NQ  <J>  <J>  t  <J>  constant matrices with respect to the ith body coordinate system defined in equation (2.32) and (2.28)  , Njj, Nj.  j  constant matrices with respect to the ith body coordinate system defined in equations (2.49) - (2.51)  pi<j>pi<j>  position and velocity vectors of any point on body i  Pj  active force vector of a point on body j  PS  l p  constant matrix associated with the shape function of the ith body and the position vector of joint definition point p  q,q  l  independent generalized coordinate vectors of the total system and the ith body  q'f  independent elastic displacement or mode coordinate vector of body i  q'jj  the ith mode coordinate of body i  q'f  nodal displacement vector of body i with respect to the body frame  q'j  the y'th element nodal displacement vector of body i referred to the body coordinate system  q'j  the jth element nodal displacement vector of body i referred to the intermediate element coordinate system  q'f  generalized modal coordinates of the ith body  QQ  generalized nonconservative external force vector of body i  Q'  generalized conservative force vector of body i  Q  generalized nonlinear force vector  m  c  g  Q  s  Q  v  l  l  Q'^,Qg  generalized elastic force vector of body i generalized centrifugal and coriolis force vector of body i mass moments of element j of the ith body Xlll  generalized force vectors shape function (of element j ) of body i with respect to the body coordinate frame the 7th mode shape of body i with respect to the body frame the jth element shape function of body i referred to the element coordinate system constant matrix of element j of body i defined in equation (2.52) kinetic energy of the ith body position vector of any point (on element j) of body i with respect to the body coordinate system under deformation position vector of any point (on element j) of body i with respect to the body coordinate system under undeformation elastic displacement vector of any point (on element j) of body i referred to the body coordinate system position vector of the joint definition point o on body i with respect to the global inertia frame position vectors of the two joint definition points on one joint with respect to the global inertia frame the strain energies due to axial forces and foreshortening potential energy of body i volume (of element j) of body i potential energy of conservative forces and strain forces the hth translation or rotation axis vector on body j partial velocity elastic displacement vector of any point on element j of body i with respect to the intermediate element frame xiv  XI, X 2  nodal vectors of element j of body i referred to the body frame  Y, Y  absolute coordinate vectors of the total system and the ith body  l  y  the y'th generalized speed  a  the stiffness proportional damping constant  P  the mass proportional damping constant  d  partial differential operator  }  7  relative translation displacement of body j  1  e  strain vector  &l  rotation angle about the /ith joint rotation axis of body j  t?#  elastic rotation angle of joint definition point R of body k  A, fj.  Lagrange multiplier vectors  %ipT\ii  components of coefficient matrices defined in equation (2.133)  1<J>  and (2.134) 7iy  component of the body path matrix of the system  p  mass density (of element j) of body i  l<J>  o  stress vector  l < J >  ]jT  summation  T  translation displacement along the hth joint axis of body j  3 h  (p  generalized orientation coordinate vector of the ith body  1  (j> ,a ,y/ J  J  J  Euleranglesofbodyy  O i  constraint Jacobian matrix of body i  Xik  component of the joint path matrix of the system  l  (o , (5 l  l  angular velocity vectors of the ith body referred to the inertia and the body frames, respectively  Q.  J  relative angular velocity of body j  XV  Superimposed Symbols differentiation with respect to time of a variable the skew symmetric matrix of a vector  Right Superscripts T  the transpose of a matrix or a vector  -1  the inverse of a matrix  xvi  Acknowledgments  I greatly appreciate the support, the trust and the advice which I have received from my supervisor, Professor A. B. Dunwoody, throughout my graduate studies. I would like to thank Professors Stanley Hutton, D. P. Romilly and Dale B. Cherchas for many useful suggestions during the development of this thesis. Professor Gary Schajer also deserves special thanks for his advice on the experiment design. I would also like to thank Professor Siegfried F. Stiemer for his kind and generous help on providing the software and the computer. I would like to express my gratitude for the two-year financial support of the University of British Columbia through UBC Academic Award: University Graduate Fellowship (St. Johns Scholarship). And the thanks also go to my supervisor for his consistent financial support throughout my graduate studies. My final but foremost thanks go to my husband, Jifang, my lovely son, Ximai, and my parents whose understanding and support made this work possible. To all of them, I dedicate this thesis.  xvii  Bibliography  148  Appendices  158  A  The Orientation Transformation Matrices  158  B  The Invariant Matrices in Finite Element Method  159  C  The Nonlinear Stiffness Matrix Due to Axial Forces  165  D  The Nonlinear Stiffness Matrix Due to Foreshortening  xviii  166  Chapter 1 Introduction 1.1 Motivation In recent years, deformable multibody systems (DMBS) have been intensively studied as a result of growing needs for the design of high speed, lightweight, precision systems. Many mechanical and structural systems in the world, such as space structures, robots, vehicles, mechanisms and aircraft which undergo large translation and rotation displacements, are made massive in order to increase rigidity or are driven slowly so that dynamic flexibility is not significant. As a result, more power is needed to drive them and lower work efficiency is achieved. An example is the 15-meter long space shuttle remote manipulator system (SRMS) of NASA which is used in the assembly of space platforms and of large communication systems. It can only move slowly due to its low natural frequency [1]. The use of lightweight materials would reduce the driving power and increase the response speed. However, the lighter members are more flexible. From the design point of view, it is necessary to be able to accurately evaluate the elastic deformations due to large and fast rotational and translational motions of a multibody system. Simulation is an important tool in mechanical design and in understanding the dynamic behavior of deformable MBS, which are represented by a number of rigid or flexible bodies connected by ideal joints and force elements. The time and money required to evaluate the dynamic responses by numerical modeling are orders of magnitude less than would be required for physical testing. However, some physical tests are still necessary to verify the simulation results to ensure that the numerical modeling is accurate and correct. Research into the dynamics of deformable MBS plays a central role in both control and simulation. In order to help engineers to design better products and design efficiently, many 1  Chapter 1. Introduction  2  researchers continue to search for better ways to describe the dynamic behavior of deformable MBS in an accurate, efficient and simple form. The work presented in this thesis is a contribution to this effort.  1.2 Preliminary Remarks Any solution scheme for obtaining the highly nonlinear dynamic response of a deformable MBS with complicated topology must incorporate three key procedures; describing the system, deriving the equations of motion and solving the equations.  1.2.1 System Description Many multibody systems have tree-like structures. If the system graph has closed loops, as shown in Figure 1.1, a tree structure is made by cutting a joint in each independent closed loop. The resulting structure is called a spanning tree.  Figure 1.1 A closed-loop deformable multibody system  A methodology should be defined to describe the system topology of a spanning tree in terms of how the bodies and joints are connected with each other. Meanwhile, a set of generalized coordinates is used to express the motion of the rigid MBS. The generetlized  Chapter 1. Introduction  3  coordinates can be either absolute coordinates which represent absolute positions and velocities of each body in the system, or relative coordinates (or joint coordinates) which express positions and velocities of degrees-of-freedom of each joint in the system. The number of absolute coordinates is larger than the number of relative coordinates because relative coordinates are independent and represent degrees of freedom of a rigid body system. These two sets of generalized coordinates will induce different formulation methods. The same pair of approaches are used for deformable MBS with the addition of extra degrees of freedom to describe the deformations of the individual bodies. Those deformations can be described by the finite element method using nodal displacements and element shape functions or by the assumed vibration modes method using vibration modes and modal coordinates.  1.2.2 Equations of Motion The derivation of equations of motion for a flexible MBS is very complicated due to the Wghly-nonlinear coupling. Although the dynamic principles used in developing MBS formulations are not new, the resulting equations have different forms depending on the system description method and the processing method. The point deserved to be noted is the solutions may be different although the same numerical method and dynamic principle are applied. The formulations can be distinguished by their explicit or implicit form. An explicit method creates equations of motion by explicitly solving some variables numerically and then substituting these solved variables to get equations for other variables. Consider a nonlinear system being represented by the equations as follows:  M (q)q\ u  M (q)q, 2l  + M (q)q  = Q,(q)  (1.1)  +M (q)q  =  (1.2)  n  22  2  2  Q (q) 2  Chapter 1. Introduction  whereM (q),  M (#)and M {q)  M (q),  u  4  21  n  are generalized inertial matrices. Q (q)  22  x  #1 Q {q) are generalized force vectors. And q =  and  is a generalized displacement vector.  2  112. The explicit equations can be derived as:  ^ =[M ( )-M (^)M- (^)M ( )]" [a(9)-M (^)M - e (^)] 1  1  11  9  =  12  1  2  2I  9  1  12  2  2  2  (1.3)  M-\q)[QM)-M (q)M- lQ {q)} n  2  2  q = M~ (q)[Q (q) - Af (q)q\) 2  2  2  (1.4)  2I  The implicit equations are developed as:  M (q)  M (q)  M {q)  M (q)_  u  2x  n  22  <7i A .  (1.5) M4)_  The explicit and implicit equations are different. The explicit equations are heavily dependent on numerical calculation. The matrices M"'(^)and M (q) 22  in the explicit equations have to  be calculated numerically. Their numerical errors could be enlarged and propagated in equation (1.3) since M (q) x2  and Q {q) 2  are the multiplication factors which could be very  large. The situation would deteriorate if there are several variables that depended on such recursive relationships or if one of those inverse matrices is ill-conditioned. In an explicit method, an assumption already has been made automatically in the formulation development. That assumption is that the explicitly-solved variables in the equation development are only related to the variables at the previous time step. In contrast, an implicit method builds up equations of motion exactly without intermediate variables being numerically solved. In other words, explicit formulations are approximate due to the dependence on numerical calculation in forming the equations. One question raised here is how much accuracy these equations represent and whether a correct solution can be obtained no matter what kinds of systems are being dealt with.  Chapter 1. Introduction  5  A formulation for flexible MBS should have the following capabilities: (1) Description of the dynamics of the individual bodies and of the interconnections between them should be straightforward. (2) The approach should be easily extendible from tree to closed-loop systems. (3) The results of the formulation should have high accuracy and stability independent of the system being modeled. (4) The computational efficiency should be high and computer coding easily implemented. (5) In the case of real time simulation for large scale systems, the formulation should be suitable for parallel computation.  1.2.3 Numerical Solution and Error The numerical solution of a set of differential equations can use either an explicit or implicit integration solver. Explicit solvers, such as Runge-Kutta, generally require fewer evaluations of the differential equations per time step, but they are conditionally stable and therefore put limitations on the time step that can be used. For nonlinear problems, ^stabilities are harder to detect [2]. Implicit solvers, such as backward difference formulas (BDF) provide stable solutions independent of the step size, although they are more computationally-intensive per time step than explicit solvers. The time step can be much larger in implicit formulations than in the explicit formulations [2]. The selection of the particular time-history integration solver to be used is dependent on whether equations of motion are expressed explicitly or implicitly and the requirements of speed, accuracy and stability. The explicit formulation method in MBS requires many inversions of matrices which are based on manipulation of results one by one. The numerical error of one matrix inversion is propagated and enlarged and may cause other matrices to become ill-conditioned or at least inaccurate. The major cause of ill-conditioning is the significant difference in the elements of the coefficient matrix. But ill-conditioning may also arise even when the physical system is  Chapter 1. Introduction  6  stable because of the way the computer operates on the numbers [2]. The uncertain numerical errors (truncation error and rounding error) in the manipulations of numbers in inverse matrices may induce unpredictable responses.  1.2.4 Nonlinear Systems and Chaos The irregular and unpredictable time evolution of many nonlinear systems is termed chaos. It occurs in mechanical oscillators such as pendula or vibrating objects as well as in other fields such as chemistry, celestial mechanics and electrical circuits [29-31]. Whenever dynamical chaos is found, it is accompanied by nonlinearity. The effect of a nonlinear term often results in a periodic solution unstable for certain parameter choices. Its central characteristic is that the system does not repeat its past behavior. The unique character of chaotic dynamics may be seen most clearly by imagining the system to be started twice, butfromslightly different initial conditions. This small initial difference can be thought as resulting from measurement error or computation error. For nonchaotic systems this uncertainty leads only to an error in prediction that grows linearly with time. For chaotic systems, on the other hand, the error grows exponentially in time, so that the state of the system is essentially unknown after a very short time. The possibility of chaotic motion would exist when thefirst-ordergoverning equations of a system have at least three variables and a nonlinear term [30]. For example, a forced vibration system might be modelled as [30]:  dy  dy  2  ~dr  +  Y^:  +  ,  3  (y - 3 0 = Msinf  ,  N  (1.6)  where the constant coefficients y,p: represent the dissipative effect and forced vibration magnitude, respectively. The three variables are  Chapter 1. Introduction  7  1*2= — 2  dt  And there is a nonlinear term y . Whether the motion is regular or chaotic depends on the 3  choice of the parameters y, fi . When y = 0.6,  = 10.0, the responses y and y both have  periodic motions, shown in Figure 1.2. The chaotic motions, shown in Figure 1.3, will happen when y = 0.06, \i = 10.0.  Time (s)  Figure 1.2 Nonchaotic behavior when y = 0.6, \i = 10.0  Time (s)  Figure 1.3 Chaotic behavior when y = 0.06, fi = 10.0  Chapter 1. Introduction  8  There is a high risk of chaotic motion in a deformable MBS since there are many variables and nonlinear terms in the equations of motion. Whether the motion is chaotic or nonchaotic depends on the values of the parameters of the systems. When the numerical errors are uncertainly enlarged due to improper formulation methods, the chaotic motion might be simulated even if the parameters of the system are in the range of nonchaotic motion. The examples given in Chapter Five will give a better demonstration for this.  1.3 Literature Review In recent years, the issue of deformable MBS has been developing in a tremendous way. The derivation of equations of motion for flexible MBS has been presented in a variety of forms. The mathematical models of such systems have been formulated using generalized NewtonEuler equations, Kane's equation, Lagrange's equations and variation principles amongst others [3-4,11-12,15-18,20,24,62]. Mainly four formulation methods have been proposed to describe the complicated nonlinear coupling between the small elastic deformations and the large rotation and translation displacements. They are the Cartesian coordinate method [35,24], relative coordinate method [16-19], recursive or order n method [20,25-28], and joint coordinate method [36-38]. There are different methods to model the flexible bodies, such as lumped masses and springs, finite elements, assumed vibration modes, Rayleigh-Ritz and component mode synthesis [3,16-17,49,61]. The methods of lumped masses and springs and finite elements require a large amount of computation time. Using assumed vibration modes effectively reduces the distributed parameter system to a discrete system. However, the choice of vibration modes is not easy. Whalen [46] showed in an experiment that the vibration frequencies and modes change while the flexible body undergoes large rotation motion. Though the topology of deformable MBS can be very complicated with tree or closed-loop structures, many formulations only deal with some special systems such as chain structures or only consider single degree-of-freedom joints. There are only a few formulations in published  Chapter 1. Introduction  9  papers which deal with spanning tree systems having different numbers of degrees of freedom of joints. Some popular programs for rigid MBS in recent years have been extended to include flexible bodies such as ADAMS and MADYMO. The influence of the small elastic deformations of a body to the large overall motion of the entire MBS is approximated or neglected [63]. Other software such as DISCOS and TREETOPS consider the nonlinear coupling  effect. In published papers, few articles demonstrate the validation, either  experimentally or numerically, of the software or formulations, especially for the small elastic deformations. Control experiments of the elastic bodies have been developed in a variety of methods [74-76]. But their dynamic models are not explicitly validated. The dynamic equations of the experimental models are mainly derived manually instead of using the formulations of MBS. Moreover, the dynamics of flexible MBS with variable kinematic structures has been studied for the cases where the constraints imposed on a MBS may change in the operating range and result in higher frequency vibration of a flexible structure [67-71]. Nonlinear elastic deformations of flexible structures (geometrical nonlinear) have drawn more and more attention due to the design requirement of high speed and heavy load [52-56]. Although all of these issues have been widely studied, disagreements still exist. The following detailed review, grouped by formulation methods, software, validation methods and geometrical nonlinear, explains why further research is necessary.  1.3.1 Formulation Methods There are many forms of formulations for deformable MBS. They can generally be divided into four groups, depending on the coordinates used for describing the motions of the system and deriving the equations of motion.  Chapter 1. Introduction  10  1) Absolute Coordinate Method The most straightforward approach is to formulate the equations of motion in terms of Cartesian coordinates and elastic coordinates [3-5]. In this formulation, the algebraic constraint equations that describe the mechanical joints are adjoined to the system differential equations of motion using the vector of Lagrange multipliers. This leads to a mixed system of differential and algebraic equations (DAE) requiring little effort to formulate. Some commercial software using this approach are ADAMS and DADS. The formulation can be obtained as:  MY +  <l = Qe + Q + c  Q =Q  (1.8)  v  <D(MO = 0 where M is an inertia matrix composed of a number of  block matrices in diagonal  connection where each block represents the inertial matrix of one body. Q ,Q ,Q e  c  v  are the  external applied force vector, conservative (gravity and elasticity) force vector and quadratic velocity (coriolis and centrifugal) force vector, respectively. Y is a displacement vector that includes position, rotation angle and elastic coordinates of each body. O is a joint constraint vector and X is a Lagrange multiplier vector. The equations of motion can be constructed in a systematic way, easily extended to closedloop systems and amenable to parallel computation. However, this approach leads to a large number of generalized coordinates, constraint equations and differential equations of motion. Therefore, the numerical computation is not efficient, although sparse techniques in matrix manipulation can be used. Also the relative coordinates are not readily available from the Cartesian coordinates since control variables are often  relative coordinates (joint  displacements or velocities) so that it is difficult to model coupled control and mechanical systems. Another disadvantage of this approach is that the resulting DAE can not be  Chapter I. Introduction  11  formulated properly as state space equations as are required by numerical solvers. Usually the second derivatives of the algebraic (constraint) equations rather than the algebraic equations themselves are used to construct the final state space equations. The procedure can be described as following: Differentiation of the joint constraints yields:  Y=A dY  where A =  dY dY  Y-2^—Y dYdt  Substitute (1.9) into (1.8):  (1.10)  dt  2  I  0  0  M  (1.9)  0 dY,  (1.11)  K  0 where  \Z = Y  dY  0 (1.12)  The solution of such a set of state space equations often drifts away from its constraints [5,6]. Many techniques and methods have been developed to tackle this problem [5-8] , for instance Baumgarte's stabilization method [8,9], the mass-orthogonal projection method [5,10-12] and the generalized coordinate partitioning method [5,13,14]. But the accuracy or computational efficiency of these methods are still under investigated. The method of Pradhan et al. [24] essentially belongs to this category although it is a variation on the absolute coordinate method. The absolute coordinates are transferred into a set of generalized coordinates through a velocity transformation. However, the number of the generalized coordinates is the same as the number of absolute coordinates. Lagrange multipliers are incorporated in the dynamic equations for uncut and cut joints.  Chapter 1. Introduction  12  2) Relative Coordinate Method The Relative Coordinate Method, was first proposed by Thomas R. Kane [15] based on Kane's equations and was used in the dynamics of spacecraft as rigid MBS. This approach defined some scalars as generalized speeds. The angular and translational velocities of each body could be described by generalized speeds. Some vectors called partial angular and translational velocities were vector functions of the generalized coordinates and time, but not the generalized speeds. The generalized forces were obtained as projections of forces and moments onto the partial angular and translational velocity vectors. The final formulations were state space equations which were represented by generalized speeds. The typical commercial software is TREETOPS. The procedure can be described as following: The translational and angular velocities of body i can be expressed as:  where n is the number of DOF of the system, y is a generalized speed and vf is a partial i  velocity. Kane's equation:  F,+F*-=0  l = l,~-,n  F^IJvf  (1-14)  dPj  (1.15)  7=1  F^lll^i-ajdntj)]  (1.16)  7=1  where N is the number of bodies, P , m , a are the active force vector, inertial mass and b  ;  7  7  acceleration vector of a point on body j, respectively. F and F* are the generalized active t  and inertial forces. This approach leads to the smallest and most strongly coupled equations of motion. The advantage of using relative coordinates is that only independent variables are included in the  Chapter 1. Introduction  13  final formulations of tree systems so that the equations can be solved efficiently. However, the derivation of the equations of motion is very complicated so that the procedure is very difficult  to extend for closed-loop problems except if the constraint equations are kinetic  rather than kinematic [18]. Unfortunately, joint constraints are usually kinematic. Also it is difficult to formulate forcing functions [34] and the method is not well suited to parallel computation. The method later was extended to deformable MBS [16-19]. One technique used in this extension is to treat the flexible bodies as a number of rigid segments and springs [16,17] so that the formulations used in rigid MBS can be directly used without modification. The stiffnesses of the springs are evaluated by the material and the geometry of the deformable body. The other technique is to develop the equations from the beginning using shape functions to describe deformations [18,19]. Only chains with one DOF joints [19] or tree structure systems [18] have been considered due to the complicated coupling. Although Singh et al. [18] have developed kinetic equations for tree-structure systems, the calculation formulations for partial velocities are contradictory because the partial velocities in the beginning were defined as vectors but in the end they became scalars (formulation (26) and (28)). Also in the formulations of the partial velocities, only the deformation of one joint definition point on each body is considered. But the deformation effect of the other joint definition point is not included (equation (28)). These two joint definition points of a joint are the points that are placed on two adjacent bodies and idealized to represent the relative motion of the joint. The situation of only considering the deformation of one joint definition point is a special case, that is one of the boundary conditions of the deformable body is clamped.  Moreover, the final dynamic equations in [18] can only be solved by implicit  numerical integration since the generalized inertia and the generalized speed are coupled and can not be explicitly separated. The final kinetic equations can not easily to be used since the  Chapter 1. Introduction  14  equations are expressed in vectors and dyadics. There were no simulation examples or validation of the formulations in [18].  3) Recursive or Order N Method Another popular method is called the Order N (or recursive) method. Commercial software using this method include DISCOS and SIMPACK. Computational efficiency in MBS simulation is an important issue especially for large scale systems. The Order N method was proposed based on the requirement of miiumizing computations. Initially, work on the Order N method dealt with rigid robot manipulators [20,21]. For an N-link chain configuration with simple revolute or prismatic joints, the total computational complexity of the joint space inertial matrix and its inverse matrix is 0(N ) 3  according to traditional algorithms. Therefore,  the computation cost grows rapidly with N. The Order N method can formulate equations of motion in a way that the number of calculations per integration step increases only linearly with the number of bodies (or degrees of freedom of the system). Instead of inverting an N*N joint space inertial matrix, the Order N method inverts operational space inertial matrices which are always 6*6 for rigid MBS [22]. The comparison of the operation numbers for different recursive algorithms can be found in [22,23]. Based on this significant advantage in operations per integration step, many researchers have been trying to develop Order N (or recursive) formulations for deformable MBS [20, 25-28]. The main idea in the derivation of the Order N method for deformable MBS is to express the kinematic relations of each body one by one from the base body to each tip body and then explicitly solve joint variables and elastic variables by inverting the operational space inertial matrices one by one from each tip body to the base body. The formulations can be obtained as follows [20]: The absolute velocity vector v, of body i is related to the absolute velocity vector v,^ of body i-1 and the relative vector v- of joint i by  v = *i v ._ +r v; J  >w  I  1  ll  (i.i8)  Chapter 1. Introduction  15  Differentiation of the above equation yields  v = /?,,_, v,_, + T|, v' + b ;  (1-19)  t  where b- = /?,,_, v + t vf M  (1.20)  a  The dynamic equations of the system can be obtained as:  Tj,R Mv = T£ R f + TpR f T  where f,f  T  T  (1.21)  r  represent generalized force vectors due to external forces, the quadratic velocity  r  forces and constraint forces at the cut joints. Equation (1.20) also can be written as where f  T? £ Rj. M  jVj  = T„ f* + T« £ R]jJ  =R f  (1.22) (1.23)  T  Starting with the final body, i = n, equations (1.19) and (1.22) are  V„ = /?„.„_, V -, + TnnK + K  (1.24)  B  TlR M v T  nn  nn  nn  n  = T V + T Rl f T  n  n  n  (1.25)  T  nnJ n  nn  nn  n n  nnJ n  ^  s  Substitute (1.24) into (1.25) and solve for v' as a function of v _,, then substitute this result n  n  into (1.24) to obtain an equation for v„ as a function of v _,. Thus equation (1.22) is reduced n  and will not include v . n  The operations are repeated next for the case i = n — \. Lastly, the results are as follows:  *,'=4>,-.+fl;+w/  >  (i.26)  v,-=A,--iVM+«/ + ^ - *  (1.27)  r  where A _ = -N'„M' R f  it  x  iM  K = M'r TZ l  (1.28) (1.29)  M; = TJM;T  (  M^M.+^RlM-A,  (1.31)  U  I. O) 3  Chapter 1. Introduction  16  Though the equations of motion have been formulated using Newton-Euler equations, Lagrange's equations or variation principles, the Order N (or recursive) method is an explicit formulation method which builds up formulations with the help of numerical calculations. The numbers of numerical inversions of inertial matrices (M\~ ) may cause some of the matrices x  to become ill-conditioned which will significantly affect the accuracy of the solutions [2]. From a nonlinear dynamics point of view, chaotic behavior may be induced due to the illconditioning of some matrices [29-31]. Thus, this method may break down in the case of stiff systems, for instance, when dealing with flexible bodies and with control or contact problems [32]. Even if it produces a solution for a stiff system, the solution may not be reliable. Moreover, the explicit method is time step size dependent. Not only accuracy but also stability of the solution is influenced by the time step. The stability conditions can be detected for linear equations using explicit integration solvers. However, it is very difficult to predict stability conditions for such numerical calculation dependent formulations with large nonlinearities. Van Woerkom and Boer [25] have pointed out that the selection of time stepsize is of considerable importance. Although a comparison of operation counts for different formulation methods in rigid MBS shows that the Order N method has fewer operations per integration step, the calculation time in a time period may not be reduced since the time step must be very small for numerical stability. Also the extent of algorithms being applied in parallel computations is limited by the topology of the systems due to the recursive relationships of the bodies.  4) Joint Coordinate Method A variation on the relative coordinate approach is called the joint coordinate method. The equations of motion for each body are first described in terms of Cartesian coordinates and Lagrange multipliers. The equations are then transformed to represent joint (or relative) coordinates using a global velocity transformation proposed by Jerkovsky [33] and developed  Chapter 1. Introduction  17  by Kim [34] and Pankiewicz [35]. The formulations in rigid MBS can be obtained as following: Using Cartesian coordinates, build up the dynamic equations for each body: T X = Ql + Ql + Ql  M'Y' +  (1.32)  The dynamic equations of the system (as in the absolute coordinate method) can be written as: MY +  & = Qe+Qc+Q  (1.33)  v  \dYj  Define a global velocity transformation: Differentiate equation (1.34) to yield  Y = Bq  (1.34)  Y= Bq + Bq  (1.35)  Multiply equation (1.33) by B and substitute equation (1.35) into equation (1.33), then T  B MBq T  + (—B) X  = B (Q  T  +Q +Q-  T  e  For tree-configuration systems,  C  V  (1.36)  MBq)  B= 0  Therefore equation (1.36) becomes  B MBq = B (Q T  T  e  (1.37)  + Q + Q - MBq) c  v  For closed-loop systems, the constraint equations of cut joints  Y) = 0  Therefore equation (1.36) becomes  B MBq T  + (—B) X T  = B (Q T  e  +Q +QC  V  MBq)  (1.38)  The global velocity transformation matrix B and its derivative matrix B in a rigid MBS can be developed by two methods proposed by Kim [34] and Pankiewicz [35].  Chapter 1. Introduction  18  This formulation has some of the advantages of both the relative coordinate and absolute coordinate approaches. It is easy to formulate force elements and constraints for cut joints. The resulting differential equations for tree-structure systems are expressed in terms; of a minimum number of variables. Kim demonstrated that this kind of formulation has very high computational efficiency [34]. Moreover, the structure of this formulation can be easily implemented in parallel computers and for closed loop systems. However, the necessary equations have only been developed for a rigid MBS. The coupling between large translation or rotation motions and small elastic motions in a deformable MBS makes the global velocity transformation very difficult to derive. Nikravesh and Ambrosio [36,37] have applied the joint coordinate method to formulate equations of rigid-flexible MBS. The flexible bodies were assumed to have rigid parts that flexible parts can be attached to. The analysis of a partially deformable body is initially approached by treating the rigid and flexible parts as separate bodies. Then the parts are connected by noting that points on the boundary of rigid and flexible parts have the same global displacements. Thus the velocity transformation formulations still can be used as in a rigid MBS. Although Pereira and Proenca [38] have said that it is possible to define a linear transformation between the vector of system generalized velocities and the time derivative of the vector of system relative coordinates, they never showed how it can be obtained. Instead, they only developed recursive velocity and acceleration relations. Thus, the most difficult problem in this approach for deformable MBS is still unsolved.  1.3.2 Software General purpose MBS software has been developed in different ways and applied in different areas. Some programs are capable of efficiently analyzing large scale complex mechanical systems. For instance, ADAMS (Mechanical Dynamics, Michigan), DADS (University of Iowa), MADYMO (Netherlands) and SIMPACK (Germany) have become popular [39,40] for  Chapter I. Introduction  19  dealing with rigid MBS. ADAMS and DADS use the absolute coordinate method to formulate equations of motion. They perform coordinate reduction by generalized coordinate partitioning which is not very efficient in computation. DADS applies four Euler parameters rather than three Euler angles as in ADAMS to describe the orientation of each body. MADYMO and SIMPACK both use the Order N formulation method. Although there have been some popular commercial programs capable of being applied to deformable MBS in recent years, the equations of motion are basically not from recently derived formulations but some assumptions are made in them. For example, a combination of explicit FEM commercial software such as PAM-CRASH or DYNA3D with rigid MBS commercial software such as MADYMO has been developed by the TNO Crash-Safety Research Centre to analyze vehicle occupant safety [41,42]. One method is to discretize flexible bodies as numbers of element nodes by commercial FEA software, then these nodes can be considered as a number of rigid bodies connected by force elements. The stiffness and damping matrices of these force elements are calculated by FEA software. ADAMS uses this method. Another method is to obtain the solution of a rigid MBS first assuming the forces coming from flexible bodies are known from the previous time step. Then the forces at the next time step are solved using explicit FEM software [41]. Not only are these approaches approximate but also the explicit FEM puts limitations on the time step that can be used. The typical time step is in the microsecond level [41]. Thus the computation efficiency is not high. Ho and Herber [43] described another method for dealing with small deformation. The solution of a rigid MBS is obtained first, when the distributed structural flexibilities are taken into consideration, the positions and orientations of the bodies must be perturbed slightly. This method is used in the program ALLFLEX. The programs DISCOS (Order n method) [64,65] and TREETOPS (relative coordinate method) [18] are mainly designed for deformable MBS. They both consider the nonlinear coupling between small deformation and the large rigid body motion.  Chapter 1. Introduction  20  1.3.3 Validation Although many formulations for deformable MBS have been developed and published recently, very few of them are verified to be correct either by experiments or by theoretical calculation. Most of them only provide formulations and algorithms. Some of them do give simulation examples to test the formulations [13,19,44] but the correctness of the simulation cannot be judged because no comparison and verification are available. Validation on the basis of total energy balance has been pursued by Caron [28] , Grewal [61] and Woerkom [25]. The judgment of correctness is based on whether the ratio of energy absolute error to total energy is small or not. The question raised here is whether this standard can be used for the indication of the validity of deformable MBS formulations or not. Generally, the strain energy in deformable bodies is very small compared to the other components of the total energy. In order to verify the small deformation as well as large translation and rotation displacement, the absolute error of total energy compared with strain energy rather than the ratio of energy absolute error to total energy could give a convincing indication. If the maximum strain energy is several orders of magnitude larger than the level of the maximum absolute error of v the total energy, then the small deformable displacements can be trusted. Otherwise, the strain energy is no greater than the noise of numerical calculation and the results cannot be trusted. Validation on the basis of experiments appears to be rather scarce in the published literature. Modal analysis has been used to identify experimentally vibration modes and frequencies at different angular positions to compare with calculation results [45,46,47]. Also, experiments based on flexible link control have been studied extensively [46-49]. The experimental methods using strain gauges have been introduced to measure the dynamic response of flexible links [50,51]. The dynamic modelling of an experimental rig in 2-D control studies mainly is usually developed manually. The developed formulations for a deformable MBS rarely is verified by experiments in published papers.  Thompson and Sorge [62,73] showed  experimental verification in 2-D. Nevertheless, the simulation results of derived formulations  Chapter I. Introduction  21  for deformable 3-D MBS have not yet been found to be compared with experiments in the published papers to date. Chun, Turner and Frisch [66] designed an experimental validation of the DISCOS software to be carried out at NASA Goddard Space Flight Center. A seven-link robot arm was to be tested. The links were rigid bodies. However, the use of an harmonic drive introduces nonlinear stiffness into the joint dynamics. Instead of modeling the harmonic driveflexibilityas a nonlinear spring between two rigid bodies, they modelled the driven link and the harmonic drive as one flexible body with one elastic mode [66]. A comparison of the experiment results and simulation results has not been found in recently published papers.  1.3.4 Geometric Nonlinearities The problem of geometric nonlinearity, another interesting topic in deformable MBS, has been drawing more and more attention recently due to large load and high speed inertia. Different methods have been proposed to tackle this problem accurately and efficiently [5256,72]. Kane et al. [72] studied the behavior of a cantilever beam built into a rigid body that is performing a specified motion of rotation and translation. The effect of the transverse displacement on the axial displacement is incorporated in the kinematic description of the deformation. The stretch in the beam as well as the transverse displacements are considered as a set of generalized coordinates. Hanagud and Sarkar [52] did not treat the stretch as a generalized coordinate, instead they discretized the axial displacement. The nonlinear straindisplacement relations were then used to derive the equations. Some papers [54-56] use the nonlinear strain-displacement relations to derive a general formulation for a general deformable MBS. A survey and simulation comparison on this topic can be found in [53]. Sharf [53] presents a precise understanding of the existing methods and how they relate to each other through an in-depth review of some publications. The simulation results using different methods show that there are large differences amongst them. However, they have not been compared with any kind of experimental results. No detailed analysis has been made to  Chapter 1. Introduction  22  explain the differences in simulation results. Shabana et al. [54] pointed out that the longitudinal displacement could occur due to two effects; one due to axial forces and the other due to the foreshortening effect as the transverse displacement causes an axial displacement. Most of the published papers [52-56] obtained strain energy with the nonlinear terms of thirdorder and fourth-order by retaining the nonlinear strain-displacement relation. This results in the controversy over which term can be neglected [54]. The disagreement on which formulation is accurate and efficient remains an open problem. The method using the nonlinear strain-displacement relation could lead to a misunderstanding of the geometrical nonlinear behavior, that is, geometrical nonlinear behavior is caused only by large deformations. However, in actuality geometric nonlinearity arise when deformations are large enough to significantly alter the way load is applied or the way load is resisted by the structure [2]. Large deformation is one important condition. The other important factor is the applied load or the structural constraint. Recognizing this point gives us a better understanding why geometrical nonlinear behavior occurs and which can be neglected for the third-order term or fourth-order term of strain energy.  1.4 Objective and Scope Based on the above literature review, the objective of this thesis is to pursue an efficient and accurate simulation solution for an arbitrary deformable MBS to meet the simulation requirements in system design and control. This goal will be approached by exploiting the potential of dynamic formalism and numerical computation. Its implementation will be verified experimentally and theoretically. The specific steps to obtain the objective are described as follows: 1) To define a new topology description for an arbitrary deformable MBS to obtain a rigorous mathematic model; then using the joint coordinate method and the new topology description to develop a general implicit formulation which is accurate, efficient, simple  Chapter 1. Introduction  23  (straightforward) , easy to use in closed-loop systems and with high potential for parallel computation. An important point is that the new formulation be stable and able to deal with any kind of stiff system. 2) To implement the implicit formulation as a general purpose program and verify it by numerical calculation. In this part, the verification standard is the combination of absolute error of total energy and strain energy. The simulation of different configurations of systems provides evidence that the software can deal with arbitrary deformable MBS. 3) To design a test rig and verify the software by comparing experiment results and simulation results. The large joint motions and the small elastic motion are measured and calculated. The joint motions simulated by ADAMS software are also compared with the experimental results. 4) To develop an explicit (or Order N ) formulation and compare it with the developed implicit (joint coordinate method) formulation. Also the simulation results of the joint coordinate method and the absolute coordinate method (through ADAMS commercial software) are compared. 5) To investigate geometric nonlinear behavior by separating the two nonlinear effects caused by axial forces and foreshortening. To compare the simulation results of the test rig with the experiment results to demonstrate the effects of geometric nonlinear on a real system. This thesis is composed of seven chapters. In Chapter 2, a detailed formulation development of the joint coordinate method for deformable MBS is presented. A new topology definition is introduced to meet the requirements of the joint coordinate method. Lagrange equations are applied to derive dynamic equations for each body independently. Then two methods are used to develop the global velocity transformation matrix and its derivative matrix. A treeconfiguration system is presented to show the structure of the transformation matrix. The advantages of this formulation are discussed in the end of this chapter.  Chapter 1. Introduction  24  Chapter 3 and Chapter 4 present the implementation and the validation of the formulation. Chapter 3 describes the capacities and the structure of the software. A two-link 3-D manipulator with one rigid link and one flexible link is used as a test example. The validation also includes using different numbers of elements for the flexible link. In order to show that the software is able to deal with different configuration systems, a three-link manipulator and an eight-link manipulator are tested. Also modal representation and the effect of small elastic deformations on the largerigidmotion are investigated. In Chapter 4, a test rig is described to verify the software experimentally. An overall description of the experimental method and devices is presented. The prelirninary experiments which include calibration experiments of the sensors and the parameter identification experiments are discussed. The test rig is tested at two different initial conditions and the simulation results and experiment results are compared. Also the simulation results from ADAMS software are compared with the experiment results. Chapter 5 involves the formulation comparison. First, an explicit formulation (Order N ) is developed and implemented in software. The software is verified numerically using a rigid MBS. Then the simulation results of three approaches, joint coordinate method, Order N method and absolute coordinate method (ADAMS), are compared. The objective is to demonstrate the advantages of the joint coordinate method. In Chapter 6, two geometric nonlinear effects are investigated and compared with the experiment results of the test rig. The purpose is to illustrate that geometric nonlinear modelling is not a purely mathematical issue. It depends on which geometric nonlinear effects dominate in real applications. Finally, the conclusions and further research work are presented in Chapter 7.  Chapter 2 New Formulations for Deformable MBS This chapter presents a detailed development of a general implicit formulation based on the joint coordinate method for deformable MBS. A new topological definition is introduced to describe a deformable system topology. Two methods are used to derive the global velocity transformation matrix and its derivative matrix. The derived new formulation has some advantages which are very important for the dynamic simulation of a complicated deformable multibody system.  2.1 Topological Description Consider a deformable system of N bodies that are interconnected by joints with 1 ~ 6 h  DOF as indicated on Figure 2.1.  Figure 2.1 Topology of a Spanning Tree  If the system graph has closed loops, a tree structure is made by cutting a joint in each independent loop. The result structure becomes a spanning tree. A method has been reported by Kim and Vanderploeg [34] that selects cuts which minimize the number of generalized coordinates and differential and algebraic equations. A base body is defined as the body fixed to the ground for a grounded system or any of the bodies for a floating system. Let the bodies  25  Chapter 2. New Formulations For DMBS  26  of the spanning tree system be labeled from 1 to N , starting at the base body or "root" and b  ended up at tip bodies or "leaves". Then the body path matrix can be defined as [34] 1 **  =  if body j is between the base body and the ith body  |0  (2.1)  otherwise  Usually joints, simplified as two joint definition points connected by ideal joint axes, are also labeled in this way [24]. However, this description can not satisfy the requirement of the joint coordinate method, which needs to distinguish between the two joint definition points due to the possible deformation of one or both joint definition points. Thus, instead of labeling joints, joint definition points are labeled from 1 to JV. in the same numbering order (from root to leaves) and connection directions are described from lower numbered points to higher numbered points. This is the new definition (contribution) made for the requirement of joint coordinate method. Therefore, the joint path matrix is defined as  1  if the join t definition poin t k is directly on the path from the base body to the ith body and connects the higher numbered body  Xik  ~  -1  if the join t definition poin t k is directly on the path from the base body(2.2) to the ith body and connects the lower numbered body  0  otherwise  For example, the body path matrix and joint path matrix for the system in Figure 2.1 are "10  0 0 0 0'  "1 0  0  0  0  0  0  0  0  0"  1 1 0 0 0 0  1 -1 0  0  0  0  0  0  0  0  1 1 1 0 0 0 1 1 1 1 0 0  1 -1 1 -1  0  0  0  0  0  0  1 -1  0  0  0  0 0  1 1 0 0  1 0  1 1 0 0  1 1  x =1  -1 1 -1  1 -1 0  0  0  0  1 -1  0  1 -1 0  0  0  0  1 -1  1 -1  The body path matrix and joint path matrix make the system topology unique.  Chapter 2. New Formulations For DMBS  27  2.2 Background for the Description of Deformable Bodies A deformable body could be discreted by the Rayleigh-Ritz approach, finite elements, component mode synthesis, lumped rigid segments and springs or assumed vibration modes [61]. However, the Rayleigh-Ritz approach and component mode synthesis are not commonly used in deformable multibody systems due to the difficulty of choosing admissible functions. The lumped rigid segment and spring method divides a deformable body into a number of rigid segments and springs according to the geometric parameters of the body. Then the formulation for rigid multibody systems, which is developed by Kane's equations, is applied. A detailed description is presented by Huston and Wang [17]. This section only deals with the methods which treat the deformable bodies as continuum structures, i.e. finite elements and assumed vibration modes. The discretization via finite elements makes the solution to the approximate problem converge to that of the real problem as the mesh size is reduced. But, for a complex problem, this results in many elastic degrees of freedom. The large number of degrees of freedom reduces the simulation speed. The assumed vibration mode method can reduce the number of elastic degrees of freedom significantly. However, the admissible functions ( vibration modes) satisfy only the geometric boundary conditions but do not, in general, satisfy the natural boundary (force) conditions. The directions of the forces at boundaries are constantly changing for the deformable bodies undertaking large rotational motions. Whalen et al. [46] have done an experiment to investigate the changes in mode shapes under different angular displacements. Usually, the mode shape functions under no rigid motion condition are chosen for the sake of convenience. This certainly results in some approximations.  2.2.1 The Position Description by the Assumed Vibration Modes Assume the elastic deformation of the deformable body / can be described by N  m  modal  coordinates q' and mode shape functions S . Two coordinate system frames are applied to 1  f  Chapter 2. New Formulations For DMBS  28  describe the rigid and elastic modes. As shown in Figure 2.2, XYZ is the global inertial frame and X'Y'Z'  is the ith body coordinate system which is rigidly attached to a point on the  undeformed body. This frame is a floating frame which represents large translation and rotation motions.  Figure 2.2 Coordinate systems for assumed mode method The elastic displacement of any point u^ on body i referred to the X'Y'Z'  is:  (2.3)  ";=^;=I% N  M  7=1  where S] is the 7th mode shape of body i (a function of u^) and q' is the 7th mode f]  coordinate of body i. The position vector P' of a point UQ on body i can be obtained as P' = r' + A'u' in which  (2.4)  u' = UQ + u' = UQ + S'q' f  f  (2.5)  where r' is the position vector of the origin of the ith body coordinate system and A' is the transformation matrix that describes the rotation of the ith body coordinate system with respect to the global inertiaframe,u' is the position of any point  on body i referred to the  ith body coordinate system under the deformation condition.  2.2.2 The Position Description by Finite Elements The position description using finite elements is similar to assumed modes. However, the development is more complicated since the shape functions usually are given with reference to  Chapter 2. New Formulations For DMBS  29  the element coordinate system rather than referred to the body coordinate systems as in assumed modes. Moreover, the element assembly and the imposition of boundary conditions need to be considered in the finite element method. Thus, a unique shape function with respect to the body coordinate system needs to be developed first. Assume a deformable body can be divided into N  e  elements. Four kinds of coordinate  system frames are used to describe joint and elastic motions for deformable bodies [3]. As shown in Figure 2.3, X Y Z and X'Y'Z'  are the global inertia frame and the ith body  coordinate system, respectively, as in the assumed mode method. In addition, intermediate element coordinate systems X^V'Z''  are rigidly attached to the body coordinate systems to  represent the initial orientation of the ith element with respect to the body coordinate system. Element coordinates X' Y Z' J  ij  are attached to a point on the y'th element of body i .  j  Figure 2.3 Coordinate systems for finite element method  The elastic displacement of any point u^ on element j of body i referred to X'Y'Z'  can be  represented by element nodal displacements through element shape functions. Therefore, the elastic displacement vector W' of a point on element j of body i with respect to X' Y' Z' J  W  S  J  T 2  S*J  qj  J  J  is:  (2.6)  Chapter 2. New Formulations For DMBS  30  where S' and q'j are the jth element shape function and the node displacement vector of J  element j referred to  XYZ. iJ  iJ  iJ  The elastic displacement vector TL'} of any point on element j of body i with respect to X'Y'Z'  is then described as:  uj = C W ij  (2.7)  ij  where C' is a transformation matrix from the intermediate element coordinate system J  XYZ ij  ij  to the body coordinate system  ij  X'Y'Z'.  For beam elements, the transformation matrix can be evaluated as [3]:  C> =  2  (2.8)  + l c  2  b'i-a'i  where  -'J  +c?  Vl-a*  y  _ JTl""J:  —  , c 'j  Vt-a* z  I'  z  1  (2.9)  (2.10)  (a ,a ,a )'  J  x  y  z  and (b ,b ,b )'  J  x  y  z  are the coordinates of the nodes of beam element j with  reference to the body coordinate system  Therefore  C" 0 0 C  u  X'Y'Z'.  0  0  0  0  0  0  C  0  0  0  0  c  IJ  u  1)  (2.11)  Chapter 2. New Formulations For DMBS  31  where q'j is the node elastic displacement vector of element j with respect to the body coordinate system  X'Y'Z . 1  The ith body coordinate system X'Y'Z' represents a unique standard for all elements on this body. Let q' be the total vector of nodal displacement vectors of body i resulting from the f  finite element discretization. The relationship between the nodal displacement vector q' of f  body i and the element nodal displacement vector q'j can be written as  q) = B?q'  (2.12)  where B[ is a constant Boolean transformation whose elements are either zeros or ones and j  serves to express the connectivity of this finite element. In addition, a set of reference conditions that are consistent with the kinematic constraints on the boundaries of the deformable body are imposed to represent unique deformation fields with respect to X'Y'Z'.  Let q' be the vector of independent elastic displacements of body i f  and B' be a linear transformation which is also a constant Boolean transformation and serves 2  to express reference conditions. Then  q' =BWf  (2-13)  f  Therefore the elastic displacement vector of any point on element j of body i with respect to X'Y'Z'  is then written as  u'j = CS'C^B^  = S"q'  f  The position vector P' of a point 11^ on element j of body i can be obtained as J  (2.14)  Chapter 2. New Formulations For DMBS  32  P = r' + Au ij  (2.15)  ij  u = u* + uj = u$ + s' q'f  in which  j  ij  (2.16)  where r is the position vector of the origin of the ith body coordinate system and A' is the l  transformation matrix that describes the rotation of ith body coordinate system with respect to the global inertial frame. u' is the position of any point ujj on the y'th element of body i J  with reference to the ith body coordinate system under the deformation condition. Equation (2.4) and (2.5) are virtually identical to equation (2.15) and (2.16), so either approach to defining deformed shapes may be applied during the remainder of the development.  2.2.3 The Orientation of the Body Coordinate System The orientation of the body coordinate system can be represented by a generalized coordinate vector (p. If Euler angles are used, <p =  the angular velocity vector of  the ith body with respect to XYZ and x V z ' is  co = G y = A'W • 1  1  fl>'' = G V  where co'and co are the angular velocity vectors referred to XYZ and X'Y'Z' l  A',G',G'  (2.17)  (2.18)  respectively.  are transformation matrices associated with the generalized orientation vector (p  [57] and are given in Appendix A.  Chapter 2. New Formulations For DMBS  33  2.3 Dynamic Equations for a Single Constrained Deformable Body The dynamic equations of motion for a single constrained deformable body can be formed using Lagrange's equations. A advantage of using this principle is that it easily allows a total energy balance check. Define Y' as the absolute coordinate vector of the ith body coordinate system augmented with the vector of generalized elastic coordinates,  T =[r ,<p ,q\l l  l  (2.19)  For a deformable body in a MBS, the body usually is connected to other bodies by joints. Those joints limit the motion of the body so that constraint forces are applied to the body, and the generalized coordinates Y' are not independent of each other. Therefore, the Lagrange equations cannot be represented only by partial differential equations but can be described as [3] d_ dlt_  dlJ_  dt dY  dY  1  (2.20)  1  in which the Lagrangian is  L =V-V L  (2.21)  where T is the kinetic energy of the ith body; V" is the potential energy including both strain energy and potential energy due to any conservative external forces; O'^A represents the constraint forces expressed by Lagrange multipliers A and the constraint Jacobian matrix  <f>' i; Y  QQ represents generalized nonconservative external forces. The following equations present in detail the development of the equations of motion for a single constrained deformable body. The formulations are suitable to both assumed modes and finite element discretization. The only difference between the two methods is that there is no assembly in the assumed mode method. To distinguish them, the terms in angle brackets are excluded for the assumed mode method.  Chapter 2. New Formulations For DMBS  34  First, consider the kinetic energy of body i.  Kinetic energy: V =< X >  T > < i >  = ~ < X > J*  7=1  ^  Differentiate (2.4) or (2.15):  P  j  =  1  l  = r + Au  ,<J>  l  = r'+A'(W  Xu  (2.22)  p^P^P^dv ^  , < J >  = r' -A'u G'(p'  +  l<J>  )  +  A ti l  i<j>  A'S q' ,<J>  f  +A'S' q'  ,<]>  <]>  >  = [l - A^ G'  AS  i<}>  i<j>  f  ] q>'  (2.23)  if where the sign ~ over a vector denotes a skew symmetric matrix. / is an identity matrix. Substitute (2.23) into (2.22), get 1  Ne  T =-Y {< l  l>M )Y  lT  l<J>  7=1  2  I where  M  i<j>  =j  (2.24)  l  -A'u'^G  p  iT  i<j>T  m  r(p  (pq>  m  symm.  -G u  i<j>  symm.  rr  i<j>  G u u G'  i<j>  m  AS  1  iT  r  rrirf (pf  m  >  r  S  i<J>T  f  j  i<j>  dvKJ>  >  \KJ> (2.25)  Chapter 2. New Formulations For DMBS  35  m 0 in which  m^ = j" p ldv ;>  l<i>  i<j>  = J p A'S  <J>  dv  i<j>  i<j>  '<j>  m  = Jp  where  iT  =-G  m  ,<;>  s  m 0  0  0 m  i<i>T  , T  / < > r  , < v >  )G ' = -A N; G •  i<j>  !  i<j>  i<J>  (2.27)  <J>  =  i<J>  A'N  (2.28)  i<j>  )G' = GH'^G  IT  Jp W S  5  (2.26)  'dv** = G (J p'Ou'^n'^dv'^  I<I>  I<J>  ><J>  = A ' j " p S dv  i<j>  " C " = J P' G u u G <J>  =0  = - A ' (J p T dv  mrip  m  i<!>  0  dv  <i>T  dv'  /V,'^ = J p^u* * 1  i<j>  = Jp  < ; >  dv  l<1>  = G ' J p'^U^S  i<}>  i<7>  =Jp  dv  T  (s;  i<J>  <i>T  s(  <J>  +s!  (M^ +5  i<j>  sf  <i>T L  , < ; >  = G l£  i<J>  iT  (2.30)  j>  (2.31)  + sys; )dv  J>  <J>  ; )dv  i<j>  ?  (2.29)  1  = N'  <j>  i<j>  +N  q'  i<J>  f  (2.32)  KJ>  ji<]>  w  A2  A3  I  7  22  symm.  (2.33)  23  I  33  KJ>  +2[j <> (u^sr+u^rs^ ;  P  ^ > 4f  +qf [J p'<> ( 5 / ^ 5 / ^ + S ^ S y )dv  i<j>  ]q>  (2.34)  Chapter 2. New Formulations For DMBS  (k,l,m= 1,2,3 and  36  k*l*m)  i;: >=- j p ^ u r u ^ d v ^ = - j  p^u^u^dv*  j  -j  (u^>s: > u^>sr>yv < %  i<j> P  j  i  j  +  -\q'f[lp (s; s: +sir i<j>  <j>T  j>  sr^dv^  7  (2.35)  (/,m = 1,2,3and l^m)  i«r  =\p {u^s: -u^sr]dv-^ i<i>  i>  (k,l,m = 1 — > 2— > 3 and k^l^m,  which means they follow the rotatory rule such  that if k - 1, I = 2 and m = 3;if k = 2, 1 = 3 and m = 1; if k = 3, 1 = 1 and m = 2)  where p'  <J>  ,v'  and w' are, respectively, the mass density, the volume and the mass of  <j>  <J>  body / or of element j of body /. And 7io = [« <J>  For beam elements, let XI'-' = [a, a a f 2  [«oi "02 "  Ql=\  0  r = k *2 0  3  p r}da ij  ij  3  3  0]  u  u ]' ,S  =\S  b bf  and define  i<!>  02  and X2' = y  03  T  2  3  X  F + ' ' ^ i? g f = X l * ' + z V [ § »? d  '  7  S  2  S] ' . 3  (2-38)  (2.39)  =i p gda  (2.40)  /<>' = J pV-rfda*  (2.41)  /j' = J p g da  (2.42)  ij  ij  ij  2  ij  Chapter 2. New Formulations For DMBS  37  (2.43) X V =l>j ^ , r)' • -, J  7  J  in which  £ *  J  J  Q' *  J=L  J =  V  L  ~  J  V  then N% = Jp^'dv' = — [X1+ X2] 7  N = jp°S dv ij  l,  } p( o u  + P  = C ($p S dv )c B?B! iJ  u  iJ  u  v  =  u  l  C'N'jC'B^B^  (2.45)  (aJtm+/a mC(iU) + 2a / (c (Jt,2)G, + C(/c,3)<2)) k=p,q 2  Oq)  u  (2.44)  j  dv  ik  w a  (2.46)  ,  Jk  J  ?  +  I l C{k,l) -mC(A;,l) + C(/c,2)^ + C(*,3)fi k=p,q\  +  I / (c (*,3)/,, +C U,2)/ +2C(*,2)c{ik,3)/^ )f k=p,q  3  s  3  2  2  ?  ^ + v ( « p C ( ^ l ) + a c(p,l))+-^c(p,l)  (2.47)  <?  + / (a C(^2)+a C(p,2))(2 +/ (fl C( ,3)+^C(p,3))e 2  2  p  9  n  p  9  s  -(c(/ a)c( ,2) + C(p,2)c(^))^+^(c(p,l)c(^3) + C(/ ,3)c(^l))!2 J  ?  J  g  + / (c(p,2)c(<7,2)/ +C(p,3)c( ,3)/ +C(p,2)c(^,3)/ + C(/ ,3)c(< ,2)/ )]'' 3  g  9  n  TC  J  ?  7?g  (/>,<? = 1,2,3, and p*q)  let /V' = \cjpSdvCB.BJ 7  = [ c ^ C 5 , B f = [/V," /Y /Y ]' 7  2  r  2  r  3  F  (2.48)  Chapter 2. New Formulations For DMBS  N« = [c{\lpt;Sdv)CB B^  38  = [CIN^CB B ]  x  X  = [#J /7[ / V j f  2  2  (2.49)  ^ = [ c ( J / p r ^ v ) c 5 5 [ = [ c / ^ C B B f = [JvJ ^ A ^ J  (2.50)  ^ = [ c ( J / p ^ v ) c 5 5 ] = [ c / A / ^ C f i B f = [/7J N  (2.51)  ?  1  2  1  2  y  1  T  2  1  NN% = [BlB C SS CB B ] T  then [\pu S dvJ 0p  p  [\ S S dv] T  P  x  2  q  i2  A ^ f  [B BfC (\pS S dv)CB B T  T  2  + C{pX)N^ + C(p,2)N„  q  = f^NN  J  p  pq  = aN  q  =  T  x  2  T  p  q  x  2  + C(p,3)N  qq  C{p,k)C{q,r)  (2.52)  (2.53)  (2.54)  kr  k=\ r=l  (p,q = 1,2,3)  where the above invariant matrices Nf, N'^ , N^„, N'^ and SS  pq  (p,q = 1,2,3) and the 3-  D beam shape function S will be given in Appendix B [3]. iJ  Therefore, Lagrange's equations can be rewritten as:  d_  = MY + l  dt  Define  l  Q' =-M Y +2 l  lif'Y'--I  dY  (2.55)  (2.56)  l  v  which are the centrifugal and coriolis( gyroscopic) force components. They are written as  i  \ y~.iT y~.iT  Q' =[Q" v  r  Q  V<P  y~.iT  (2.57)  Chapter 2. New Formulations For DMBS  where Q  l  = -[m r  + m (p + m ^ / )  l  vr  l  rr  r(p  = &'A N a'+A N G q>'-a' ,  39  ,  ,  ,  t  A' N q' + A'b'G' <p  ,  l  t  l  f  b'=N'q'  (2.58)  (2.59)  f  ^;'=<X>^v;  (2.60)  <  7=1  Ne  N = < £ > AT l  (2.61)  7=1  Qv<p =-[G l G iJ L  +G  L  yf " V =  G  i _ l  P  d I  k=\  I^G^y-G' 1'^q'f  LT  W  (2.62)  7  (2.63)  2  +  qn  kk(-i\  2  dqj  ,  |  d I  /,m=l l-±m  lm-i-i  (2.64)  dqj-  (2.65) 7=1  Ne  V  =< Y > 7  i<J>  (2.66)  7=1  ^ = 2 1  " A +2 < £ > j p (k,l,m= 1,2,3 and  , < j >  _ _ (5/ ' s/ < 7  k*l*m)  > 7 ,  < j >  _ _ +S* S* > T  J >  yv  ,v  (2.67)  Chapter 2. New Formulations For DMBS  dl lm dq'  <1 > Jp  40  far s i r + ^ > s r ytv**  i<J>  f  < X >J^ { s i ^ s ^  +  s^ s; )dv j>T  <j>  i<j  if  (2.68)  (l,m= 1,2,3 and /*m)  Potential Energy:  V' = V + V l  (2.69)  l  c  s  where V and V ' are the potential energy due to conservative forces and strain forces. They l  c  5  are obtained respectively as following:  a) Potential energy due to conservative forces Assume F is a conservative force vector of the ith body which is applied at position P l  i<j  c  The generalized force vector Q is developed by virtual work as l  c  dw[ =  SP  i<J>  F^ 8P T  i<j>  (2.70)  = 8r + SA'u l  i<J>  = Sr - u'^G'&p' 1  + AS  Sq'  i<J>  + AS  f  ;<>  (2.71)  <Vf  Substitute equation (2.71) into equation (2.70), to get 5w = [F ' l  c  - F- U G'  T  T  C  F- AS ] b\p  i<j>  T  i<j>  l  (2.72)  Sq'  f  = [F;  where  u'  <J>  - F?U  T  = Au  ,<J>  G  I<I>  1  F; AS T  I<J>  ]  (2.73)  (2.74)  Chapter 2. New Formulations For DMBS  41  Therefore the potential energy due to conservative forces is  r  v: =  -W: Y T  I  -[K  A'S ]  -  T  i<j>  (2.75)  b) Strain energy The virtual work due to internal elastic forces can be written as [3] Ne  Ne  Swl =< X > <5<  =-  J>  <^>ja  7=1  r'<7> _  Se dv  i<J>T  i<j>  i<J>  (2.76)  7=1  j?i<j>gi<]>  where a' , e' <J>  <J>  (2.77)  and E'  <J>  are, respectively, the stress, strain vectors and the matrix of  elastic coefficients. The strain displacement relation can be written as  «j>  £  =  «J>U'<J>  D  where D'  =  D'^S'^q'f  (2.78)  is a spatial differential operator. In the case of small strains and rotations, the  <J>  differential operator reduces to i< j>T dx  D'  <J>  =•  o 0  dy dz  2  L  o 4 40  ay dx dz 0 2 ^ 0 A dz dx dy  (2.79)  Chapter 2. New Formulations For DMBS  42  Substitute (2.77) and (2.78) into (2.76),  Ne  Ne  '  •'• W = - l l < X > J { ' D  <i>S U i >  ) E(D S i<J>  )dv 6q' = -qf < £ > JC^Sq)  i<j>  i<j>  f  J=I  (2.80)  j=i  For the finite element method, let  Ne  K = l  f  K i = B -, l  f  1  7=1  1  Ne  Ne  (2.81)  ^  where £J? is an element stiffness matrix. Its formulation for beam elements will be given in Appendix B.  Thus  bw^-q'JK^Sq)  =-Y K'SY iT  1  =-[r  iT  0 0  0  ' 8r '  <p qf] 0 0  0  Sep'  iT  1  (2.82)  0 0 4.  The generalized elastic force vector is  r.iT  0 0 0 (Tl 0 0 0 <P If \ 0 0 4 .  [ iT  iT  Qs =-[  r  The strain energy  (2.83)  0 0  0  vj = fif Y = [ iT iT ;Tl 0 0 -V <P Qf J 0 0  0  l  (' r <f>'  1  (2.84)  i  3f.  Therefore the Lagrange's equation (2.20) for one single deformable body with constraints becomes  Mf  + &I l=Q +Q +Q  l  where  i  Q =[Q'  i  l  0  T 0 r  i  v  i  c  Q$ \  T  f  +Q  i  s  Q  (2.85)  (2.86)  Chapter 2. New Formulations For DMBS  43  The generalized friction forces or driving forces Q are formed in the same way that the l  0  generalized conservative forces were. As a result, the dynamic equations for a constrained deformable MBS with N bodies can be b  written as MY + 0ln  = Q +Q +Q +Q v  c  s  (2.87)  o  where M is a block diagonal inertial matrix comprising the assembly of each body's inertia matrix; <& H represents the constraint forces amongst bodies. T  Y  y = [y" Y 17  2T  ••• Y " f N  T  (2.88) (2.89)  2.4 Velocity Transformation for DMBS This section is a contribution. Consider a deformable MBS with a spanning tree structure, a part of it can be shown in Figure 2.4. On the path from the /th body to the ith body, the joint definition points are noted as T, S, R, Q, P, O. Attach joint coordinate systems to each joint definition point noted as X^Y^Z^ in which m represents the joint definition points  A  mm m  —r  (m = 1,2,• • •,Nj) and n represents the bodies (n = \,2,~-,N  b  )  Figure 2.4 Relative motion among bodies  Chapter 2. New Formulations For  44  DMBS  a) Angular Velocity From Figure 2.4, the angular velocity of the ith body coordinate system can be written as:  co' =m +0) + Q - col) j  p  where Q  op  X' Y Z L  0  L  0  (2.90)  J  op  is joint relative angular velocity which describes the angular velocity of frame  with respect to frame X YJZ . J  0  co' and co' are the relative angular velocities of the  J  P  p  P  points P & O with reference to X'Y'Z  1  0  & X'Y'Z',  they can be expressed in terms of the  elastic coordinates of body j and i as follows [11]:  = A PS' q'  f  (2.91)  G)' = A'PStfj  (2.92)  0J  J  J  p  P  0  in which PSQ or PS I  is a constant matrix associated with the shape functions and the  position vectors of the joint definition points referred to the ith or the ;th body coordinate systems. The matrix PS will be given in Appendix B. Repeating the procedure of equation (2.90) for body j, k,l,---, and the base body,  CO=co +C0 +CL -co j  k  K  J  R  QR  Q  co =co +a>j+Q. -cos k  l  ST  (2.93) (2.94)  Substitute equation (2.93) and (2.94) into (2.90) and repeating this procedure to the base body, yields: co^^^+tx^PS^j]  (2-95)  Chapter 2. New Formulations For DMBS  where K  and % are  45  body and joint path matrices, CI is defined as the relative angular  t n e  J  velocity between bodies which includes the base body, then  if body j is a floating base body  co  J  h  (2.96)  otherwise  where V/and 6{ are the /ith joint rotation axis vector of body j and the rotation angle about this axis. Let VJ be unit vectors along the rotation axes of body j that are rigidly attached to the origin of the intermediate frame X Y„Z n  n  m  (X Y%Z for body j in Figure 2.4) and defined in this k  m  k  R  R  frame. Then  (2.97) where  is the transformation matrix that defines the orientation of the intermediate joint  coordinate system X Y Z k  with respect to the coordinate system X Y Z  k  R  R  k  R  k  k  of body it". This  transformation matrix can be written in terms of the elastic rotations.  For small rotation [27],  0  -#  A = / + tf  0  R  3  -02  3  ^1  #  2  -tfj 0  (2.98)  .  where I is a 3 x 3 identity matrix. And  K = PSkf  (2.99)  Redefine Q as J  h  (2.100)  Chapter 2. New Formulations For DMBS  46  in which, if body j is a floating base body, let  J ei = p,  d = y/ J  e{=a\  j  6  V: = G i, V> = G j, Vi = G'k s  where (j> ,a ,y/ J  J  j  are Euler angles of body j;  j  ( 2  i,j,k  "  1 0 1 )  are unit vectors representing the  orientation of the global inertial frame XYZ.  ••  fl»'=X^[XW+5i^P5^] j=l  h  (2-102)  k=]  Equation (2.17) can be written as  0'=GV  (2.103)  where & is a matrix associated with cp and will be given in Appendix A. l  Substitute (2.102) into (2.103), then  9'=i *uG'EeiVj!+%z A PSiq ] J  d  7=1  J  lt  (2-104)  f  jt=l  h  b) Translation Velocity From Figure 2.4, the position vector of the ith body coordinate frame can be written as:  r' = r +UT+TS-US+U l  + RQ-U Q + U P + PO-U J  R  J  0  (2.105)  Therefore, the translation velocity of the ith body coordinate frame can be represented as:  Chapter 2. New Formulations For DMBS  r =(r + TS+ RQ+ PO) + {u  47  f  l  l  + u -ii k  T  J  R  + u -it^) J  Q  P  Where  (2.106)  (2.107)  u = (O x u + A Ssq k  k  k  s  k  k  (2.108)  s  Repeating the procedure of equations (2.107) & (2.108) for the joint definition points T,R,Q,P,0, and substituting equation (2.90),(2.93),(2.94),(2.95),(2.108), etc. into equation (2.105), and applying the body & joint path matrices, yields:  Asi-^K^XiX^psi l=j  where  J+-  k=l  ifXik=-  11 =  (2.110)  & Xik = -1  J  (2.109)  k=l  y represents the relative translational motion of joint definition points which is defined as: j  if j is a floating base body  r  ^- h^h  otherwise  T  (2.111)  where V/ andx[ are the /ith joint translation axis vector and relative displacement along the axis.  if j is a floating base body  r  Therefore, y' =  £^i^h  +  'h^h )  T  otherwise  (2.112)  Chapter 2. New Formulations For DMBS  Redefine y  as  j  48  y =I(f^V^+T^/)  (2.113)  ;  h  Where if body j is a floating base body, let  r/ = { = r / = 0  (2.114)  x  where x ,y ,z J  J  are the absolute velocity components of the floating base body with respect  j  to the global inertial frame. Again, let VJ be unit vectors along the translational axes of the joints that are rigidly attached to the origin of the intermediate frame X Y^Z K  and defined in thisframe,then according to  K  R  R  equation (2.97),  Vi=A A V ! k  R  where  A A V^ k  k  R  + A A VJ  k  R  = AW A V k  k  k  R  Then according to [39],  R  = A W (A A A V  }  k  h  k  k  k  k  R  A W (A )'  k  k  R  (2.116)  j h  = a>  1  R  equation (2.116) becomes  Therefore  (2.115)  k  l  (2.117)  R  A A*V = CO * x V k  7  A  J h  Vl = (co + co ) x vj k  k  R  (2.118)  (2.119)  Substituting equation (2.93) into the above equation, yields  V* = -yi (co +co - &) j  J  K  (2.120)  Chapter 2. New Formulations For DMBS  49  where K is the number of the joint definition points on body j which makes % = -1. Then jK  C0 = J  K  from (2.102), get  A PSlqj  (2,121)  oo  (2.122)  j  Substituting equation (2.113),(2.120),(2.121),(2.122) into equation (2.109), equation (2.109) can be written as:  7=1  Ih  (2,123)  h  N:  where  ^ = 1 ^ - 1 ^ ^ i=j  h  k=l  ^=lLX*^Sl-^K&xX k=\  (2.124)  Xik"k+^hyi h  )  + ^W  l=l\  (2,125)  h  m=l  h  c) Velocity Transformation In general, it is difficult to develop dynamic equations in terms of independent generalized coordinates for deformable multibody systems. Dynamic equations in terms of Cartesian coordinates can easily be obtained. The joint coordinate method transforms Cartesian coordinates to independent generalized coordinates based on a global velocity transformation. The vector of Cartesian velocities of the system with N bodies is written as b  Y =  [Y ,Y ,-,Y » ] 1T  2T  N  T  T  (2.126)  Chapter 2. New Formulations For DMBS  50  A vector of independent generalized coordinates of the ith body, q', is defined by the joint relative coordinates and elastic coordinates. That is: q =[T ,0 ,q ITJlT  (2.127)  f  The vector of independent generalized coordinates for the system is written as  • IT -IT -b q ,q ,---,q ° N  T  17" (2.128)  Therefore the global velocity transformation matrix B expresses the relationship between Y and q. That is Y = B(q)q  (2.129)  The time derivative of equation (2.129) yields an acceleration transformation equation.  Y = B(q)q +B(q,q)q  (2.130)  The matrices B and B are derived from equation (2.104) and (2.123). For example, a 6-body system is shown in Figure 2.1, the structures of matrix B and B are expressed as  2l  H  3l  H  A\  H  «51 -^61  0  0  0  0  0  22  0  0  0  0  H  0  0  0  42  #44  0  0  "55  0  65  "66-  H  32  H  52  0  0  62  0  0  H  H  H  0  0  0  0  0  0  0  0  0  0  0  2\  H  22  0  3\  "32  "33  A\  "42  "43  0  "51  "52  "44 0  0  0  "55  0  -"61  "62  0  0  "65  "66-  H  H  H  Chapter 2. New Formulations For DMBS  51  In which  0  G'Vt-G'V*  G'XA,,D/  0  0  R  , where I R = / when i = j {R = 0 when i ^ j  (2.131)  (G V/;+CV/;)...(G V ;...G V ; ,  ,  ,  /  /  k=\  0  where  =  % u[  Kg  - £n C  Xik  (2.133)  n  i n  n=7+1  (2.134)  %=£^[ *-fe'- ") /- ^] £  c  D  c  (2.135) 1=j  k=l  n= j+\  N •  iij = 1  x [Et-(tm-c )D -(s -c )b(-CJD -C D n  ik  k=\  J  n  k  ul  J  j  K  J K  (2.136)  l\ is defined in equation (2.110). U=u h  U=u h  D-  APS  D = APS  E -  AS  E = AS  (2.137)  Chapter 2. New Formulations For DMBS  52  The superscripts of matrices U,C,D,E  (or U,C,D,E)  represent body numbers, the  subscripts of matrices U,D,E (or U, t>, E) represent the joint definition points. Note that the formulation can be simply applied to rigid multibody systems by eliminating matrices D and E and their derivative matrices. For higher computational efficiency, the formulation can be further developed as the following recursive algorithm:  (2.138)  (Nj l u' +c  or  (2.139)  +l  Xik  k  U=i  2.5 Dynamic Equations for DMBS For the whole system, the dynamic equations can be written as equation (2.84). Substitution of equation (2.124) and (2.125) into equation (2.84), yields  B MBq + (® B) p: T  T  Y  = B (Q + Q + Q + Q - MBq) T  v  C  S  0  (2.140)  For the constraint equations: SHK.t,Y) = 0  (2.141)  ® SY = 0 =*® B5q = 0 Y  (2,142)  Y  For a tree configuration system, q is an independent vector. Thus Ovfl = 0  (2.143)  The dynamic equation for a spanning tree system is:  B MBq = B (Q T  T  V  +Q +Q +Q C  S  0  MBq)  (2.144)  Chapter 2. New Formulations For DMBS  53  Equation (2.144) is a pure differential equation. It can be easily solved by any of a number of different numerical methods. For a closed loop system, the constraint equations for cut joints can be written as  ®*(t,Y) = 0  (2.145)  d®* =& bY = & BSq = 0 Y  (2.146)  Y  The dynamic equation for a closed-loop system is:  B MBq + (& B) fl T  T  Y  = B (Q T  V  +Q +Q + Q - MBq) C  S  0  (2.147)  o*(r,r) = o  Equation (2.147) is a differential-algebraic equation. For closed-loop systems, the dynamic equations will always be differential-algebraic equations as we found in the absolute coordinate method. Unlike the absolute coordinate method, only algebraic equations for each cut joint are needed to be constructed. The reduced number of algebraic equations will certainly lead to higher computational efficiency.  2.6 An Alternative Method for Deriving the Velocity Transformation The velocity transformation can be obtained in an alternative way which is described in paper [35] for rigid MBS. This section extends and completes the method for deformable MBS. The part of the topology description of the spanning tree system is still expressed by the body path matrix and bodies are labeled from "root to leaf. Each joint is labeled in the same way so that it agrees with the body number. In addition, the position vectors of the two joint definition points of each joint with reference to the body coordinate system are noted as u (j) and « o - 0 ) in which j is the number of the body or the joint. u (j) and u _(j) are 0+  Q+  0  Chapter 2. New Formulations For DMBS  54  represented whether the two points are connected with high or low numbered body respectively. A branch of a tree system is shown in Figure 2.5.  Figure 2.5 A branch of a tree system  ,'_L'T  Assume q =[T  niT _ « T l „ _ j , , i \JT .AT JT\  d  r  q' j and Y' = [r  <p" q'j ]  f  T  are the relative and absolute vectors  of the ith body. The absolute velocity of the base body can be written as:  (2.148)  01  For the second body, the absolute velocity can be expressed as the combination of the absolute velocity of the first body and the relative velocity between the first body and the second body, i.e. Y = J,Y +J q 2  ]  (2.149)  2  n  Substituting equation (2.148) into equation (2.149), yields  •I J\J\Qq  •2  r  I  +J{2<i =1^1^10  - [ • M i o ^12]  •2 -  20%2  J  \2\  J  •2  (2.150)  Chapter 2. New Formulations For DMBS  55  In the same way, the absolute velocities of the rest of the bodies are  Y — [J J  20 ^23]  2  with  JnO  [ ^ n - l *^(n-l)0  <?02 •3  (2.151)  ^30^03  —  J(n-l)n  (2,152)  ]  Therefore, the velocity transformation matrix B for a tree-configured system is obtained by assembling J  (n = l,---,N ) according to the body path matrix. For example, for the system  n0  b  shown in Figure 2.1, the B matrix is  J  I0  20  J  30  J  B=  40  J  0 0 00 0 00 00 0 50  . ^60  a n  (2.153)  oo7 o oo7  /50  where ^50 = [^50 ^50]  0 0 0 0  60  ^ ^60 [- 60 -^o] ^ =  partitioned into two parts which  /  are  associated with the DOF of the first two bodies and the rest of the bodies. Differentiation of equation (2.152) yields  n0  J  -  n-\J(n-l)0  J  +  n-\J(n-l)0  J  (2.154)  (n-\)n  J  Thus the derivative of matrix B also can be obtained by assembling j  n0  body path matrix in the same way that matrix B was.  according to the  Chapter 2. New Formulations For DMBS  The matrix J  N 0  56  is composed of three matrices shown in equation (2.152). The first one J _N  represents the motion due to the inboard directly-connected body; The second one („_i)o 7  represents the motion due to the other inboard bodies from the base body to the inboard directly-connected body. This matrix is obtained recursively. The third one J( -\) N  N  reflects  the motion due to the joint between two connected bodies. They can be developed as follows. Consider the two connected bodies shown in Figure 2.6. Body i  Figure 2.6 Two adjoining bodies The position vector of the ith body coordinate system is  r'=r  + A (ul  J  + S^qj) - A''(4. + S _q ) + T  J  (2.155)  l  0  f  Differentiation of the above equation and substitution into equation (2.113) yields  r' =r +U M)co -Zw(i)fi>' + A'S^qj - A'S _q J  i  0  0  + f  T  TO  < 2  156  )  h  Also, the angular velocity can be written as  = 0) +Q' + A'PS^f J  ~ A'PS^q'f  Substitute equation (2.157) into (2.156), then  (2.157)  Chapter 2. New Formulations For DMBS  57  r' =r + (U _ - U )CO + « _Q' + (u^A'PS^ J  + A'S^)qj  J  0  0+  0  +2 W +XX T  -(U^APSi+ASi)^ According to equation ( 2 . 1 2 0 ) ,  (2.158)  V = -V (co' + COQ_ - o!) l  l  (2.159)  h  h  Substitution of equation ( 2 . 1 5 7 ) into (2.159) yields  V^-V^CoU^PS^qj)  (2.160)  Substitute (2.17), ( 2 . 1 0 0 ) , ( 2 . 1 0 3 ) , ( 2 . 1 6 0 ) into ( 2 . 1 5 7 ) and ( 2 . 1 5 8 ) , then  r =r + uo.-Uto.-I.rM l  J  p q> +u _'Lvl ef\+XVf\z J  J  ,  0  h  J  i  h  h  h  f  + "o- -  A PS  q>' = G (G ( i  J  0+  j  i P  q  J  0+  f  - ( M _ A ' P V + A'S _)<7/  (2.161)  APS^q))  (2.162)  0  0  J  l>  and  + AS  1  +£ ^  + A PS^qj j  -  h  Therefore  f  G  j  "o- -  r  ¥  h  A'PS^  "o- -  )  h  )  A PS,  &G  0 0  GXr-GVi  r  J  9  if  _ AS  -u^A'PS^  +  o+  0  0  J  J  J  j  + A S^  n  -G APS l  (2.163)  n  i)  Chapter 2. New Formulations For DMBS  Let/ = n, j = n-l  (n =  N ), then b  u,  w_ 0  58  G  A ft!  «-i  f  A - PS n  J  7( n - l ) n  y  fe  M  A^PS,  (2.164)  04-  0  0- lr K  l  J  h  0  i!  + A - S^ n  0+  G"G"  v  l  U  0 - Ir V  0  GYJ-GX"  0  0  -U _A"PS _ - A"S _ 0  0  0  -G"A"PS  (2.165)  n  For a system with a grounded base body n = 1, 7 = 7 I0  01  which can be evaluated from  equation (2.165). For a system with a floating base body n = 1, 7  ]0  can be developed from  equation (2.100), (2.103) and (2.113). Thus  0 0  GX-GX  (2.166)  0  where V^, and V \ are defined in equation (2.114) and (2.101). 1  h  2.7 Discussions The joint coordinate method developed above has been shown clearly to have some advantages. First, the final dynamic equations are expressed by independent generalized coordinates so that the computation is efficient compared with other methods expressed by redundant generalized coordinates. Second, the use of both absolute and relative coordinates makes the modelling of all kinds of force elements, control elements and joint constraints very  Chapter 2. New Formulations For DMBS  59  convenient. Moreover, the inertial matrices, the centrifugal and coriolis ( gyroscopic) force vectors, elastic force vectors and external generalized force vectors associated with absolute coordinates are formed in a systematic way independently for each body. Thus they are capable of being applied in a parallel computation implementation for large-scale system simulations. The velocity transformation matrix developed in section 2.4 can be constructed recursively from two different directions. This again makes the parallel computation implementation of velocity transformation feasible. The formulation method also allows us to extend the formulations easily into closed-loop configuration systems. Although there are some dependent generalized coordinates in the final dynamic equations for closed-loop systems, the number of dependent generalized coordinates is very small compared with the absolute coordinate method. And the joint coordinate method does not need to model joint constraints and solve differential-algebraic equations for tree topology MBS. The last but not the least is the formulations developed by joint coordinate method are implicit. This will ensure the numerical stability no matter what kinds of systems are dealt with, whether stiff or nonstiff systems. Deformable MBS do need this security to get accurate and efficient solutions. The velocity transformation matrix developed in section 2.4 has been decomposed into the terms that are associated with the properties of joint definition points and joints.  The  properties of joint definition points and joints which are reflected in matrices U,C,D,E or 0,C,D,E  are able to be evaluated independently for each joint and joint definition point. The  velocity transformation is then obtained by manipulating those matrices according to the body path matrix and joint path matrix. This not only makes it possible to implement parallel computation to a great extent but also reduces the possibility of redundant calculations. The velocity transformation matrix developed in section 2.6 also has the capability of being implemented in parallel computation because matrices  and J^ -\) n  n  can be calculated  independently for each body and joint. Although J _ and •/(„_])„ have simple forms and the n  x  Chapter 2. New Formulations For DMBS  60  topology description seems to be compact, the computation time for constructing velocity transformation matrices may not be less than that of the formulations in section 2.4 due to the numbers of factors for forming J  N  0  and its assembly procedure. An absolute comparison  between these two velocity transformation matrices is beyond the scope of this thesis. The main contribution in this chapter is the development of global velocity transformation which is described in section 2.4 and 2.6. The new topology definition or joint path matrix makes the derivation of section 2.4 feasible.  Chapter 3  Numerical Validation of the General Purpose Software 3.1 Introduction This chapter describes the capacities and the structure of the developed general purpose software. In order to verify the software, a two-link 3-D manipulator with one link rigid and the other link flexible is used as a test example. The tracked variation of the total energy of the system is compared to the strain energy to demonstrate that the strain energy is much larger than the numerical noise level of the total energy. The validation also includes using different numbers of elements for the flexible link. A three-link manipulator and an eight-link manipulator are tested to show that the software is able to deal with different configuration systems. Moreover, a modal representation of the flexible link and the effect of the small elastic deformations on the large rigid motion are investigated.  3.2 Software Implementation A general purpose program for time-domain simulation of deformable MBS based on the joint coordinate method developed in section 2.4 has been programmed into FORTRAN. The program deals with arbitrary configuration systems with 1-6 DOF ideal joints and with rigid or deformable bodies. Internal force elements such as springs and dampers, and external driving force elements such as PD controllers, concentrated forces and moments are also included in the software. The constraints for spherical joints, universal joints, cylindrical joints, revolute joints and translational joints are modelled to meet the requirement of dealing with closedloop systems. This code is not designed to solve constraint drifting problems which may occur in closed-loop systems although it does handle closed-loop systems. The input parameters include body and joint path matrices, the properties of each rigid or elastic body (although 61  Chapter 3. Numerical Validation of the General Purpose Software  ^  62  begin ^ \  input body number  A  7  join t definition poin t number  external&in ternal force number  call FLEX subroutine input flexible body properties  t Aall RIGID subroutine f  input mo dal coordinate number eigenvalue & eigenvector analysis] X  input rigid body propertied elastic coordinate reduction  /  call RELA subroutine  7  input join 1 properties  input body & joint path matrix join t definiion poin t position cut  joir^utmbgj^^g^^^^^  Stnput external & in ternal al^f X  force elements  z —4  initial condition  call DDASSL  call RES to form  •  dynamic equations  output  I  e n d  -)  to form Lagrange's equations  ^ for each body  t  V  I m velocity transformation trix&derivative matrix  final dynamic  integration solver  DAE  ^wmaimtraTnT equations  Figure 3.1 Flow chart of the developed software  equations  Chapter 3. Numerical Validation of the General Purpose Software  63  only beam elements are included currently), the positions of the joint definition points, the directions and number of DOFs of each joint, the positions and the types of internal and external force elements, initial conditions, etc.. The eigenvalue and eigenvector analysis of each elastic body is also incorporated. The eigenvalue analysis can be used for elastic coordinate reduction to increase computational efficiency. Users may choose how many modal coordinates they want according to the application. The numerical integration solver used is the code DDASSL which is the most widely-used solver in commercial MBS simulation software. DDASSL, which is based on backward difference formulas, is designed to be used for the solution of differential equations and differential-algebraic equations. A detailed introduction of DDASSL can be found in [58]. The total code includes more than 5000 lines excluding the DDASSL code. A flow chart of this program is shown in Figure 3.1.  3.3 Total Energy Validation for Two Different Examples This section presents a basic validation of the formulations and the code. The software and the dynamic formulations are general. The total system energy is tracked for a simple system. The system shown in Figure 3.2 is a two-link manipulator with one link rigid and the other one flexible. The kinematic parameters of the two links are listed in Table 3.1. The joint inertia has been included in the moment of inertia in y direction of link 1.  link 2  Figure 3.2 A two-link manipulator  Chapter 3. Numerical Validation of the General Purpose Software  64  , Link 1 (rigid) . steel, mass = 2.067 kg, length = 0.231m  aluminum, mass = 0.2007 kg, length = 1.02m  mass center: x = 0.0m, y = 0.0m, z = -0.1155mdensity p = 2770 kg In?, A = 7.104x 10~ m 5  c  c  c  moment of inertia: I^ = 0.04022 kgm  2  square tube 12.7mmxl2.7mmxl.6mm  I =0mS84kgm  2  modulus of elasticity E = 0.69 x 10 Pa 11  yy  l  =3.752 xXQi^kgm  2  u  modulus of rigidity G = 0.26 x 10 Pa 11  Table 3.1 The parameters of the two links  The joints are both revolute joints and are regarded as rigid. No control torques or friction is applied at either joint. At the tip of the second link, a concentrated mass m =05kg t  attached. The initial conditions are 0,(O) = L5rad and d (0) = -\.0rad 2  is  with no initial  deformation or velocity. The finite element method is used for modelling the flexible link as a Bernoulli-Euler beam with clamped-free boundary conditions. The elastic link is expressed by one beam element in this simulation. Thus the system has eight degrees of freedom in total, two rigid joint variables and six elastic variables. The six elastic variables are the translational and rotational deformations of the link tip. The system is a grounded chain system with three joint definition points. The body path matrix and joint path matrix are  n  o"  1 1 1  o"  x = "-1 0 -1 1 -1  The expected motion of the two joints would be harmonic if the second link were rigid and the motions were small. The use of an elastic link should not significantly affect the motion of two joints. The total energy of the system is set to be zero at the rest position.  Chapter 3. Numerical Validation of the General Purpose Software  65  Small structural damping is considered in the elastic link. The structural damping is usually modelled as Rayleigh or proportional damping to form a damping matrix [C] which is a linear combination of the stiffness and mass matrices, that is [59]  [C] = a[K]+B[M]  (3.1)  where a and /3 are, respectively, the stiffness and mass proportional damping constants. For structures that may have rigid-body motion, it is important that the mass-propoitional damping not be excessive [59]. Assume  is negligible, then  [C] = a[K]  (3.2)  The generalized structural damping force is  Q =-aKq d  (3.3)  f  The followings are the simulation results with structural damping a = 0.0003 and without structural damping. Their differences can not be distinguished in some figures due to the small value of the structural damping.  0.015  Figure 3.3 Elastic displacement along the axis direction of joint two (Z ) 2  Chapter 3. Numerical Validation of the General Purpose Software  Time (s)  Figure 3.4 Elastic displacement perpendicular to the axis direction of joint two (Y  z  Time (s)  Figure 3.5 Angular displacements of both joints (upper-joint 1; lower-joint 2)  Chapter 3. Numerical Validation of the General Purpose Software  1.0055  Time (s) Figure 3.6 Total energy  — with damping - - without clamping  0.03  0.025  — 0.02  <B 0.015  LU c CO  55 0.01 0.005  2  Time (s)  3  Figure 3.7 Strain energy  67  Chapter 3. Numerical Validation of the General Purpose Software  68  Time (s)  Figure 3.8 Kinetic and potential energy  The above simulation results show that the motions of the joints are as expected. The total energy in the conservative system (without damping) is a constant. The maximum absolute error of the total energy is about SxlO^Nm. The maximum strain energy is around 0.028 Nm, which is much larger than the maximum absolute error of the total energy. Comparing these two values lends credence to the accuracy of the simulation. The relative error of total energy should not be a standard to evaluate the accuracy of the simulation results although reference [25,28] both showed very small relative error. The total energy in the nonconservative system (with damping) shows that it indeed continues to decrease due to structural damping. The addition of structural damping makes the high frequency vibration die out (shown in Figure 3.4). Thus the numerical stability is improved. All the results show that the simulations are stable. Another check on the validity of the software is to replace the one beam element of the elastic link with 5 beam elements. The comparison of energies for the case of one element and  Chapter 3. Numerical Validation of the General Purpose Software  69  five elements are illustrated in Figure 3.9-3.11. The differences are hardly discernible in Figure 3.10 and 3.11. However, Figure 3.9 does show a small difference which indicates that the total energy for the one element case is larger than the total energy for the five elements case. This is consistent with the energy being a minimum for the exact solution since the fiveelement representation is more accurate than the one element representation.  —  E 1.0045^ 5^ o 1.004 LU 15  jS 1.0035 1.003 1.0025  0  2  3 Time (s)  Figure 3.9 Total energy comparison  one element ive elements  Chapter 3. Numerical Validation of the General Purpose Software  0.03 r  Y~ one element five elements  0.025 E 0.02 a» <j 0.015[ LU c CO ft 0.01 0.005 0  l  2  3 Time (s)  Figure 3.10 Strain energy comparison  Time (s)  Figure 3.11 Kinetic and potential energy comparison  70  Chapter 3. Numerical Validation of the General Purpose Software  71  Another example is a three-link manipulator with two revolute joints and one translational joint, two rigid links and one flexible link, shown in Figure 3.12. The properties of the links .are listed in Table 3.2.  Figure 3.12 A three-link manipulator  Link 1' (rigid) '  Link 2 (rigid)  steel, mass = 2.067 kg  Link 3 (flexible)  steel, mass = 2.067 kg  length = 0.231 m  aluminum, mass = 0.2007 kg  length = 0.462 m  mass center:  length = 1.02m  mass center:  density:  p = 2770 kg 1 m? A = 7.104x10~ m  x = 0.0m, y = 0.0/n, z = -0.1155m  x = 0.0m, y = 0.0m, z = -0.23 lm  area:  moment of inertia:  moment of inertia:  square tube:  c  c  c  1^ = 0.04022£gm  c  c  c  2  / „ = 0.04022 kgm  lyy =0.08884 kgm  1^ =0.08884&gm  2  1  4  zz  2  modulus of elasticity  I = 3.752 x 10"" kgm 4  a  2  12.1mm x 12.1mm x 1.6 ram  2  I = 3.752 xlO" kgm  5  2  E=0.69 x l o" Pa  modulus of rigidity G = 0.26 x 10 Pa  Table 3.2 The parameters of the three links  11  Chapter 3. Numerical Validation of the General Purpose Software  72  There are no control torques, friction orflexibilityat the joints. A mass m = 05kg is t  attached to the tip of the elastic link. The structural damping a = 0.0003 is included in the flexible link. The initial conditions are d =\5rad, x  Tj=0.0m, 6 =-h0rad 2  with all other  variables equal to zero. The simulation results are shown in Figure 3.13 ~ 3.17. The energy check and the responses show that the results are reasonable.  < -2.51 0  •  •  •  0.5  1  1.5  ..  1 2  Time (s)  Figure 3.13 Angular displacements of joint 1 and joint 3 (upper-joint 1; lower-joint 2)  0.14  Time (s)  Figure 3.14 Translation displacement of joint 2  Chapter 3. Numerical Validation of the General Purpose Software  Time (s)  Figure 3.15 Elastic displacement of link 3 (perpendicular to the axis direction of joint 2)  Time (s)  Figure 3.16 Elastic displacement of link 3 tip (parallel to the axis direction of joint 2)  Chapter 3. Numerical Validation of the General Purpose Software  14  ^1.0052 E >5 1.005  without damping with damping  \  CD  w 1.0048 CO  o l _  1.0046 0  0.5  *0.015  1 Time (s) i  1.5 — without damping - - • with damping  >  CO  c  ^ 0.005  'co 0.5  1 Time (s)  1.5  Figure 3.17 Total energy and strain energy  The software can handle tree topology systems as easily as chain topology systems,. The following example demonstrates the program applied to a system with PD control of a treeconfigured MBS. An eight-link manipulator, shown in Figure 3.18, is operated under a PD control law to reach a new configuration. The manipulator is composed of two flexible links and six rigid links. Six revolute joints, one translational joint and one spherical joint interconnect the links. The parameters of the links are given in Table 3.3. All the motors and PD control gains Eire the same. The motors have inertias of 2.0 kgm 1 rad each. The motors have mechanical and 2  electrical resistances of 0.5 kgm 1 rad I s. The gear ratios are all 1.0. The proportional and 2  derivative gains are 100 Nm/rad and 30 Nrn/(rad/s), respectively. The initial joint positions and final required joint positions are given in Table 3.4. The structural damping coefficient for the two flexible links is 0.00025. The simulation results are shown in Figure 3.19 ~ 3.24.  Chapter 3. Numerical Validation of the General Purpose Software  75  Figure 3.18 An eight-link manipulator  ' Link 1 & 2 (flexible)  •  aluminum, length = 1.0m  , Link 3-8 (rigid) steel, mass = 0.06125kg, length = 0.1m  density p = 2110 kg 1 rn, A = 7.104x10~ m  mass center: x = 0.05m, y = 0.0m, z = 0.0m  square tube 12.7 mm x 12.7mm x 1.6mm  moment of inertia: 1^ = 7.6563 X 10~ kgm  5  2  7  modulus of elasticity £ = 0.69xl0 Pa  I  modulus of rigidity  1^ = 2.0455 x 10" kgm  11  = 2.0455 x 10" kgm. 4  yy  4  G = 0.26xlO Pa H  Table 3.3 The parameters of the eight-link manipulator  joint 1  joint 2  joint 3  joint 4  (rad)  (rad)  (m)  (rad)  initial  1.3  0.4  0.0  final  0.4  0.6  -0.02  joint 5  ''joiht:6;- joint 7  joint 8;  (rad)  (rad)  (rad)  (rad)  0.0,0.0,0.0  0.0  1.0  0.0  1.0  0.5,0.5,0.5  0.3  0.5  0.7  0.4  Table 3.4 The initial and final positions of the joints  Chapter 3. Numerical Validation of the General Purpose Software  Figure 3.20 Displacements of joint 5 ~ 8 vs time (s)  76  Chapter 3. Numerical Validation of the General Purpose Software  77  0.15 E  Time (s)  Figure 3.22 Tip elastic displacement of link 1 in the z direction of the body coordinate system  Chapter 3. Numerical Validation of the General Purpose Software  78  0.06 Q.  * -0.04' 0  1  1  1  1  2  1  1  3  4  5  Time (s)  Figure 3.23 Tip elastic displacement of link 2 in the x direction of the body coordinate system  _ 0.04 r E. a. 0.03 CM  Time (s)  Figure 3.24 Tip elastic displacement of link 2 in the y direction of the body coordinate system  Chapter 3. Numerical Validation of the General Purpose Software  79  The results show that the elastic deformations of the two flexible links not only affect the tip positions of the manipulator but also affect the control of the joint positions. An effective control method needs to be provided to reach accurate positions.  3.4 Modal Representation The finite element method requires a large number of elastic coordinates. In many applications, the number of elastic coordinates is much larger than the number of joint degrees of freedom. In order to reduce the number of elastic degrees of freedom, the modal representation method can been employed. A few of the lowest frequency normal vibration modes and modal coordinates are used to represent the elastic deformations. This approach may lead to some errors due to the inability of normal vibration modes to account for local deformation effects induced by concentrated loads or constraint forces [38]. However, this method may have enough accuracy for some applications and significantiy decreases the number of degree of freedom required. In the method presented below, the elastic deformation of a flexible body i is represented by a linear combination of the component modes B multiplied by time-dependent generalized l  m  modal coordinates q'^, i.e.  (3.4)  q'f =  The modal matrix B' whose columns consist of linearly independent deformation modes is m  obtained by solving the eigenvalue problem of each deformable body only once. If the body i is assumed to vibrate freely about a reference configuration, then  m' q' K' q' =0 ff  f+  f  f  (3.5)  Chapter 3. Numerical Validation of the General Purpose Software  Define q'j = a'e ' ja  80  in order to solve the eigenvalue problem for equation (3.5).  [A-}-co /n^]a'=0  (3.6)  2  Solving the eigenvalue problem yields the eigenvalues (co^)  and eigenvectors a[ (  k -1,..., nj), where « / i s the number of DOFs of the elastic coordinates of body i. Let n be the lowest frequency modes, which are more likely to be excited, and discard the m  other higher modes, then  q' * [a, a • • • a„ ] q^ = B^q^ f  2  (3.7)  m  If n « n , the number of elastic coordinates are reduced significantly. m  f  The implementation of this coordinate reduction in the software requires little change in the formulations. The only change in the formulations is S of equation (2.14), i.e. lk  S  ik  = CSCBBB ik  ik  ik  ik  i  i  2  m  (3.8)  and the final values for elastic deformations are modal coordinates. For the previously-introduced case of the two-link manipulator with one link rigid and the other link flexible, the pre-processing of the software gives six eigenvalues and eigenvectors for the second link represented by one beam element. The frequencies and the modes respectively are shown in Table 3.5. The first two frequencies and modes in Table 3.3 indicate that two first-order bending modes in two perpendicular planes. Previous simulation results (section 3.2) show that only the first order mode is excited for the case of a two-link manipulator. Therefore, retaining only two modal coordinates should be sufficient to represent accurately the elastic deformation of the second link. The simulation results demonstrate that although the solutions are slightly  Chapter 3. Numerical Validation of the General Purpose Software  81  Freq. (Hz)  3.39  3.39  58.8  58.8  440.5  6604.4  Modes  0  0  0  0  1  0  -0.6653  0  0  0.0484  0  0  0  0.6653  0.0484  0  0  0  0  0  0  0  0  1  0  -1  1  0  0  0  -1  0  0  -1  0  0  Table 3.5 Eigenvalues and eigenvectors of a two-link manipulator  different numerically for the six elastic coordinate case compared with the two modal coordinate case, their joint motions exactly agree with each other, as shown in Figure 3.25. In other words, the modal representation method can reduce the number of elastic coordinates for more efficient computations while maintaining high accuracy in some applications.  c  < .2.51 0  1  1  .  .  2 3 Time (s)  1  1  4  5  Figure 3.25 Comparison of the six elastic coordinate case and the two modal coordinate case  Chapter 3. Numerical Validation of the General Purpose Software  82  3.5 Flexibility Coupling Effects Coupling between the large rigid body displacements and the small elastic deformations exists as the result of the finite rotations. This coupling is represented by non-linear terms in the mass matrix and the coriolis and centrifugal force components. The inertia tensor of the deformable body defined in the body coordinate system is not a constant matrix since it depends on the body deformation. The deformation may have a large effect on joint motions. This section gives a quantitative impression of the importance of this coupling effect through an example. The two-link manipulator discussed in section 3.3 is employed for this purpose. The response of the manipulator to an initial disturbance is compared between two cases; the one discussed previously in section 3.3 and a second one in which both links are treated as rigid. The comparison of simulation results for the two cases demonstrates the flexibility effect on the joint motions. The simulation results under the same initial condition as section 3.3 are shown in Figure 3.26 and 3.27.  T3 CO  1.6  — rjgid-flexible - - rigid-rigid  1  2  3  4  5  Time (s)  Figure 3.26 Comparison between rigid-rigid and rigid-flexible manipulator (joint 1)  Chapter 3. Numerical Validation of the General Purpose Software  83  Time (s)  Figure 3.27 Comparison between rigid-rigid and rigid-flexible manipulator (joint 2)  The deformation effect on the joint motions is most clearly revealed in Figure 3.26. The magnitude difference indicates that a part of the system total energy has been used for overcoming elastic deformation. The magnitude of the difference depends on the magnitude of deformation of the flexible link. The small difference in Figure 3.27 shows thai the deformation in the perpendicular direction of joint 2 is very small due to the rigid motion in this direction. Moreover, Figure 3.26 and 3.27 show that the frequency is shifted. The small difference in frequency is caused by the effect of small elastic deformation on the inertia matrix, gyroscopic and centrifugal forces. The comparison results show that the elastic deformations have a large effect on the joint motions and ignoring elastic deformations results in modelling errors.  3.6 Summary Numerous examples of flexible MBS simulations have been presented, based on the formulation developed in Chapter Two. Two flexible manipulators were used for the  Chapter 3. Numerical Validation of the General Purpose Software  84  validation of the newly developed formulations and the program. The absolute errors in total energy and the strain energy were compared to ensure that the small elastic deformations are not buried within the numerical noise. The results for the two manipulators show that the total energies are conserved within the simulations and the strain energies are several orders larger than the absolute errors of the total energies. These demonstrations support the accuracy of the formulations and the code. One and five element FEA discretization of a flexible link within a MBS have been compared. The example of the eight-link manipulator under PD control has demonstrated that the program can deal with tree-configured systems. Modal reduction was discussed and comparisons presented. The results show that modal reduction is able to reduce the number of elastic coordinates without compromising accuracy. The effect of elastic deformations on the joint motions is investigated numerically. The simulation results show that the magnitudes and the frequencies of joint motions are affected by the elastic deformations.  Chapter 4  Experimental Validation 4.1 Introduction The validity of any numerical modelling is greatly strengthened when compared and confirmed against experimental data. Several experiments in 2-D DMBS have been conducted [62,73]. Unfortunately, no 3-D experiments to verify the simulation modelling of deformable MBS have been found in the published literature. This chapter covers the description of an experimental apparatus, a series of experiments performed with that apparatus and comparisons with the numerical simulations performed with the numerical modelling previously described.  4.2 Physical Description The experimental apparatus is shown in Figure 4.1. It corresponds to one of the numerical models described in the previous chapter. The test rig consists of two links, joined by a revolute joint and attached to the ground by another revolute joint. The revolute joint which attaches the first link to the ground is implemented through a pair of precision ball bearings mounted one above the other. Thefirstlink, made of a solid steel rod, is constrained to rotate in a horizontal plane by thefirstjoint. The second joint is also implemented through another pair of precision ball bearings aligned with the axis of the first link. The second link is constrained to rotate in a vertical plane, perpendicular to the axis of thefirstlink. The second link is a hollow, square, aluminum tube of sufficiently small bending stiffness to permit significant flexible movement. A mass (0.604kg) is attached to the tip of the flexible link. The objective of the design was to produce enough strain at the jointed end of the flexible link to be measurable. The resolution of strain gauges is on the order of one micro m/m. The noise 85  Chapter 4. Experimental Validation  86  level of a good circuit is around ten micros. Therefore the strain signal needs to be at least one hundred micros. Simulation calculations were performed using the above-described software to check the strain level. The link parameters are shown in Table 3.1.  Figure 4.1 The test rig  4.3 Instrumentation Two high precision potentiometers (PRECISION MIL STYLE RV4 with 7/8" long shaft standard bushing linear taper and 5 KQ., manufactured by Electrosonic ) were mounted at the joints to measure the angular motions of the rigid link relative to ground and of the flexible link relative to the rigid link. The four precision strain gauge compacts (CEA-13-062UT-350, manufactured by Micro-Measurements), each of which consists of two mutually perpendicular strain gauges, are mounted on four sides of the flexible link near the second joint. The four strain gauge compacts are wired to form two full bridges (8-Channel Strain Gauge Signal Conditioning SC-2043-SG, manufactured by National Instruments) to measure the strain  Chapter 4. Experimental Validation  87  signals and hence the beam curvature in two different directions. The wiring of one direction is shown in Figure 4.2. The SC-2043-SG is a signal conditioning board with amplifier gain of 10 that interfaces a National Instruments Data Acquisition Board (DAQ) directly to the strain gauge transducers. The potentiometers share the DC power with the strain gauges. Together, the signals from the potentiometers and strain gauges are connected to the DAQ card (LabPC-1200, manufactured by National Instruments), which in turn is read by the computer (486DX33) through the code programed by my supervisor, Dr. Dunwoody. The sampling rate is 200 Hz. The whole system is shown in Figure 4.3.  Figure 4.2 Strain gauge wiring  Computer  'train Gauges 'otentiometer  Figure 4.3 Test devices  Chapter 4. Experimental Validation  88  4.4 Calibration Experiments Before any experiments could proceed, the potentiometers and strain gauges were calibrated. Both potentiometers use the same calibration method. Take potentiometer 1 as an example. A line was drawn on the ground and the rigid link respectively to form two datum lines. When the two datum lines meet, the relative angle was noted as zero degree and the voltage was recorded. Different lines which represent - 6 0 ° , - 3 0 ° , 30° and 60° angles relative to the datum line of the ground also were drawn on the ground. The voltages were recorded when the datum line of the rigid link meet those angular lines. The following are the results of the calibrations. The stars represent experiment data points and the straight line is a least-squares linear interpolation.  Angle (deg)  Figure 4.4 The calibration of potentiometer 1  The calibration equation for potentiometer 1 is d = (1.27- Vj)* 100.0 deg . Where V is the x  recorded voltage in volts.  l  Chapter 4. Experimental Validation  89  Angle (cleg)  Figure 4.5 The calibration of potentiometer 2  The calibration equation for potentiometer 2 is 6 = (V -1.26) *115.0 deg . Where V is the 2  2  2  recorded voltage in volts. The strain gauges were calibrated using the following method. The flexible link was placed horizontally and clamped at the end on which the strain gauges were mounted. Different forces were applied at the other end. The magnitudes of the strains caused by the forces were calculated and the voltages were recorded.  Chapter 4. Experimental Validation  90  The calibration equation for the strain gauges is st= v/36.7- 0.00007 m/m. where v is the recorded voltage in volts. In addition, the structural damping of the flexible link was evaluated experimentally. One end of the flexible link was clamped and an accelerometer was attached onto the tip of the other end of the link. The tip was given an initial displacement 0.02m, then released. The method is demonstrated in Figure 4.7 and the test result is shown in Figure 4.8.  Figure 4.7 Damping measurement method  Time (s)  Figure 4.8 Tip acceleration of the experiment  Chapter 4. Experimental Validation  91  As demonstrated in equation 3.2, the structural damping matrix is equal to the damping coefficient times the stiffness matrix. Therefore, the governing equation for the vibration of the flexible link with one end fixed is:  [M]{y}+a[AT]{y}+ [K]{Y}= {0}  (4.1)  The flexible link was discretized into three beam elements. The flexible link was given an initial displacement of 0.02m and released. The resulting tip acceleration is shown in Figure 4.8. The same conditions were simulated using different values of a until the simulated response matched the measured response. The damping coefficient a was found to be 0.00025.  Time (s)  Figure 4.9 Tip acceleration of the simulation  4.5 Comparison Between Experiment and Simulation Results This section covers the experiments performed with the test rig under different initial conditions, the modelling of the test rig and the comparison between the experimental and simulated results.  Chapter 4. Experimental Validation  92  4.5.1 The Modelling of the Test Rig The test rig is similar to a manipulator without control motors at the joints. For many manipulators in the physical world, jointflexibilityis significant. In addition to the torsional flexibility of the gears, joint flexibility is caused by effects such as shaft windup and bearing deformation [60]. In this test rig, bearing deformation is the main cause of joint flexibility. Torsional springs are used to represent the joint flexibilities of the two joints. Moreover, damping forces exist at the two joints in the real test rig. They can be better modelled as viscous dampers. The relationship between joint definition points, torsional spring and damper is illustrated in Figure 4.10. Link 2_ Link!  Jo int Definition Po int s Figure 4.10 loint damping and flexibility  Therefore, the test rig can be modelled as two bodies with one rigid and the other one flexible, connected by two revolute joints with two rotational springs and dampers. The tip mass can be considered as a point mass attached at the tip of theflexiblelink. The test rig can be modelled as:  B MBq + C q + K,q = B (Q + Q + Q +Q T  T  d  where C and K D  T  V  C  S  0  MBq)  (4.2)  are the viscous damping and torsional stiffness matrices. Most of their  elements are zeros except the elements that represent  the damping and stiffness  corresponding to the joints. The viscous damping and torsional stiffness were determined by trial and error. Values were chosen to ensure that the trajectories of the modelled joints  Chapter 4. Experimental Validation  93  agreed with the experimental results under one test condition (0, = 0 ° and 6 = 3 0 ° ) as 2  closely as possible. For this test rig, the damping and stiffness were found to be: C\ = C] = 0.098 Nm I {rad I s), K] = K = 0.1 Nm I rad . 2  t  The modelling parameters for the links of the test rig is listed in Table 3.1. The equivalent principal moments of inertia for link one not only include the inertia of the link itself but also the inertia of the two joints. The tip mass is 0.604 kg.  4.5.2 Experiments Performed and Comparison with Simulations The experiments were performed starting with two different initial conditions. No forces were applied after release in either case. In the first experiment, the joint angles were given initial values of 0, = 0° and 6 = 30°. The joint displacements and flexible link strain of the 2  simulation and the experiment are shown in Figure 4.11 ~ 4.13. They have demonstrated a very good agreement between simulation and experiment results. There appears to be a small zero-sliift problem with the experimental results of the joints, shown in Figure 4.12. The friction of the potentiometers or the joints is probably the main cause of this problem. 40  Time (s)  Figure 4.11 Response of joint 1 at initial condition 0, = 0° and 6 = 30° 2  Chapter 4. Experimental Validation  -40 -2  94  L  2  4 Time (s)  6  10  Figure 4.12 Response of joint 2 at initial condition 6 = 0°and 6 = 30° l  2  x10 • experimenfl simulation  3 2  h c  '5 0 CO -1  -3  0  2  4 Time (s)  6  8  10  Figure 4.13 Jointed-end strain of the second link at initial condition t9, = 0° and 0 = 30 ( Z ) H  2  2  Chapter 4. Experimental Validation  95  A frequency domain comparison of the strain signal between the simulation and the experiment is illustrated in Figure 4.14 and Figure 4.15. From Figure 4.14 and Figure 4.15, we can clearly see that there are several principal frequencies in the bending response of the elastic link. They are scattered around 1.1 Hz, 2.2 Hz, 4.0 Hz shifting to 3.4 Hz, and 5.5 Hz. The essential reason for this behavior is the resultant of the nonlinear coupling between the large joint motions and the small elastic deformations. A detailed explanation can be given with the help of the study of a nonlinear vibration system, a single degree of freedom pendulum with a large angular displacement. Hagedorn [77] showed that the frequency of free oscillation of this pendulum is not a constant as with linear vibration but is a function of the amplitude. Also, in the forced oscillation case, an essential difference between linear and nonlinear systems is the fact that the latter present not only periodic oscillation with the excitation frequency but also periodic oscillations with other frequencies [77]. Superharmonic oscillations, oscillations with the frequencies being multiples of the excitation frequency, or subharmonic oscillations, oscillations with the frequencies being fractions of the excitation frequency, often can be observed. The test rig is much more complicated than that pendulum. However, the mechanisms are similar. The motion of the flexible link of the test rig can be deemed as a nonlinear pendulum. The action of the first link on the flexible link can be thought as a forced vibration, background oscillation. The excitation frequency which can be found in Figure 4.11 or 4.12 is about 0.55Hz. Therefore, according to the forced vibration result of Hagedorn [77], the superharmonic vibrations could be found at the frequencies 1.1 Hz, 2.2 Hz and 5.5 Hz. Moreover, the free vibration frequency of the flexible link depends on the angular amplitude of the flexible link. The angular amplitude of the flexible link decreases with time due to the joint friction. Thus, the free vibration frequency shifts from 4.0 Hz to 3.4 Hz in which 3.4 Hz is the free vibration frequency of the flexible link under no large angular displacement. This also indicates that the assumed vibration mode modelling method has some approximation and may cause errors in thefinalresponses. The conclusion can be drawn that  Chapter 4. Experimental Validation  96  the responses agree very well also in frequency domain after comparing Figure 4.14 with Figure 4.15.  Chapter 4. Experimental Validation  97  The second experiment used initial conditions of 0, = 0° , 6 = 45°. The time and frequency 2  domain comparisons between experiment and simulation are shown in Figure 4.16-4.20.  Chapter 4. Experimental Validation  98  Chapter 4. Experimental Validation  99  From Figure 4.18, digitization noise was present in the strain signal of the second experiment result but not of the first experiment. This was because the gain of DAQ was chosen differently. From Figure 4.16 ~ 4.20, we can see that the experiment and simulation results do not agree with each other very well, especially in Figure 4.19 and 4.20. The angular displacement error of the second joint which may be caused by the calibration experiment or not enough initial angle probably is the reason. Figure 4.20 ~ 4.24 give the comparisons of the experiment results under the initial conditions of 0, = 0° and 0 = 45° and the simulation 2  results under the initial conditions of 0, = 0° and 0 = 40°. They show good agreement. The 2  regularity in the frequency domain still can be explained as in thefirstcase (initial conditions 0, = 0 ° and 0 =30°). 2  Chapter 4. Experimental Validation  Chapter 4. Experimental Validation  -40' -2  >——  0  101  1  2  1  4  i  i  6  8  10  Time (s) Figure 4.23 Response of joint 2 -4  x  _6i  -2  10  1  0  1  2  1  4  i  i  I  6  8  10  Time (s)  Figure 4.24 Joint-end strain of the flexible link ( Z ) 2  Chapter 4. Experimental Validation  102  4.6 Comparison with the Simulation Results of ADAMS The test rig was also modelled by using the commercial software program ADAMS (version 8.2). The rigid link, the two revolute joints and the two rotational springs and dampers were easily built up by creating bodies, constraints and force elements. Flexible bodies can be modelled by incorporating FEA or by importing some matrices of the modal representation. For simple shape flexible bodies, ADAMS designs force elements, such as BEAM and FIELD elements, to represent elastic deformations and structural damping. These force elements are massless. Therefore, the inertias of the flexible bodies are considered as the rigid inertias. This means that ADAMS ignores the elastic coupling effect in inertia matrices and deems the inertias of the flexible bodies to be constants. The flexible link was modelled by a BEAM element. The tip mass was modelled as a point mass connected to the flexible link by a spherical joint. The joint simulation results of ADAMS under the same initial conditions were also compared with the experiment results, shown in Figure 4.25 ~ Figure 4.28.  Chapter 4. Experimental Validation  103  Chapter 4. Experimental Validation  104  Compare Figure 4.11-4.12 and Figure 4.16-4.17 with Figure 4.25-4.28, we can see that the joint displacements of the developed program and ADAMS are different. Thefirstpeak of the angular displacement of joint 1 obtained by ADAMS is almost 10 degree larger than that of the experiment and that of the developed program. The modelling difference between ADAMS and the developed program is the cause of this. ADAMS uses a lumped mass method and ignores the inertia coupling between rigid and elastic motion. The elastic deformations obtained by ADAMS can not be compared with the experimental results directly because strain signals were measured in the experiments but elastic displacements were obtained from ADAMS. An indirect method is to use the elastic translational and rotational deformations obtained from ADAMS to get strain signals. Figure 4.29 and Figure 4.30 show the comparison results under the two initial conditions.  Chapter 4. Experimental Validation  105  -4  x 10  Time (s) Figure 4.29 Strain signal comparison at initial condition 0, = 0°and 0 = 30° 2  x 10 1.5  I  --| I  -2  1  1  1  1  0  2  1  1  r  i  i  i  1  4 Time (s)  6  8  10  Figure 4.30 Strain signal comparison at initial condition 0, = 0°and 0 = 45° 2  Chapter 4. Experimental Validation  106  The above results show that the strain achieved by ADAMS is larger than that by the experiment. Also, Figure 4.31 ~ 4.34 show the frequency domain results of ADAMS and the experiment under the two initial conditions. There are also differences in the frequency components and magnitudes due to the modelling method and experiment errors.  Frequency (Hz)  Figure 4.31 Spectrum of the strain signal by ADAMS at 0, = 0°and 0 = 30° 2  0  2  4 6 Frequency (Hz)  8  10  Figure 4.32 Spectrum of the strain signal by experiment at 0, = 0°and 0 = 30° 2  Chapter 4. Experimental Validation  107  4.7 Summary The experiment and simulation results have been compared under two test conditions for a test rig. The angular displacements of the two joints and strain signal of the flexible link were  Chapter 4. Experimental Validation  108  measured and simulated. The good agreements in both time and frequency domain demonstrate that the formulations are accurate. The frequency domain responses also illustrate that the elastic vibration frequencies of the elastic link change with the amplitude of the motion of the flexible link. The nonlinear coupling between the joint motions and the elastic motions is displayed distinctly in the frequency domain. The angular displacements and the strain signal of the test rig as simulated by the commercial program ADAMS are also compared with the experiment results.  Chapter 5  Simulation Comparison of Different Formulation Methods  5.1 Introduction This chapter presents the simulation comparison between different formulation methods. An Order N formulation is developed and implemented in software according to the method Keat [20] proposed. The accuracy and efficiency of the Order N method are compared with that of the joint coordinate method developed in Chapter 2. Moreover, the simulation results of the joint coordinate method and the absolute coordinate method (ADAMS) are also compared. The chaotic behavior of nonlinear systems is illustrated.  5.2 Recursive or Order N Formulations In recent years, most of the dynamic equations for deformable MBS have been derived by recursive or Order N method due to the its low operation count per time step [2022,24,26,27]. Although the Order N method is regarded as the most efficient method when comparing the number of operation per time step, the Order N method is essentially an explicit method which needs very small time steps to ensure the stability for flexible MBS. Therefore, it may not be faster than other methods which are implicit when the comparison is a fixed period of time to simulate. Moreover, the numerical truncation error and round-off error from inverting the inertia matrix of each body will be enlarged and propagated to the next body. The numerical errors may become significant when large systems are involved or poorly conditioned inertia matrices exist. Thefirsttask in this chapter is to demonstrate a detailed formulation of the recursive or Order N method using finite elements. In order to compare both formulations in the same conditions, the developed velocity transformation matrix and its  109  Chapter 5. Simulation Comparison of Different Formulation Methods  110  derivative matrix are used to derive a new recursive or Order N formulation as was done by Keat [20]. The main idea of the Order N method is to evaluate the relative and elastic variables explicitly by inverting the inertial matrix for each body sequentially [22]. In section 2.4 of Chapter 2, it has been shown that the velocity transformation matrix B and its derivative matrix B are always lower triangular matrices which are composed of many block matrices H  IS  and H no matter what kind of topology the system has. That is TJ  fi = [//y]  and  = [#,;,•  Then from equation (2.87)  where  (5.2)  T = Q + Q +Q V  c  (5.1)  + Q  s  (5.3)  0  T  Multiplying equation (5.2) by B systems yields  and combining equation (2.130) and (2.143) for tree  1 ^^ [ Y = Bq + Bq Nb  _  .  .  Nb  = 1  I H MY T  ]  n  H  or  ]  (5.4)  „  .  Hj r  J  t  H  (5.5)  Y=iH q +iH q j  j  ij  ij  J=\  7=1  In order to explicitly solve for the relative and elastic variables, start with the final body, i = N ,i.e., b  ijT NbNb  H  t, NbyNb _ TJT Nb ~ NbNb r  M  Y  M  l  (5.6)  Chapter 5. Simulation Comparison of Different Formulation Methods  Nb  111  Nb  = 1 H qJ + I H qJ  Y  Nb  Nbj  (5.7)  NbJ  7=1  7=1  Equation (5.7) may be substituted into equation (5.6) to solve for q , Nb  ..Nb  :  q  ,NbU ,  (H NbNb M" H  Nb-\  T -r^Nb „T .Nb NbNb^ NbNb '  I H q  J  t  H  NbNb  H  M  NbJ  \)=\  Therefore, q  is a function of q  Nb  Nb  reveals that Y  Y  Nb  —M  ,q  Nb  1  m  (5.8)  7=1  ,---,q . Substitution of equation (5.8) into (5.7) X  is also a function of q ^, q ~ , • • •, q .  N b  1  1  Nb  + I H q*  Nb  (H  T  " NbNb\  n  I ~ H (H  l  v  l  n  NbNb)  \  2  r  1  T  Nbm  2  1  j.  Nb  NbNb  H M  Nb  NbNh  U  T  1  M H )  T  NhNb  M U Nh  NbNb  Nb  Nb-l  Nb  7=1  7=1  Nb  (5.9)  NbNb  The same procedure may be repeated in turn for every other body, ending with the root body. Finally, the general equations can be obtained as i-l  ,  Nb  Miq^F^lT^+lP^q ••k , v n -k k=l k=l Nb  where  (5.10)  Nb  M; = ^H^M'H,  - £ ( / , ) (M*)' 7,  (5.11)  7=i+l  7='  Nb  Nb  (5.12) k=j  k=7+1  Nb  ^ = I ^ 7=<  r  7 - / :  (5.13)  Chapter 5. Simulation Comparison of Different Formulation Methods  112  Nb  (5.14) j=M  Nb  (  T =  (5.15)  Y,HlM>H +t jk  ik  Nb  (5.16)  Nb  P =~ ik  J,HlM'H  jk  +  P i k  (5.17)  Nb  (5.18) To solve for the joint accelerations, equation (5.11) ~ (5.18) must be solved in order from i = N to / = 1. Then equation (5.10) is computed sequentially from i = l to i = N . The b  b  procedure is analogous to inverting a matrix by LU decomposition, but without pivotting. From equation (5.11) ~ (5.18), we can see that many inertia matrices (M*) need to be inverted and the inverse matrices are multiplied by other matrices to form other inertia matricesfromone body to the next. This is a common characteristic of Order N or recursive methods. As stated previously, the accuracy of the solution may suffer as a result of the repeated inversions.  5.3 Simulation Comparison The following comparison is to demonstrate problems in explicit formulation methods when numerical simulation proceeds.  Chapter 5. Simulation Comparison of Different Formulation Methods  113  5.3.1 Numerical Validation of Derived Order N Formulations The Order N formulation was subjected to similar validation procedures as the previous implicit formulation. The matrices H ,H ,M ,T i  ji  ji  i  are the same matrices as in the implicit  formulation. If the system has only one body, the Order N and implicit formulations are same. Therefore, the simulation results should be same if the same implicit numerical integration software (DDASSL) is used. A rigid link with tip mass attached was tested. The two formulations produced identical results. Next, a two-rigid-link mechanism with two revolute joints similar to the example described in Chapter Three was tested. The properties of the mechanism are listed in Table 5.1. The initial conditions were 0, = 1.5rad,0 = -1.0rad. The 2  tip mass on the second link was 0.5kg. The results of the Order N formulation compared with the implicit formulation are shown below. The implicit numerical integration solver DDASSL was used in both implicit and Order N formulation simulations.  iM&"Link  f (rigidX  '' •''''pS^f;|S ; ;  Link 2 (tigid)  mass =1.76 kg  mass =0.49 kg  length = 0.2 m  length = 0.2 m  mass center:  mass center: x = 0.1m  x = 0.1m  y = 0.0m  y = 0.0m  z = 0.0m  z = 0.0m  c  c  c  c  c  c  moment of inertial:  moment of inertial:  I =3.92xl0~  4  xx  kgm  2  I „ =2.45x10' kgm 5  Iyy = 0.0236kgm  1^ =0.0065 kgm  l  l  2  2  =0.0236 kgm  2  zz  2  = 0.0065kgm  2  zz  Table 5.1 The parameters of a two-rigid-link mechanism  Chapter 5. Simulation Comparison of Different Formulation Methods  2  -2.5 0 1  ' 1  1  2  Time (s)  • 3  ' 4  1 5  Figure 5.1 Angular displacements of the two joints (upper-joint 1; lower-joint 2)  —  0.855 h  implicij  - - exphcil  E 0.85 [ >. » CD C LU  « 0.8451  .o 0.84 k / 2  Time (s)  3  Figure 5.2 Total energy  Chapter 5. Simulation Comparison of Different Formulation Methods  115  The main reason for the small differences between the results of the Order N and implicit formulations is the error propagation from matrix inversion, which can be indirectly proved using two different inversion methods, as demonstrated in the next section. In this case, the differences are very small compared with the kinetic and potential energies. These tests indicate that the coding of the Order N formulation does not contain errors. The formulation of equation (5.10) ~ (5.18) do not distinguish between rigid and deformable systems. Therefore, no separate tests are required using a deformable system.  5.3.2 Numerical Comparison The comparison between explicit and implicit formulations proceeds under the same conditions (properties and initial conditions) as in Chapter 3 with the structural damping of the flexible link being a = 0.0003. The simulations were conducted using three different combinations. One was the implicit formulation developed in Chapter 2 combined with the implicit numerical integration solver DDASSL. These results are abbreviated as Im-Im; The  Chapter 5. Simulation Comparison of Different Formulation Methods  116  second one was the Order N formulation derived in this chapter combined with the implicit numerical integration solver DDASSL. These results are noted as Ex-Im. Thefinalone was the Order N formulation combined with the explicit numerical integration solver DIVPRK. These results are labeled as Ex-Ex. DIVPRK is a Runge-Kutta-Vernerfifth-order& sixthorder explicit integrator. The step size was 1 x 10" second in this case. DDASSL is a variable 6  -step-size solver based on the backward difference formulas. In this case, the step size changed from 1 x 10" to 0.01 second. The three combinations were simulated at the same 6  accuracy requirements, which were an absolute error limit atol =10" and a relative error limit 7  rtol = 10~ . The results are shown below. 4  Time (s) Figure 5.4 Angular displacements of the two joints (upper-joint 1; lower-joint 2)  Chapter 5. Simulation Comparison of Different Formulation Methods  0.03 — —  Im-lrr Ex-lrr Ex-E>  w -0.01  -0.02  0  1 Time (s)  Figure 5.5 Elastic displacements of second link tip along the axis direction of joint 2  117  Chapter 5. Simulation Comparison of Different Formulation Methods  118  Im-lrr Ex-lrr Ex-E>  2  3 Time (s)  Figure 5.7 Total energy  0.12  r  — Im-lrr — Ex-lrr — - Ex-E>  0.1 | 0.08 \ >> g 0.06 [  LU  I 0.04 \ CO  0.02  ,l/A. 2  3 Time (s)  Figure 5.8 Strain energy  Chapter 5. Simulation Comparison of Different Formulation Methods  •1.5  1  \  1  1—  cn 1  o0-5  CO g  *  0  -j \J 2  — — >v  119  Im-lrr Ex-lrr E x  E >  3 Time (s)  2  3 Time (s)  Figure 5.9 Kinetic and potential energies  The above results (Figure 5.4 ~ 5.6) show large differences not only in elastic displacements but also in rigid-body motions. The total energy and energy component check give us a clear picture which one is the most accurate. The tested system is nearly a conservative system since the energy absorbed by structural damping is very small. The total energy shown in Figure 5.7 for the Ex-Im and Ex-Ex cases do not show conservation of total energy. The maximum deviations in total energy for both Ex-Im and Ex-Ex results are of the same order of magnitude as the maximum strain energy shown in Figure 5.8. Therefore, the elastic displacements of both Ex-Im and Ex-Ex formulations are suspect because they are within the noise levels of these formulations. Moreover, the computational times for these three combinations were diverse. The Im-Im simulation took less than ten minutes on a PC to simulate the dynamic responses for thefirst5 seconds compared with more than two hours for the Ex-Im simulation and nearly thirty hours for the Ex-Ex simulation. The Order N formulation did save calculate time per time step. The Ex-Im simulation needed many iterations per step and the Ex-Ex simulation needed a very small time step (about 1 x 10~*  Chapter 5. Simulation Comparison of Different Formulation Methods  120  second in this case). The Im-Im simulation calculated less than 4000 points to reach 5 seconds but the Ex-Im simulation calculated more than 16000 points and the Ex-Ex simulation calculated 5 millions points. From the above simulation results using the Ex-Im combination, we can see that the Order N formulation has convergence problems as shown by the many iterations required. Inaccuracies in the matrix inversions could be the cause. The coupling between the large joint motions and very small elastic motions may result in the matrices M* being ill conditioned. The inversion of matrices M* can cause roundoff errors which are enlarged by J and propagated to M*. The main danger of ill-conditioning is not that the jt  solution may fail, but that it may succeed yet produce a solution whose errors are serious but not large enough to make it obvious that something is wrong [59]. The Order N formulation method is very sensitive to errors in matrix inversion. Different inversion methods or different software can produce different truncation errors and roundoff errors. These errors are enlarged step by step so that the dynamic responses can become significantly different. The following example demonstrates the large differences which can be caused by different inversion software. The physical modal is still a two-link manipulator whose parameters are described in Chapter 3. The initial conditions are the same as the above example. The Order N formulation method with the implicit numerical integration solver DDASSL is investigated using two different inversion methods, DLINRG and MTINV. Both methods use double precision. DLINRG computes an LU factorization of the coefficient matrix and the inverse matrices of L and U respectively. Then the final inverse matrix can be represented by the multiplication of the inverse matrices of U and L. MTINV employs GAUSS elimination to compute the inverse matrix. All the simulation conditions are the same except for the two different inversion methods. The simulation results are shown in Figure 5.10-5.12. The large differences between these two cases indicate that the Order N formulation method is sensitive to the method of inverting matrices. In other words, small numerical errors of inverse matrices induces large differences in the solutions. The same test was done with the implicit  Chapter 5. Simulation Comparison of Different Formulation Methods  121  formulation using DLINRG and MTINV. The results show that there is no difference in the figures of the two cases.  1.06  0.98' 0  1  1  •  1  2 3 Time (s)  • 4  1 5  Figure 5.10 Total energy  Time (s)  Figure 5.11 Angular displacements of the two joints (upper-joint 1; lower-joint 2)  Chapter 5. Simulation Comparison of Different Formulation Methods  122  0.025 0.02 .§. 0.015 CO  •£  CD o TO  Jj" 0.005 a o  I  0  LU  -0.005 -0.01 0  1  2  3  4  5  Time (s)  Figure 5.12 Elastic displacement of second link tip along the axis direction of joint 2  5.3.3 Comparison with ADAMS The simulation results of the test rig obtained by ADAMS have been compared with the experiments and the results of the developed program in Chapter Four. However, the differences in dynamic responses was not clearly illustrated due to the flexibility and damping of the joints of the test rig. This section covers a comparison of the simulation of the test rig without joint flexibility or damping achieved by the developed program and ADAMS. The purpose is to analysis the difference. The simulations of the above model were under the initial condidtionfl, = 0° and 8  2  -- 30°.  The results are shown in Figure 5.13-5.16. The legend "simulation" on the figures indicates the results obtained by the developed program. The results show that there are large amplitude differences not only between elastic displacements but also between joint angular displacements. However, the frequencies are the same. The reason can be understood by comparing the modelling and formulation methods.  Chapter 5. Simulation Comparison of Different Formulation Methods  123  ADAMS uses the absolute coordinate formulation method with generalized coordinate partitioning to eliminate the redundant variables and produce pure differential equations containing only independent generalized coordinates. The method was described by Shabana [3]. The algebraic equations of all joints in the system are constructed and equations and variables are partitioned into two parts, one of which has the number of the degree of freedom of the system. The partitioning method calculates the determinant of the coefficient matrix to find a set of independent coordinates. The efficiency is a problem when large systems involved. Moreover, the implicit fomulation becomes an explicit formulation after the dependent coordinates are represented by independent coordinates since the terms involved dependent coordinates are explicitly expressed by independent coordinates. However, the formulation is different from the Order N formulation because the inertia matrix for the entire system is inverted once rather than for one body at a time as in the Order N method. Therefore, the numerical inverse errors do not propagate from one body to the next. But the numerical instability can be seen in Figure 5.15 which shows the elastic displacement amplitude changes a little. ADAMS models the flexible link through a force element, wliich is a massless beam element in this case. The inertia of the flexible link is considered to be the inertia of a rigid beam, ignoring the elastic inertia and rigid-elastic coupling inertia. The flexible link is considered to have a lumped inertia. The elastic and damping forces represented by the beam force element are applied directly to the lumped inertia. This is the main factor that causes the difference between the displacement magnitudes.  Chapter 5. Simulation Comparison of Different Formulation Methods  Time (s)  Figure 5.13 Angular displacement of joint 1  Time (s)  Figure 5.14 Angular displacement of joint 2  124  Chapter 5. Simulation Comparison of Different Formulation Methods  0.02  2  3 Time (s)  Figure 5.15 Tip elastic displacement along the axis direction of joint 2  x10'  2  3 Time (s)  Figure 5.16 Tip elastic displacement perpendicular to the axis direction of Joint 2  Chapter 5. Simulation Comparison of Different Formulation Methods  126  5.4 The Chaotic Behavior in Simulation Results due to the Formulation Method The simulation of the Order N formulation method may become chaotic due to the propagated and enlarged numerical error. The unique character of chaotic dynamics may be seen most clearly by imagining the system to be started twice, butfromslightly different initial conditions. If one chaotic mechanical system is started at initial conditions x and x + e respectively, where e is a very small quantity, their dynamic responses will diverge from each other very quickly. This phenomenon, which occurs only when the governing equations are nonlinear, is known as sensitivity to initial conditions. The simulations have been rerun with slightly different initial conditions, 0, = 1.5 rad and 6 = -1.02 rad , to test the implicit and 2  Order N formulations. The Im-Im case shows that the differences in responses between the two simulations can hardly be seen in the graphs. The response differences for the Ex-Im case are much larger as shown in Figure 5.17 ~ 5.20. The irregular divergence in elastic displacement (Figure 5.18) and total energy (Figure 5.19) show that chaotic behavior does exist in the Ex-Im case. Usually, if governing equations which have at least three independent variables are nonlinear, the possibility of chaotic behavior exists. The condition for chaos to occur depends on the parameters of the system. For the system with the parameters we investigated above, the dynamic responses are not chaotic since the implicit formulation simulation has not shown any irregular or unpredictable motion. Chaotic behavior occurred with the explicit formulation method, demonstrating that the explicit formulation method has a problem.  Chapter 5. Simulation Comparison of Different Formulation Methods  Time (s)  Figure 5.17 Angular displacements of the two joints (upper-joint 1; lower-joint 2)  0.02  Time (s)  Figure 5.18 Elastic displacement of the flexible link tip  127  Chapter 5. Simulation Comparison of Different Formulation Methods  2  3 Time (s)  Figure 5.19 Total energy  0.06  —  i  —  1  1  -1.00 - 1.02 0.05  ;  E 0.04 z  •  <B 0.03  LU c co  C75 0.02  0.01  /u  2  3 Time (s)  Figure 5.20 Strain energy  128  Chapter 5. Simulation Comparison of Different Formulation Methods  129  5.5 Summary A detailed Order N formulation has been developed based on the approach of Keat [20]. The newly developed explicit formulation was programmed and verified using two conservative rigid systems. The simulation of a rigid-flexible system was compared between the joint coordinate method and the Order N method. The results show that the joint coordinate method is more accurate and efficient. The Order N method is an explicit formulation method requiring very small time steps in a flexible MBS to ensure numerical stability. Moreover, the error caused by inversion of the equivalent inertia matrix of each body is enlarged and propagated to the equivalent inertia matrices of the other bodies. The resulting errors may make the calculations of the small elastic deformations unreliable, especially for large-scale systems. A simulation comparison was also made between the joint coordinate method and the absolute coordinate method incorporating coordinate partitioning through the commercial software package ADAMS. The frequencies of the responses were the same, but the magnitudes of the joint displacements and the elastic displacements were different. The most probable explanation is the simplified modelling of the flexible bodies within ADAMS. A simulation using the Order N formulation displayed chaotic behavior for a nonchaotic system.  Chapter 6  Geometrical Nonlinearities 6.1 Preliminary Remarks In many mechanical system applications, geometric nonlinearities may have a negligible effect on the dynamic response of the system, while in some other applications, where lightweight mechanical systems operate at high speeds and under heavy loads, the effect of geometric nonlinearities can not be ignored [54]. The issue of geometric stiffening has been a topic of many recent publications dealing with deformable MBS [53]. The controversy over the nature of geometric stiffening, the debate on the correct approach to model it, and the seeming incongruity between existing methods has motivated some review papers [53,54]. These papers strive to provide a better understanding of how geometric stiffening is incorporated and what approximations are made in its derivation. Although simulation results have been compared for different methods, no conclusions have been drawn as to which method or approximation is the most reasonable. In this section, the two most common geometric nonlinearities are discussed. The two nonlinear effects are discussed in order to understand them better. Simulation results including the effects are compared with experiment results from the test rig.  6.2 The Two Geometric Nonlinear Effects Geometric nonlinearities arise when deformations are large enough to significantly alter the way load is applied to or the way load is resisted by the structure. Two effects have been discussed. The first is due to axial forces such as externally applied forces or inertial and centrifugal forces due to large rotational and translational movements during high speed operation. The other one is a foreshortening effect due to a transverse displacement which 130  Chapter 6. Geometrical Nonlinear  131  gives rise to an axial deformation. The most common approach to incorporate geometric nonlinearity is to retain nonlinear higher-order terms in the strain-displacement relationship [53,55,56]. Green-Lagrange strain defines a strain measure as:  du 1 fdu"\2 'AO \dx) + K.dx)+ 2  2  dv 1 (du^  = — + -  dz 2  z  Y  =  / x y  1  fdv^  T  =  l^y)  (6.3)  +  dz  2dy  2dx  1 dw  1 du  1 'du du dv dv  1  1 <9v  (6.2)  i  ifdudu^dvdv^dw dw^ 2 dx dy dx dy dx dy  2  x z  y  2  1  dw 1  £,  (6.1)  dx 2 dz  1 dv 2 dz  1  Idw 2 dy  1  1-  2 \dx  dz  dx dz  \ du du (  1-  2 ydy dz  (6.4)  dw dw^  (6.5)  dx dz  dv dv  dw dw  dy dz  dy dz  + — — +  (6.6)  where u, v, w are the elastic displacements in three orthogonal directions. The initial terms in equations (6.1) ~ (6.6) are the customary engineering definitions of normal and shear strain. The added terms, in parentheses, become significant if displacement gradients are not small. Compared with the initial terms, the terms in parentheses for e ,e ,y y  z  xy  ,y ,y xz  yz  are higher-  order of magnitude and can be neglected in beam applications. Usually — is several orders  dx  Chapter 6. Geometrical Nonlinear  132  dv ^ dw dv dw smaller than ^ dx- or dx but may be the same order as (—) dx' and (—) dx . Equations (6.1) ~ 2  2  (6.6) can be expressed as  £ = Du +  D,q N q  where e = [e  e  x  e  y  0  0  „d 2— dy  0  0  0  d_ dy d_  d  „d 2— dz  dz  f  y ]  T  (6.8)  yz  (6.9)  0  dx  d  0  dx d dy  d  0  t  y„  xy  0  'dx  D=  f  y  z  (6.7)  T  f  dz  £>. = [l 0 0 0 0 0]  u =[u  v  f  w] =[s? S T  dN v  C,B\,B ,qf 2  (6.10)  T  dN,  dx j v  d  x  J  2  Sl1cB B q =[Nl x  dN V *  2  f  dN  3  )  v  dx j  are defined in Chapter 2.  Thus the strain energy can be expressed as  N  T  N Jq T  2  3  f  = Nq  f  (6.11)  (6.12)  Chapter 6. Geometrical Nonlinear  U=-je Ee  133  dV  T  = \l t\{DN  + D,q N^E{DN  T  +  T  f  = \q][\{DN)  E{DN)dV  T  \q]\[{DN)'  l(DN) E(DN)dV T  dN , T  Sri  T  T  t  f  „  " dN,  N  = \q B B C T  T  f  T  2  T  x  f  E(DN)^Vq  }  If  tiVq,  + ^[\( *<lf){ *<lf) dv]qf N  f  + (D.q N, f  T  f  T  t  If  E(D q N ) t  D q N,)q dV  (K, + K  T  ga  +  K^CB^q,  (6.13)  (6.14)  where e is Young's modulus. The first term in equation (6.13) comes from the linear components of the straindisplacement relationship and will produce the linear stiffness matrix. The second and third terms are the strain energies caused by axial forces and foreshortening effects respectively.  Chapter 6. Geometrical Nonlinear  134  6.2.1 Nonlinear Effect Induced by Axial Forces Consider a one-dimensional plane beam undergoing transverse vibrations while subject to an axial force, as shown in Figure 6.1. The total strain energy includes not only the strain energy due to axial extension and transverse bending but also the strain energy due to axial forces. Large transverse bending results in shortening in the direction of the axis. If there is no axial forces, no extra axial strains are produced. Large axial forces resist the shortening in the direction of the axis and cause axial strains. The strain energy produced by the axial force P is equal to the negative of the work done by the axial force at the axial displacement.  +P  dx Figure 6.1 Beam vibration under axial forces  From Figure 6.1,  dA =dx-cos9dx  = —sm 0dx~ —0 dx = dx 2 2 2{dx )  (6.15)  (6.16)  Note that P = eA— , thus dx  strain energy = U „ =-PA a  eA , du( dw  = —J 2  dxydx,  dx  (6.17)  Chapter 6. Geometrical Nonlinear  135  which is the same form as the second term of equation (6.13) when applied to one dimensional problems.  6.2.2 Foreshortening Nonlinear Effect  A  -  dx Figure 6.2 Foreshortening effect  Consider a one-dimensional plane beam simply supported at both ends, as shown in Figure 6.2. When it vibrates, axial strain is induced since the two ends are fixed. Let a small lateral displacement w = w(x) take place. Each differential length dx is increased to a new length ds because the distance between supports is not allowed to change. From Figure 6.2,  ds = dx.1 + \dx j  'dw^  if c5 =  (6.18)  is very small, then  (l 5 )  = l+-c5-- 5 2 8  1/2  2  +  (  (  ds ~ dx\ 1 + v  Ifdw^  —  Adx,  +  (6.19)  (6.20)  Chapter 6. Geometrical Nonlinear  136  The strain due to the foreshortening is  ds-dx  1 'dw*  dx  (6.21)  ydx j  Thus the strain energy due to the foreshortening is  l 9 &i Ugf =—J eAe dx = — J e  dx  f  ydx j  (6.22)  6.3 Comparison of Simulation and Experiment Results The two different types of geometric nonlinearities in MBS are the third-order nonlinear term and the fourth-order nonlinear term in the strain energy equation 6.13. To date, only simulation results have been used in discussions about which formulation or approximation should be used. This section gives a comparison between the simulation results of the two derived formulations and experimental results using the test rig. The nonlinear stiffness matrix K due to axial forces is expressed as ga  dS  r  dx  The terms  If  B'BjC  dx  1  dS  T 2  2  dx  dx  |  c95  r 3  dx  d§  dS  ( 3  dx  2  l  2  dx  q B fi, C and CB B q f  T  2  f  dS  2  dx  |  dSj dx  dS ^ }  dx  CB,B q ^2  f  dx  IdV  (6.23)  — - express the strain induced by axial forces. dx  Usually, this strain term is very small compared with other deformation terms. For simplicity, the average strain  u  2  —  I  stiffness matrix. Thus K  w,  (with axial end displacements u^andi^) is used to develop the is simplified as  Chapter 6. Geometrical Nonlinear  137  ^ d S dS =y(" -"i)J dx dx  dS dS dV dx d. T  r  s a  K  2 +  2  K  3  2  3  (6.24)  will be given in Appendix C.  The nonlinear stiffness matrix K  which involves the foreshortening effect is defined as  a  dS  2  dSf  dx  dx  ds^)  T  dx  [A, A  2  -A,  2  d$l  dS  dS*  dx  dx  dx  2  d§ ) 3  (6.25)  dx  where [A, A ••• A ] = q B B C^ T  2  12  T  s  T  2  x  Kg will be given in Appendix D. From equation 6.14,  = B BjC {K T  + K + K^CB&q,  T  2  f  (6.26)  gv  CB B q  q)BlB[C  T  where  +Q  ga  x  2  f  Q = |  (6.27)  gv  q]B\B[C  CB B q  T  1  2  f  in which the derivatives of K with respect to the elastic coordinates are neglected. ga  Then the equation (2.144) becomes  B MBq = B (Q T  T  V  "o  c  s  0  Q -MBq) g  (6.28)  0 "V"  0  where Q' = - 0 0  0  g  0 0  +Q +Q +Q +  'fg_  K  (6.29) i s .  Chapter 6. Geometrical Nonlinear  138  (6.30) In order to compare with experiment results, simulations were performed under the same conditions as the experiments conducted in Chapter 4. Figure 6.3 ~ 6.6 show the simulation differences in bending strain between the formulations with only linear stiffness and including the geometric nonlinear due to axial forces or the geometric nonlinearity due to the foreshortening effect and axial forces. Figure 6.7 ~ 6.10 show the comparisons between the experimental results and the simulations including either the geometrical nonlinearity due to axial forces or the geometrical nonlinearity due to the foreshortening effect and axial forces. The responses were investigated using two different initial conditions 6 = 30° and 6 = 45°. 2  2  -4  4  X 10  — linear stiffness • • • nonlinear including only axial forces  3  _2i 0  i  i  2  4  i  i  I  6  8  10  Time (s) Figure 6.3 Comparison between nonlinear and linear under initial condition 0 = 30 2  Chapter 6. Geometrical Nonlinear  x 10 stiffn©ss nonlinear including axial forces and foreshortening  ' linGcir  E c CO  -1  -2  4  6 Time (s)  10  Figure 6.4 Comparison between nonlinear and linear under initial condition d  2  Chapter 6. Geometrical Nonlinear  Time (s)  Figure 6.6 Comparison between nonlinear and linear under initial condition 0 = 45 2  4  x 10 simulation including only axial forces  3 2 1 c  'CO  CO  0 -1 -2 -3  2  4 Time (s)  6  10  Figure 6.7 Comparison of experiment and simulation including only axial forces d = 2  Chapter 6. Geometrical Nonlinear  -2  0  2  4  6  8  10  Time (s)  Figure 6.8 Comparison of experiment and simulation including both nonlinearities 9  x10  -2  0  2  4  Time (s)  6  8  10  Figure 6.9 Comparison of experiment and simulation including only axial forces 6-  Chapter 6. Geometrical Nonlinear  X 1 0 "  4  L  4  -2  142  0  2  4  6  8  10  Time (s)  Figure 6.10 Comparison of experiment and simulation including both nonlinearities 6 = 45° 2  The above simulation results show that there are large differences between the results using the two different geometrical nonlinearities. The comparisons of experiment and simulation results demonstrate that the strain response is better described by the geometric nonlinearity due to axial forces than due to foreshortening. The geometric nonlinearity due to foreshortening overestimates the nonlinear effect. The reason is that in this application both ends of the flexible link are not fixed but are free to move. However, if a large tip mass were used that limited the movement of the flexible link, then the foreshortening effect might have been significant. The selection of which kind formulation of geometrical nonlinear depends on the applications themselves. Moreover, the nonlinear simulation results including axial forces only were not significantly better than the linear simulation results compared with the experimental results in the case of the test rig.  Chapter 6. Geometrical Nonlinear  143  6.4 Summary A detailed formulation of the geometric nonlinearities has been derived using norilinear strain-displacement relations. The resulting strain energy includes a second-order term, a third-order term and a fourth-order term. The second-order term yields the linear stiffness term. The third-order and fourth-order terms corresponds with the two different nonlinear effects; the nonlinear effect caused by axial forces and by the foreshortening effect. The simulation and experimental results using the test rig under two different initial conditions show that the fourth-order nonlinear term in the strain energy expression overestimates the nonlinearity. The experimental results show that there is little foreshortening effect because the ends of the flexible link are not constrained. The nonlinear simulation results including axial forces only were not significantly better than the linear simulation results compared with the experimental results.  Chapter 7  Summary and Conclusions  A general and efficient implicit formulation for arbitrary deformable MBS has been developed. The derivation of the new formulation is based on the joint coordinate method through defining a new topological matrix. The new formulation has the following advantages: (i) The formulation is implicit so can ensure stability for any stiff systems; (iii) The formulation has a minimum number of variables and equations for high computational efficiency; (iii) The formulation is constructed independent of the numerical solution technique so that it can have high accuracy and numerical stability; (iv) It can be easily extended from tree-configured systems to closed-loop systems; (v) The dynamic equations of the individual bodies are formed independently and the global velocity transformation matrix and its derivative matrix can be constructed in different directions so that the equations can be easily solved using parallel computation for the simulation of large-scale systems. The new formulation has been incorporated into a general time-domain simulation program. The total energy has been verified for a 3-D two-link manipulator with one link rigid and the other one flexible and a three-link manipulator. The importance of strain energy within the total energy balance has been stressed. The strain energy must be several orders of magnitude higher than the absolute error of total energy. Otherwise, the small elastic displacements can not be trusted. The finite element analysis part of the program was verified by comparing results while using one element andfiveelement models of the flexible link. An eight-link: treelike manipulator connected by six revolute joints, one translational joint and one spherical joint was employed as a test example to demonstrate that the formulation and the code can deal with tree-configured MBS systems. A modal representation was also proved to be effective in reducing the number of elastic coordinates. 144  Chapter 7. Summary and Conclusions  145  The effect of small elastic deformations on the joint motions was found to be significant. Not only the amplitudes of the joint motions but also the frequencies were different for a two-link manipulator when one of the links was considered as rigid rather than flexible. Nonlinear coupling was identified between the joint motion and the elastic motion. A 3-D two-link test rig with one link rigid and the other one flexible was designed to verify the formulation. The experimental results of the joint motion and small elastic displacements were found to be in good agreement with the simulation results. The results were compared in the time domain and in the frequency domain. The frequency of the flexible link vibration mode was found to change with the amplitude of the joint motions. For the test rig, joint flexibility needed to be considered as well as structural flexibility. The experiment results were also compared with the simulation results achieved using the commercial software package, AD AMIS. A comparison of the results showed that the magnitudes of the joint displacements obtained by ADAMS were larger than that of the experiments. The comparison of the strain signal between ADAMS and the experiments also indicated that the strain signal was larger than the experimental results. An Order N formulation was derived following Keat [20] for further comparison. The Order N formulation was very efficient in a single time step but was much slower when simulating a given time interval because of the much smaller time steps required. The Order N formulation required sequential inversion of a numbers of inertial matrices. Numerical errors from one inversion may be propagated and enlarged by subsequent inversions. As a results, reducing the time step may not be helpful for increasing the accuracy. In deformable MBS, the coupling between small elastic displacements and large rotation and translation displacements may induce matrices to be ill-conditioning resulting in chaotic behavior through error propagation. A simulation comparison using the Order N method using two slightly different initial conditions demonstrated chaotic responses for a nonchaotic system.  Chapter 7. Summary and Conclusions  146  The developed joint coordinate formulation was also compared with the absolute coordinate formulation incorporated coordinate partitioning through the commercial software ADAMS. The simulation results showed that the frequencies were same but the magnitudes of the joint motions and the elastic displacements were different. The results achieved by ADAMS were larger. The reason for this is probably due to the approximation method that ADAMS uses to model flexible bodies. ADAMS models flexible bodies by discretizing the flexible bodies into lumped inertias and uses force elements to represent the elastic stiffness and the damping forces. The force elements are applied directly to the lumped inertias. The numerical instability caused by coordinate partitioning can be seen from the simulation results of the elastic displacements. Geometric nonlinear behavior was investigated to demonstrate the nonlinear capabilities of the formulation. Two different nonlinear effects produced by axial forces and foreshortening were examined. The simulation results were compared with experiment results and it was found that the formulation that includes foreshortening overestimated the nonlinear effect for this test rig. The reason is that the flexible link is assumed to be fixed at each end when calculating foreshortening, which was not the case in the test rig. The controversy over which one should be neglected is not a mathematic problem. It depends on what kind of application is being investigated. The main contributions of this thesis can be summarized as follows: (i) A general and implicit formulation of deformable MBS based on the joint coordinate method has been successfully developed and implemented into a program. The formulation is consistent across rigid and flexible bodies. No approximations are required for modelling flexible bodies beyond those used in finite element applications. The newly-developed formulation has some advantages compared with other methods. (ii) The numerical and experimental verification of the newly-developed formulation for the 3D rigid-flexible MBS has been conducted. The importance of strain energy was stressed in the  Chapter 7. Summary and Conclusions  147  numerical validation. Good agreement in the time domain and the frequency domain between the simulation and experimental results supports the correctness of the formulation. (iii) A simulation comparison has been conducted between the joint coordinate method, the Order N method and the absolute coordinate method (ADAMS) to investigate accuracy and efficiency. The results show that the Order N method is neither efficient nor accurate for deformable MBS due to the very small time step required and the propagation of numerical errors. The study demonstrates that Order N method may induce chaotic behavior for an unchaotic system. The comparison also shows that numerical instabihty and inefficiency may occur in ADAMS. Also, ADAMS has larger joint and elastic motions due to the approximate modelling method of flexible bodies. (iv) Geometrical nonlinear effects for deformable MBS have been investigated. Two nonlinear effects caused by axial forces and the foreshortening effect were distinguished and compared to each others and to the experiment results. The results show that which effect should be considered depends on the application being investigated. Further research should concentrate on the following aspects: (i) To conduct more experimental work to verify the two geometrical nonlinear effects. (ii) To investigate an efficient method to solve the constraint drift problem for closed-loop MBS. (iii) To study effective and reliable control methods for the control of flexible MBS.  Bibliography [1]  Wu, Sijun and Cetinkunt, Sabri, 1995, 'Techniques in Discrete-Time Position Control of Flexible One Arm Robots," Control and Dynamic Systems, Vol. 70, pp. 291-351.  [2]  Cook, R. D., Malkus, D. S., and Plesha, M . E., 1989, Concepts and Applications of Finite Element Analysis, John Wiley & Sons.  [3]  Shabana, A. A., 1989, Dynamics of Multibody Systems, John Wiley & Sons, New York.  [4]  Shabana, A. A., 1991, "Constrained Motion of Deformable Bodies," International Journal for Numerical Methods in Engineering, Vol. 32, pp. 1813-1831.  [5]  Bayo, E. and Ledesma R., 1996, " Augmented Lagrangian and Mass-Orthogonal Projection Methods for Constrained Multibody Dynamics," Nonlinear Dynamics, Vol. 9, pp. 113-130.  [6]  Andrzejewski, Th., Bock, H. G., Eich, E. and Von Schwerin, R., 1993, " Recent Advances in the Numerical Integration of Multibody Systems," Advanced Multibody System Dynamics, edited by Schiehlen, W., Kluwer Academic Publishers, Dordrecht, pp. 127-151.  [7]  Petzold, L. R., 1993," Computational Challenges in Mechanical Systems Simulation," Computer-Aided Analysis of Rigid and Flexible Mechanical Systems, edited by Pereira, M. F. O. S., and Ambrosio, J. A. C , Kluwer Academic Publishers, Dordrecht, pp. 483499.  [8]  Ostermeyer, G., 1989, " On Baumgarte Stabilization for Differential Algebraic Equation," Real-Time Integration Methods for Springer-Verlag, pp. 193-207.  mechanical System Simulation,  Bibliography  [9]  149  Bae, D. and Yang, S., 1989, " A Stabilization Method for Kinematic and Kinetic Constraint Equations," Real-Time Integration Methods for  mechanical System  Simulation, Springer-Verlag, pp. 209-231. [10] Bayo, E., Jalon, J. G. D., Avello, Alejo and Cuadrado, Javier, 1991, " An Efficient Computational Method for Real Time Multibody Dynamic Simulation in Fully Cartesian Coordinates," Computer Methods in Applied Mechanics and Engineering, Vol. 92, pp. 377-395. [11] Shabana, A. A., Hwang, Y. L., and Wehage, R. A., 1992, " Projection Methods in Flexible Multibody Dynamics, Part I: Kinematics, " International Journal for Numerical Methods in Engineering, Vol. 35, No. 10, pp. 1927-1939. [12] Shabana, A. A., Hwang, Y. L., and Wehage, R. A., 1992, " Projection Methods in Flexible Multibody Dynamics, Part II: Dynamics and Recursive Projection Methods," International Journal for Numerical Methods in Engineering, Vol. 35, No. 10, pp. 1941-1966. [13] Lieh, Junghsen, 1994, " Separated-Form Equations of Motion of Controlled Flexible Multibody Systems," Transactions of the ASME, Journal of Dyanmic Systems, Measurement, and Control, Vol. 116, pp.702-712. [14] Haug, E. J., Yen, J., 1989, " Generalized Coordinate Partitioning Methods for Numerical Integration of Differential-Algebraic Equations of Dynamics," Real-Time Integration Methods for mechanical System Simulation, Springer-Verlag, pp.97-126. [15] Kane, T. R. and Levinson, D. A., 1980, " Formulation of equations of Motion for Complex Spacecraft," Journal of Guidance and Control, Vol. 3,No. 2, pp. 99-112.  Bibliography  150  [16] Huston, R. L., 1991," Computer Methods in Flexible Multibody Dynamics," International Journal for Numerical Methods in Engineering, Vol. 32, pp. 1657-1668. [17] Huston, R. L. and Wang, Y., 1993, " Flexibility Effects in Multibody Systems," Computer-Aided Analysis of Rigid and Flexible Mechanical Systems, edited by Pereira, M. F. O. S. and Ambrosio, J. A. C , Kluwer Academic Publishers, Dordrecht, pp.351376. [18] Singh, R. P. et al, 1985, " Dynamics of Flexible Bodies in Tree Topology-A ComputerOriented Approach," Journal of Guidance, Control, and Dynamics, Vol. 8, No. 5, pp.584-590. [19] Surdilovic, D. and Vukobratovic, M., 1996, " One Method for Efficient Dynamic Modeling of Flexible Manipulators," Mech. Mach. Theory, Vol. 31, No. 3, pp. 297-315. [20] Keat, J. E., 1990, " Multibody System Order n Dynamics Formulation Based on Velocity Transform Method ," Journal of Guidance, Control, and Dynamics, Vol. 13, No. 2, pp. 207-212. [21] Rosenthal, D. E., 1990, " An Order n Formulation for Robotic Systems," The Journal of the Astronautical Sciences, Vol. 38, No. 4, pp. 511-529. [22] Lilly, K. W., 1993, Efficient Dynamic Simulation of Robotic Mechanisms, Kluwer Academic Publishers, Boston. [23] Stelzle, W., Kecskemethy, A. and Hiller, M., 1995, " A comparative Study of Recursive Methods," Archive of Applied Mechanics, Vol. 66, pp.9-19.  Bibliography  151  [24] Pradhan, S., Modi, V. J. and Misra, A. K., 1997, " Order N Formulation for Flexible Multibody Systems in Tree Topology: Lagrangian Approach," Journal of Guidance, Control, and Dynamics, Vol., 20, No. 4, pp. 665-672. [25] Van Woerkom, P. T. L. M. and Boer, A. D., 1995, " Development and Validation of a Linear Recursive " Order-N " Algorithm for the Simulation of Flexible Space Manipulator Dynamics," Acta Astronautic, Vol. 35, No. 2/3, pp. 175-185. [26] Kim, Sung-Soo and Haug, E. J., 1988, " A Recursive Formulation for Flexible Multibody Dynamics, Part I: Open-Loop Systems," Computer Methods in Applied Mechanics and Engineering, Vol. 71, pp. 293-314. [27] Changizi, K. and Shabana, A. A., 1988, " A Recursive Formulation for the Dynamic Analysis of Open Loop deformable  Multibody Systems," Journal of Applied  Mechanics, Vol. 55, pp. 678-693. [28] Caron, Mathieu, 1996, " Planar Dynamics and Control of Space-Based Flexible Manipulators with Slewing and Deployable Links ," MA. Sc. Thesis, University of British Columbia, Canada. [29] Baker, Gregory L. and Gollub, Jerry P.,  1996, Chaotic Dynamics, Cambridge  University Press. [30] Edward Ott, 1993, Chaos in Dynamical Systems, Cambridge University Press. [31] Basband Neil S., 1990, Chaotic Dynamics of Nonlinear Systems, John Wiley & Sons, New York. [32] Schwertassek, R., 1993, " Reduction of Multibody Simulation Time by Appropriate Formulation of Dynamic System Equations," Computer-Aided Analysis of Rigid and  Bibliography  152  Flexible Mechanical Systems, edited by Pereira, M. F. O. S., and Ambrosio, J. A. C , Kluwer Academic Publishers, Dordrecht, pp. 447-482. [33] Jerkovsky, W., 1978, " The Structure of Multibody Dynamics Equations ," Journal of Guidance and Control, Vol. 1, pp. 173-182. [34] Kim, S. S. and Vanderploeg, M. J., 1986, " A General and Efficient Method for Dynamic Analysis of Mechanical Syatems Using Velocity Transformations," ASME Journal of Mechanisms, Transmissions, and Automation in Design, Vol. 108, pp. 176182. [35] Pankiewicz, E., 1993, " NUBEMM A Special Multibody System As A Part of A Modern Simulation Concept in the Automotive Industry," Multibody Computer Codes in Vehicle System Dynamics, edited by Kortum,W. and Sharp, R.S., Swets & Zeitlinger, pp. 99-103. [36] Nikravesh, P. and Ambrosio, J., 1991, " Systematic Construction of Equations of Motion for Rigid-Flexible Multibody Systems Containing Open and Closed Kinematic Loops," International Journal for Numerical Methods in Engineering, Vol. 32, pp. 1749-1766. [37] Ambrosio, J. A. C , 1996, " Dynamic of Structures Undergoing Gross Motion and Nonlinear Deformations: A Multibody Approach," Computers & Structures, Vol. 59, No. 6, pp. 1001-1012. [38] Pereira, M. S. and Proenca, P. L., 1991, " Dynamic Analysis of Spatial Flexible Multibody Systems Using Joint Coordinates," International Journal for Numerical Methods in Engineering, Vol. 32, pp. 1799-1812.  Bibliography  153  [39] Roberson, R. E. and Schwertassek, R., 1988, Dynamics of Multibody Systems, Springer-Verlag. [40] Schiehlen, W., 1990, Multibody Systems Handbook, Springer-Verlag. [41] Hoffmann, R. et al, 1990, " Finite Element Analysis of Occupant Restraint System Interaction with PAM-CRASH," SAE paper 902325. [42] Lupker, H. A., 1993, " Comnined Finite Element and Multibody Techniques for Vehicle Occupant Safety Studies," Topics in Applied Mechanics, edited by Dijksman, J. F. and Nieuwstadt, F. T. M., Kluwer Academic Publishers, pp. 277-284. [43] Ho, J. Y. L. and Herber, D. R., 1985, " Development of Dynamics and Control Simulation of Large Flexible Space Systems," Journal of Guidance, Control, and Dynamics, Vol., 8, No. #, pp. 374-383. [44] Meijaard, J. P., 1996, " Validation of Flexible Beam Elements in Dynamics Programs," Nonlinear Dynamics, Vol. 9, pp21-36. [45] Sinha, S. C , Benner, J. W. and Wiens, G. J., 1992, " Experimental Verification of Component Mode Techniques for A Flexible Multibody System," ASME Dynamics of Flexible Multibody Systems: Theory and Experiment, AMD-Vol. 141/Dsc-Vol. 37, pp.65-75. [46] Whalen, Mark and Sommer III H. J., 1990, " Modal Analysis and Moda Shape Selection for Modeling an Experimental Two-Link Flexible Manipulator," Robotics Research, DSC-Vol. 26, pp. 47-52. [47] Fraser, Anthony R. and Daniel, Ron W., 1991, Perturbation Techniques for Flexible Manipulators, Kluwer Academic Publishers, Boston.  Bibliography  154  [48] Figueiredo, J. and Costa, J.M.G. Sa da, 1996, " Elastic-Link Manipulators: Modelling, Simulation and Experiments," International Journal of Robotics and Automation, Vol. 11, No. 1, pp. 13-21. [49] Kirk, C. L. and Junkins, J. L. (Editors), 1990, Dynamics of Flexible Structures in Space, Springer-Verlag, Berlin. [50] Ladkany, S. G., 1991, " The Dynamic Response of a Flexible Three-Link Robot Using Strain Gages and Lagrange Polynomials," Fourth World Conference on Robotic Research, Sept. 17-19, Pittsburgh. [51] Juang, Jer-Nan, Horta, L.G. and Robertshaw, H.H., 1986, " A Slewing Control Experiment for Flexible Structures," Journal of Guidance, Control, and Dynamics, Vol. 9, No. 5, pp. 599-607. [52] Hanagud, S. and Sarkar, S., 1988, " Problem of the Dynamics of a Cantilever Beam Attached to a Moving Base," Journal of Guidance, Control, and Dynamics, Vol. 12, No. 3, pp. 438-441. [53] Sharf, I., 1995, " Geometric Stiffening in Multibody Dynamics Formulations," Journal of Guidance, Control, and Dynamics, Vol. 18, No. 4, pp. 882-890. [54] Mayo, J., Dominguez, J. and Shabana, A. A., 1995, " Geometrically Nonlinear Formulations of Beams in Flexible Multibody Dynamics," Journal of Vibration and Acoustics, Vol. 117, pp. 501-509. [55] Bakr, E. M., 1996,  " Dynamic Analysis of Geometrically Nonlinear Robot  Manipulators," Nonlinear Dynamics, Vol. 11, pp. 329-346.  Bibliography  155  [56] Noej, A. T., et al, 1997, " Formulation of Geometrically Nonlinear Problem in the Spatial Reference Frame," International Journal for  Numberical Methods in  Engineering, Vol. 40, pp. 1263-1280. [57] Torby, B.J., 1984, Advanced Dynamics for Engineers, Holt, Rinehart and Winston, New York. [58] Brenan, K.E., Campbell, S.L. and Petzold, L.R., 1989, Numerical Solution of InitialValue Problems in Differential-Algebraic Equations, North-Holland. [59] Cook, Robert D., 1994, Finite Element Modeling for Stress Analysis, John Wiley & Sons, Inc, New York. [60] Spong, M. W. and Vidyasagar, M., 1989, Robot Dynamics and Control, John Wiley & Sons. [61] Grewal, Anant Kiran Singh, 1994, " A Study of Flexible Space Structures: Dynamics and Control," Ph.D Thesis, University of British Columbia, Canada. [62] Thompson, B. S. and Sung, C. K., 1984, " A Variational Formulation for the Nonlinear Finite Element Analysis of Flexible Links: Theory, Implementation, and Experimental Results," Transaction oftheASME, Vol. 106, pp. 482-488. [63] Johnke, M. , Popp, K. and Dirr, B., 1993, " Approximate Analysis of Flexible Parts in Multibody Systems Using the Finite Element Method," Advanced Multibody System Dynamics, edited by Schiehlen, W., Kluwer Academic Publishers, Dordrecht, pp. 237256.  Bibliography  156  [64] Chun, H. M., Turner, J. D. and Frisch, H. P., 1991, " A Recursive Order-N Formulation for DISCOS with Topological Loops and Intermittent Surface Contacts," Advances in the Astronautical Sciences, Vol. 76, pp. 493-511. [65] Chun, H. M . , Turner, J. D. and Lupi, V., 1991, " Distributed Parameter Multibody Dynamics Modeling," Advances in the Astronautical Sciences, Vol. 76, pp. 513-526. [66] Chun, H. M., Turner, J. D. and Frisch, H. P., 1989, " Experimental Validation of Order (N) DISCOS," AAS 89-457, pp. 1341-1358. [67] Khulief, Y. A. and Shabana, A. A., 1986, " Dynamics of Multibody Systems; with Variable Kinematic Structure," ASME Journal of Mechanisms, Transmissions, and Automation in design, Vol. 108, pp. 167-175. [68] Chang, C. W. and Shabana, A. A., 1990, " Spatial Dynamics of Deformable Multibody Systems with Variable Kinematic Structure," Journal of Mechanical Design, Vol. 112, pp.153-167. [69] Wu, S. C. and Haug, E. J., 1990, " A Substructure Technique for Dynamics of Flexible Mechanical Systems with Contact-Impact," Transactions of the ASME, Vol. 112, pp. 390-398. [70] Khulief, Y. A. and Shabana, A. A., 1987, " A Continuos Force Model for the Imapct Analysis of Flexible Multibody Systems," Mech. Mach. Theory, Vol. 22, pp. 213-224. [71] Khulief, Y. A. and Shabana, A. A., 1986, " Dynamic Analysis of Constrained System of Rigid and Flexible Bodies with Intermitten Motion," Transactions of the ASME, Vol. 108, pp. 38-45.  Bibliography  157  [72] Kane, T. R. and Ryan, R. R., 1987, " Dynamics of a Cantilever Beam Attached to a Moving Base," Journal of Guidance, Control, and Dynamics, Vol. 10, pp. 139-151. [73] Sorge, K., Bremer, H. and Pfeiffer, R, 1993, " Multi-Body Systems with Rigid-Elastic Subsystems," Advanced Multibody System Dynamics, edited by Schiehlen, W., Kluwer Academic Publishers, Dordrecht, pp. 195-215. [74] Oakley, C. M . and Robert H. Cannon, Jr., 1990, " Anatomy of An Experimental TwoLink Flexible Manipulator Under End-Point Control," Proceedings of the 29th Conference on Decision and Control, pp.507-513. [75] De Luca, A., Lanari, L., Lucibello, P., Panzieri, S. and Ulivi, G., 1990, " Control Experiments on a Two-Link Robot with a Flexible Forearm," Proceedings of the 29th Conference on Decision and Control, pp. 520-527. [76] Ballhaus, W. L. and Rock, S. M., 1992, " End-Point Control of a Two-link Flexible Robotic Manipulator with a Mini-Manipulator: Initial Experiments," pp.2510-2514.  ACC/TP13,  Appendix A  The Orientation Transformation Matrices  The transformation matrices A' ,G',G',G'  of body i used in equation (2.17),(2.18),(2.103)  are expressed by the generalized orientation vector (Euler angles) q>' = [0 a y/J as T  folio wings respectively:  cos 0 cos y/ — sin 0 cos 0 sin yr• cos (j) sin yr - sin 0 cos t9 cos yr  A = sin 0 cos \ff + cos 0 cos 0 sin yr• sin 0 sin i/A + cos 0 cos 6 cos sin 6 sin y/  sin 0 cos  0  cos 0  sin 0 sin 0  G' = 0  sin 0  - sin 0 cos 0  1  0  cos0  sin 0 sin yi  cos I/A  0  - sin yt  0  G' = sin 0 cos i/f  cos0  .  0  sin 0 sin 0  - cos 0 sin 0 (A.l) cos 0  (A.2)  (A.3)  1  - cos 0 sin 0 cos 0 cos 0 sin 0 G' =  sin0  cos 0 sin 0  sin 0 sin 0  sin 0  - cos 0  (A.4)  0 0  158  Appendix B The Invariant Matrices in Finite Element Method The shape function S and other invariants A>7, N^N^, N'{ , SS^ (p,q = 1,2,3), PS and IJ  IJ  K'} used in equation (2.48)~(2.52), equation (2.81) and equation (2.91) for a 3-D beam element j of body i are expressed as followings:  6(^-^)7?  0  0  l-3<f+2<f  0  0  l-3cf+2<f  0 ( l - 4 £ + 3cf)/g  0  (-l + 4§-3$ )/u  (§-25 + £ >  0  5  0  0  6(-§ + 5 >I  3<f-2<f  0  6 H +f)?  0  l  S  =  IJT  2  (-£ + 2<f -<f)/  2  3  0  3f  (-2£ + 3§ )/c;  0  (2^-3^ )/77  (-<f+cf)/  2  2  (B.l)  -2?  0 0  PS  ijT  =  B?BfC  iiT  0  0 0  0  12(S-£ )//  12(-^ ^ ) / / +  0  2  2(1-£) 0 0 2(1-41 + 3^) 0 0 2 ( l - 4 £ + 3<f) 0 0 0 0 0 12(4-^)// 0 0 2£ 0 0  i2(-^+r)// 2(-2£ + 3cf) 0  s 0 2(-24 + 3^ ) 2  159  (B.2)  Appendix B. The Invariant Matrices in Finite Element Method  m (jpSdvj^  g  m  0  0  ml  ~2  ~ 12  /V|.=(Jp^v)' =  0  20  0  ' e 2  0  0  0  6  3m  0  20  Ok  T  0  0  2 0  iQc 2  7V;{ =(Jp^v)" =  0 0  0  2  «  u  IQ 2  0  0  2  0  0  0  0  0  ' e  2  ' e 2  2  0  12  2  12  2  0  0  /a. 2  s  0 0  ' e,  12  12 ml  7m  ' e,  ml  20  3  20  0  0  2  0  2  0  3  ~20 0  0  0  2  12  ' e„ 2  0  2  2  2  12 ?  0  0  12  2  0  ' e, 2  0  2  2  0  T  2  20  0  0  " IT  m  0  2 7m  2 0  0 m/  '2,  0  0  30  0  m  12 ml  0  0  ~2  0  ' e,  30  ml 6  0  2  s  12  m  0  IT  0  -'a,  y  /n/  m  2  m  0  0  m  2 3m  0  0  T  0  y  0  lQ  2~  160  0  0  0  J \  0  0  (B.6) 2  / fi 2  2  0  12  2  12  12 t  0  Appendix B. The Invariant Matrices in Finite Element Method  m J Q£ 2 e/ 2 0  6/ / ?  5  6y  s  5 0 2  If. 10 0/ 2 6l l  12 m ~6 .91 2  0  If. 10 hf 10  hf 10  7m  01  0  2  20  1% 2 lQ  0 0 0 0  'e 2  0  ;  3 2  15 Qf 12  m 3 2  6/ / ~5~ * V  6V  2 0  5 0  5 0  2/ 12 Qf 12  V  2  0  2  0  21/  3  10 V 10  10  15 21 / <K 15 Qf 12 V 10  0  Jf  '/  21/  0  5 _6/y 5 0  5 0  Qf 12 Qf . 12  0  0  f  .91 2 0  0  symmetric  5 0  V 10  Qf 12  0  161  0 0  ?  10 V 10 0  If. 10 0  2  V 30  If 30 V 30  3  '/  3  30  ml  0  20  1%  0  10  2  VI 10 V 10  10 V 10 2  'e  0 0  2  0  2  3m  0  20  'a,  0  2  ?  6  0  15 2/ / 15  3  27/ 15  ml ~30  0  2  21/  10  'e 2  0  g  2 0  0  2 0  0 0  0  10 0  0  0  2 0  'e  0  0  2 0  0  ?  10 0  2  0 0 0  10 10 3m 20  0  _ ^ 2  0  _ ^ 2 0  0 0  • '? > ;  10  0 0 0 0  i 12 12 6 2  0 0 0 0  0  0  s  10 lQ 10 7m  0  10  0  0  0  ml  0  30 ?Qn 10  'e  0  n  20 2  0  12  <e,  0  0 0  2 0 12  0 0  ?  10 0  0  60  0  12  0  ' e, 60  60 ml  2  0 0  3  0 0  2  2  0  60  2  ~20 10  'e 2  0 0 0  2 0 10  0 0  0  10  0  2 0  i  0  ?  0  10 0  0  0  0  0  12  3  0  n 12  0  'X  12  (B.8)  Appendix B. The Invariant Matrices in Finite Element Method  162  = ss« = (Jps7s^v)' r  0  7  7m  0  ' e,  _mj_  2  20  20  / /.  0  0  0  0  0  0  0  0  0  0  0  0  3m  ml  20~  ~ 30~  2  10  0  2  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0 0  2  !£<.  0  0  2 0  10 0  10  12  0  0  0  0  0  0  0  0  0  0  0  0  0  0  12  ml  ~ To ' e,  6  ~20  2  2  2  10  < '„  2 0  2 0  10  12  10 7m  12  20  3  60 1% 60 ml 20~  2  10  Z /„  0 0 0 0 0 0 0  2  2 0  2 0  10 0  z /.  2  2  Jll.  0 0  2 0  10 3m  0  0  10  10 0  / e,  0  0  0  0  0  0  0  0  0  0  0  0  2 0  10 0  0 0  (B.9)  3  J%  12  10  60  1^1  12  10  60  10  0  12  10  0  0  12  SS^ =(jpS S dvJ T  2  2  2  0 o 0  13m 35 0  0  7/ a 2  0 0 0 0 0 0  20 0 11m/ 210 0 9m 70 0  3/ G,  0 0 0  symmetric  n  3 0  0  1%  0  20 0  0  —— 11 20  3/ e  0  0  0  0  0  2  c  2  0 0 0  20 0 13m/ 420  0 0  6 0  0  z G,  0 0  3  30  ml  1  0  105 0 13m/ 420 0  ^% 30 0 ml  0 0 0  35 0  0  7/ e  140  0  2  0  2  0  13m  0  ?  20 0 11m/ 210  0 0 0  3 0 1% 20  0 ml  2  s  0  105  (B.10)  Appendix B. The Invariant Matrices in Finite Element Method  SS^SSg 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  =  0  T  P  i  2  0  35 0  20 0  -iiz  20 0 11m/  0  0  {\ S S,dv)  0 13m  0  0  163  3 " 0  0 11m/  0  210 0  1%  0  0  0  0  0  0  0  ml  0  0  0  0  0  0  2  210 0 9m  20 0  70 0  20 0  105 0 13m/ 420 0  20 0 13m/  6 " 0  30 0  1%  ml  420  30  140  5  0 9m  3/ e„  0 13m/  70 0  20 0  420 0  3/ e, 2  20 0  s  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  2  0  0  0  20 0 13m/  0 2  6  0  30 0  1%  ml  30 0  0  7/ e„  140 0 11m/  35 0  20 0  210 0  20 0 11m/  2  - -3' «" 0  0 0 0 0 0  20 0  s  1%  210  0  2  420 0 13m  11%  0 0  1%  -tl  0  0  ml  2  20  105  0  (B.ll)  0 0  0  0  0  0  0  0  0  0  0  0 0 0  13m 35  >\  20  ml  2  11m/  —210 0  20 0  0'  0  0 0  symmetic  3 105 0  0  0  0  0  0  0  0  0  0  9m  3' e„  13m/  20  420  0  0  0  0  — 70  2  lL  13m 35  l'h  3  0  0  0  0  0  0  20  30  ml  1  13m/ 420 0  30 0  140 0  0  0  0  0  20 11m/ 210 0  m/ 20 0  z  105 0  (B.12)  Appendix B. The Invariant Matrices in Finite Element Method  164  ea  T  12c/  c  0 12f^  0  0  0  0  0  0  0  symmetric  Ip  P  6el  GJ I 0  4 ^  P  (  0  0  P  ea  ~T  0  0  12e/  ;  0  0  Zp  0  0  0  0  0  0 0  0  6eL 0  0  0  0 _GJ_  4ell  (  P 0 0 6^  12^,  0  0  0  P  r  0  6e/„  0 0  0 6e/  T ?  P 0  0 0  0  P 2ell  r  12e/  5  /p 0  12W, Ip  0  0  0  0  6eL  GJ  (B.13)  ~T P  0  P  0  where p , / , a are the mass density,length, and cross sectional area of the element j of body y  y  y  i. e' and G' are the modulus of elasticity and the modulus of rigidity of the element j of body J  J  i. Ql, Ql, /J, I , 1%. are defined in equation (2.39)~(2.43). J is the polar moment of inertia iJ  of the cross section.  Appendix C  The Nonlinear Stiffness Matrix Due to Axial Forces The geometric nonlinear stiffness matrix K' of the element j of body i due to axial forces is J  g a  expressed as  as ds  ds, ds  dx dx  dx dx  T  r  2 f  2  0  5 0  0  0  0  0  0 0  la 10 0 6a 5  0  0  0  0  0 0  dv  6a  0  0  3  0 la 10  6a 5 0  symmetric  la  2l a 2  0  10  15 2l a 2  0  0  0  0  0  0  0  0  0  6a 5 0 la 10 0  la  0 -'  3  (v  ?  )  15 0 la  "To  10 0 la  0 0  0  0  0  0  6a  T  0 0  2  0 0  0  30 0  fa __  30  165  0 0  0 la  "To  6a ~  0 la  2l a  To  15  0  0  (C.l)  2  2l a 2  15  Appendix D The Nonlinear Stiffness Matrix Due to Foreshortening The geometric nonlinear stiffness matrix  of the element j of body i due to foreshortening  is expressed as  ds  (d§l  dS* dS,  2 |  dx dx  0  k  0  ^32  - A  2  1  d§T dS,  ]  2  dS! dS. dx dx  dx dx  dV  iJ  22  &  0  symmetric ^43  ^44  ^"52  k 5j  0  **  ^62  ^63  k  0  0  0  -k  0  ~k  0  dx dx  A  0  0  0  [A,  42  ~^-32  —k  0  —k  ~k$2  —  ~ ^53  42  A. 4  —^43  42  0 -i  -*33  }2  kf,S  M  0  22  (D.l) 0  k2  0  -*«  0  —&  0  _£-  6  3  - ^ 4 4  0  -*«•  -/C54  0  km  ^113  ^114  ^-115  ^•116  0  0  k\22  ^123  *124  ^125  ^126  0  _  ^82  71  _  ^93  92  102  —  ^"104  ^113  _  ^114  ^123  ~^"124  ^103  ^ U 2  -k  _  _  *I22  where Halt .2 18a/ / *22 = — (A -A ) + — ( A - A 35 2  8  2  8  w ) ( A + A 6  1  )  2  ^2=^(A2-A )(A3-A )+^[(A -A )(A +A 8  6(2/ *42=  +/)/ 5p  9  3  9  6  \  1  +  3a/ / o — ( 2 _ 2 A  2  )-(A -A )(A +A  A  2  (2/ + / ) /  3  ( A  2  - A  8  ) ( A  4  - A  1  0  )  V +  ^ _ lOp  S  166  '  2  2  6  8  T\ )  6 / „ /, ^ L _ (  +  5  1  1  4  A  0  6  4  _  A  I  0  ^ )  )]+^-(2A A -A 5  4  (A - 1 )(A +  A  A I 2  )  6  1  1  A  1  2  )  Appendix D. The Nonlinear Stiffness Matrix Due to Foreshortening  * 5 2 = ^ ( A - A ) ( A 3 - A ) +^ - [ ( A 8  9a/ ,  2  9  v  2  a  - A  8  6  2  18a/ / — ^ - ( A 2 - A  9  .  33=  5 3 =  =  k  m  ^  l5~^  8  2  +—(A -A )A  8  2  8  2 18a/ .  "  s )  N  - ^ - (  3 -  A  3  3  9  ~^ ^^ 9  -(A -A ) 35 "~' "  9fl/2  3  v  2  +  9  ,  2  A  9 ) (  5  +  8  2  + A  6  ] + ^(A +A  5  5  1  1  )(A +A 6  1  2  1 2  ]  )  1 2  1  2  A  1 2  n  2  )  2  5  )A +(A -A 6  6  1 2  )A  6  5  4  . \ 3a/ /., „ -\ 6/ / ) +—(A ,-2A ) + ^ - ( A - A  1  1 1  )A  1 0  )  2  2  2  4  1  1 0  )  2 2  ^ (A - A )(A + A „ )  +  4  10  5  ~~^)^-6 ~(A2  3  A ) A ] - - ^ - ( A + A )(A + A,,)  _  8  5  6  Mi(A -A )A -^1[(A A )A ( A _ )  8  )A +(A -A  ">  -  3  1 1  1 0  J + -J2— ( A - A  1 2  3  5 + A  4  ]+^[(A +A  1 2  2  2  6  9  6 +  2  9  -2(A3-A )A  1 1  1 0  9  "  (A, - A ) ( A - A ) + ^ [ ( A 9  A  ~^ ^ ^^~[^  35  v /  ) A  2  ^  +—[(A  1 2  4  " ^ "140^  3  8  (A - A J A - A ) - ^  / ? ) / 3  - A  3a/ / , \ 6/ / , +—(2A +A )-^L_(A -A )  x  .v.  2  3  ~3~5~"^  =113  2  72a/ , ~35~^  *« - ^  k  (A -A )  8  6  w ) ( A  ^I2=-^(A2-A8)(A3-A )+^[(A -A )A  ^ 1 2 2 = —  + ( A  2  )  8  6  4  2  x2 2  ) A  9  3  72a// ^82 = - — ( A  - A  12a/ , \ al , , —(A -A )A +—(A +A )  2  ^ 6 2 = — l 2-A8) A  3  167  5+  n  v /  3  "  n  5 +  5  140 " ' ' """''" ' "  n  LV  -A )A 9  1 2  ;  5  v  A | i  5  A i |  8  N  ]_^:(  5  A 4  _  A i o  )  2  "11/--11J  -2(A - A ) A ] - - ^ [ ( A 2  12  6  + A , ) A +(A - A ) A ] 2  5  6  1 2  U  Appendix D. The Nonlinear Stiffness Matrix Due to Foreshortening  6/ / k =-f-(A -A ) + 5p  / — [/,(A -A )(A A )-/ (A -A )(A A )] 5p  J  2  i  2  35?K<» 7 A  +  =~  (2/ V  ^  =^  +/J/  8  6 +  5  n  lOp  y  u  (/ +2/ V  =-^"  + 2 /  7  lOp  9  2  lOp  ^.24 = ^M ^  (A3-A )(A -A 4  ^  8  4  3  )  (A -A XA -A 3  9  (A - A,XA 2  4  2  2  . \  / 5  *.» = - ^ " ( A , - A X A 9  4  4  5  +A „ ) - ^ ( A  4  A  2  5  3  -^XAs  - A )  2  1 0  [(A - A „ X A + A „ ) - ( A - A X A 9  6  1 0  3  9  12  X4A  0  2//,.  »H 3  al^_ 140  +  1 1  6  (A -A  +  + A  5  (A - A , X 4 A  5  30p  "^^^0^  3  3  1 1  X4A -A )  10  4  fy^  9  =  5 +  1 0  4  A  .  (A -A  - ^ ( A , - A ) - ^ - ( A , - A . X A , + A„) ^ ( A  *" " ^ ^  9  (A - A X 4 A - A )  /  L 0  4  ;  o  )-( "t ^ 30p  1 0  - A )+  4  .; 30p U  +  ^ 30p  1 0  . , 2 a/ ,.  6a/ ,.  1 0  2  +  3  5  (A - A )(A - A ) + ( \*  ^ lOp  s  A  2  6  1 2  'A «)]  '>,A )-(/,A A +  +  4  \  6/ / -^-(A _A ) 5p  4  2  44  168  . - A )  -A ) 5  - A )  I2  n 2 2  10  4  1 1  6  a/V., 1 4A +-A -A A v J 2  +  j  2  q  1  5  +A )]+|1(12A A -2A.A,, - 2 A A , d  S  +||i(A , + 2 A 2  + A, )] + | ^ ( A 2  I I  A  1 2  l l  +A A )  6  IL  A - A ) 2  S  M  +A A 6  U  +A A  al - A ) +^(12A  S  - | A  1 2  A  S  5  *« ~ ( A 35  2  "A ) 8  2 +  ^ ( A  2  - A,XA  6  +  A  1 2  ) ^ - ( A +  2  4  1 0  2 6  +A  2 2  -2A A„) 6  12  Appendix D. The Nonlinear Stiffness Matrix Due to Foreshortening  K N 6  al = \40 ^  al  169  5  "  2  A  g  )  (  A  s  +  n ) - (  A  / ^.26 = ^ " ( A - A ) ( A + A 4  2  8  A  - A )(A + A )] +  3  9  / / )-^-(A  6  1 2  A„A  210  +A A  1 2  5  +A A „ - - A A  1 2  6  5  6  5  6  1  2  - A  4  2 2  1 0  ) +  A  1 2  , n A + 2A A 6  3  I 2  A 2 ^  --A  6  J  1  a/ |^( 210' 5  ^ i i i i  —  (A,-A,) - ^ - (  35  —  9  A 5  -A ,) 1  +  - ^(A -A, ) 15p T  4  2  0  +  2  A  +  12A  2 1 +  2A A ) 5  N  3  3  3a/ ( A  9  2  G  a/ ^ ( A  3  4  3  *m* ~  3-A )(  a/ (A - A )(A - A ) - ^ - ( A - A )(A - A 35 140  3a/ ""1211  A  2  - A ) 8  2 +  9  6  1 2  )+^ - ( A 210  2/„/ - A ) +^ - ( A  2  5  5  2  - A )(A 8  6  1 2  +A A  +A A„ + A„A  1 2  al i^(A  6  I 2  )  5  4  - A  1 0  )  2 +  2  + 2A A 6  + 12A ) 2  12  2  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0089233/manifest

Comment

Related Items