"Science, Faculty of"@en . "Mathematics, Department of"@en . "DSpace"@en . "UBCV"@en . "Ullrich, C."@en . "2009-05-26T23:43:47Z"@en . "1998"@en . "Master of Science - MSc"@en . "University of British Columbia"@en . "This thesis is concerned with the determination of contact models for linear elastic\r\nsolids which are both accurate and computable in a real time environment. We\r\nreview previous models based on restitution coefficients and penalty functions, along\r\nwith their efficiency and accuracy. A new model is proposed which is both rapidly\r\ncomputable and based on easily measurable material properties. In addition, we\r\nconsider the role that friction plays in contact modeling and examine the various\r\nfriction based effects that each model can produce. We conclude by outlining a\r\nsimulation system in which this model can be used."@en . "https://circle.library.ubc.ca/rest/handle/2429/8295?expand=metadata"@en . "5252785 bytes"@en . "application/pdf"@en . "Contact Mechanics using Green's Functions Interactive Simulated Environments by C . Ullrich B . M a t h . , University of Waterloo, 1995 A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F Master of Science in T H E F A C U L T Y O F G R A D U A T E S T U D I E S (Institute of Applied Mathematics) we accept this thesis as conforming to the required standard The University of British Columbia August 1998 \u00C2\u00A9 C . Ullrich, 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 for financial gain shall not be allowed without my written permission. Department of , A / l The University of British Columbia Vancouver, Canada DE-6 (2/88) Abstract This thesis is concerned wi th the determination of contact models for linear elastic solids which are both accurate and computable in a real t ime environment. We review previous models based on restitution coefficients and penalty functions, along wi th their efficiency and accuracy. A new model is proposed which is both rapidly computable and based on easily measurable material properties. In addit ion, we consider the role that friction plays in contact modeling and examine the various friction based effects that each model can produce. We conclude by out l ining a simulation system in which this model can be used. ii Contents Abstract ii Contents iii List of Tables viii List of Figures ix Acknowledgements xi 1 Introduction 1 1.1 Contact Models 2 1.2 Simulation 3 1.3 Overview 3 2 Contact Model Review 5 2.1 Goals of Contact Model ing 5 2.1.1 Real Solids 6 2.1.2 Rigid Bodies 9 2.2 Contact Model Classification 10 2.3 Impulse Response Rigid Models 10 iii 2.3.1 Newton 10 2.3.2 Poisson 11 2.3.3 Stronge 11 2.3.4 Discussion 12 2.4 Force Response Rigid Models 12 2.4.1 Penalty Models 12 2.4.2 Bilinear 13 2.4.3 Harmonic Penalty 14 2.4.4 Hunt and Crossley 15 2.4.5 Three Dimensions 17 2.4.6 Discussion 17 2.5 Quasi-Static Cont inuum Models 18 2.5.1 Linear Elastic Solids in Equil ibrium 19 2.5.2 Signorini Problem 20 2.5.3 Flexibility Matrices 22 2.6 Impact Loss in Elastic Solids 24 2.7 Summary 25 3 Contact Response Maps 26 3.1 Goals 26 3.2 Contact M a p Model 27 3.3 Dynamic Linear Elastic Solids 27 3.4 Boundary Elements 28 3.5 Contact M a p Determination 31 3.5.1 Linear Systems Theory 35 3.6 Simulation Using Contact M a p s 36 iv 3.6.1 Convolution Integration 37 3.6.2 Complexity 40 3.7 Discussion 40 3.7.1 Relation to Viscoelasticity 41 3.7.2 Work 42 3.8 Summary 42 4 Numerical Experiments 44 4.1 Introduction 44 4.2 Consistency 44 4.3 Gauging Energy Loss 45 4.4 Real Contact M a p s 49 4.5 Summary 56 5 Friction 57 5.1 Introduction 57 5.2 Historical Friction 57 5.2.1 Rotational Friction 58 5.3 Tribology 59 5.3.1 Physics of Contact 59 5.4 Impulsive Collisions with Friction 60 5.5 Solution Existence and Uniqueness 61 5.6 Contact Sticking 63 5.7 Force Response Rigid Collisions with Friction 65 5.8 Discussion 67 v 6 Contact Map Simulations with Friction 68 6.1 Introduction 68 6.2 Previous Work 69 6.3 Rigid Body Simulations 70 6.4 Mode l Parameters and Observations 72 6.5 Contact M a p 77 6.6 Conclusion 84 7 Conclusions 85 Bibliography 87 Appendix A Boundary Element Review 94 A . l Boundary Elements and Finite Elements 94 A .2 Boundary Integral Equations 95 A.2.1 Fundamental Solutions 95 A.2.2 Residual 95 A . 3 Numerical Considerations 96 A . 3.1 Boundary Elements 97 A . 4 Discussion 98 Appendix B Simulation Framework 99 B . l Introduction 99 B.2 Design Goals 100 B.3 Design Overview 101 B . 3.1 Simulator 101 B.3.2 Experiment Design 102 vi B.3.3 Universe Design 103 B.3.4 Model Edit ing 103 B.3.5 Model Design 104 B.3.6 Numerical Editor 105 B.3.7 Summary 105 B.4 Object Components 105 B.4.1 Display 105 B.4.2 State Registry 106 B.4.3 Relationship Managers 107 B.4.4 Engine 1 \u00C2\u00B0 8 B.4.5 T h i n g Database 109 B.4.6 T h i n g 109 B.4.7 Clocks H O B.4.8 D a t a Streams H O B.5 Object Interaction H O B.6 IS Simulator Implementation H O B.7 IS Simulator Improvements/Future Work 114 vii List of Tables 4.1 M a x i m u m force errors between contact map and harmonic penalty model 45 4.2 B a r material properties 51 6.1 Exper imental ly determined collision parameters 72 v i i i List of Figures 2.1 Schematic of Hertz contact scenario 8 2.2 Bilinear force model, showing a sample collision path 14 2.3 Fraction of energy lost during single impact collision. In this case l = n = l 16 2.4 Schematic of Signorini problem configuration 20 3.1 Boundary configuration for contact map determination 32 3.2 Traction response curves (step and impulsive) of a steel bar to a step displacement of 1 0 - 6 m normal to one end. T h e impulse response was computed using eq. (3.16) but not scaled by 1/T 34 4.1 T h e step response curves used to gauge energy loss. T h e steady state behaviour is constant at 1500 N / m in all cases 46 4.2 Derived impulse response curves used to gauge energy loss 47 4.3 T h e computed restitution coefficient for each step response, as a func-tion of the transient peak size 48 4.4 T h e bar used to determine contact maps 50 4.5 Derived impulse maps resulting from 18 deg. impact 51 4.6 Derived impulse maps resulting from 30 deg. impact 52 ix 4.7 Derived impulse maps resulting from 42 deg. impact 52 4.8 Derived impulse maps resulting from 54 deg. impact 53 4.9 Derived impulse maps resulting from 66 deg. impact 53 4.10 Derived impulse maps resulting from 78 deg. impact 54 4.11 Derived impulse maps resulting from 90 deg. impact 54 4.12 Derived impulse maps resulting from step displacement at each angle. 55 6.1 Bar experiment configuration 69 6.2 Post impact total energy 73 6.3 Impact duration 74 6.4 Post impact normal tip velocity 75 6.5 Post impact normal center of mass velocity 76 6.6 Post impact energy fraction due to friction only 78 6.7 Total kinetic energy during collision of bar for various impact angles. 80 6.8 Total kinetic energy of bar for various impact angles 81 6.9 Impact duration for various impact angles 82 6.10 Percentage of energy lost to internal vibrations 83 B . l IS simulator design overview I l l B.2 Experiment design window 112 B.3 A running simulation shown from the user view 113 B.4 T h e scene navigation bar 114 B.5 A n example of a live t h i n g property editor. While the simulation is running, the selected t h i n g s properties may be edited 115 x Acknowledgements I would like to thank my advisor D r . D . K . Pai for introducing me to this prob-lem and for numerous enlightening discussions. D r . Pai's passion for interactive simulation served as both motivation and inspiration. I would also like to thank my parents for providing support at those critical moments, to allow me to pursue this research. Finally, I would like to thank my thesis committee, D r . Ascher and D r . Peirce for taking the time to read my thesis. C . ULLRICH The University of British Columbia August 1998 xi Chapter 1 Introduction This thesis could be classified as research into dynamics for interactive virtual reality ( V R ) . V ir tua l reality is a large and growing field with many diverging applications and research directions. Although many aspects of V R are criticized for their lack of utility, two p a r t i c u l a r l y relevant and important subfields are interactive dynamic simulation and computer aided rapid prototyping. T h e second application has sig-nificant impact on the world in which we live. Unti l recently, interactive V R research was dominated by computer graphics issues related to speed and visual quality. Now that most modern computers by default have little difficulty displaying virtual environments, we are left to consider the properties and dynamics of virtual spaces in more depth than ever before. In this thesis, we focus on algorithms, issues and problems related to accurate and rapid dynamics. Specifically, we are concerned with the fidelity of simulated dynamics of in-teracting virtual bodies and in particular their collision properties. M a n y current V R systems provide collision behaviour modeling but either make no attempt to 1 achieve physical accuracy ([38], [17]), or provide full blown continuum modeling at the cost of non-interactivity ([12], [33]). In this thesis, we seek a middle ground solution in which physics are reason-ably accurate, but at the same time all calculations can be performed at or near real-time speeds. 1.1 Contact Models Models of contact between solids date back to at least Newton. Unti l recently, most contact models were of the impulsive type in which the collision is assumed to have zero duration and result in discontinuous changes in velocity. Such models are useful in situations where angular momenta can be largely ignored, but encounter a host of problems in complex collisions involving friction. Recently, collision models have been proposed which are based on visco-elastic continuum theory. These models are collectively known as penalty models because the construct collision forces as functions of the penetration depth of two colliding objects. T h e forces tend to resist or penalize deep penetrations, hence the name. Though these models are closer to the actual physics of contact they are heuristic and have parameter determinacy problems. A third type of collision model simply uses the full set of partial differential equations (PDE's ) governing elastic solids. This formulation results in highly ac-curate solutions describing the behaviour of colliding bodies. This method suffers from severe complexity issues however and real time solution of these equations is possible only on extremely powerful computers. T h e model introduced in this thesis attempts to leverage the accuracy of the P D E approach by precomputing as much information as possible. This pre-2 computation results in a set of functions which describe the characteristic behaviour of a solid under penalty like conditions. Given these characteristics (called Green's functions) it is possible to construct accurate solutions in real time using convolution integrals. We demonstrate the efficacy of this approach in sample ideal collisions and in an actual real-world configuration. T h e latter demonstration involves frictional effects which we discuss in some detail. T h e introduction of friction into contact mechanics presents a host of unique problems and forms an area of active research in this field. 1.2 Simulation Finally this thesis considers the object oriented design of an interactive simulation environment. It is the intention of this author to produce a flexible environment in which dynamical models of contact and a variety of other object properties can be tested with minimal effort. A simple prototype which has already been implemented gives a reference for this discussion. 1.3 Overview T h e thesis begins with a review of the usual contact algorithms including impulsive, linear force and quasi-static continuum models. We deliberately do not consider issues related to tangential forces such as friction, which we reserve for a later chap-ter. Following this, we introduce the contact map algorithm for generating Green's functions and using them in a simulation environment. Chapter 4 presents some artificial and real simulation results from the contact map algorithm. In chapter 5, 3 we review the current knowledge of tangential surface forces, pr imari ly friction. We then present some numerical results comparing the contact model wi th friction to actual experimental results. In the appendices,, we review the mathematics of the boundary element method and discuss the proposed simulation environment. 4 C h a p t e r 2 Contact Model Review 2.1 Goals of Contact Modeling When two solids come into contact, many different forces act on each body. T h e forces arise from local effects near the contact region such as elastic response and friction but also include global effects which are functions of geometry, mass dis-tribution and others. Currently, a full theoretical understanding of all these effects is incomplete and experimental data are very limited. In addition, simulation of contact using more than one or two of these properties is usually computationally expensive. A primary goal of contact modeling is to produce an acceptable tradeoff between the cost of computing contact forces and the accuracy of the results. In this section, we discuss the theory of contact and some of the more popular models that have been proposed to resolve this process. 5 2.1.1 Real Solids Realistic modeling of solids which undergo deformation and respond to collision related forces is typically done with the elastic continuum model and various exten-sions for thermal, plastic and cracking effects. We postpone a detailed analysis of the elastic continuum model of solids for the moment, and consider here only the processes which occur near the surface of contacting bodies. Hertz Model T h e first model we consider was proposed by Hertz over a hundred years ago. We present this model in a separate section because it introduces the basic approach that many contact models take. Full details of this model are available in many books, for example [36]. T h e Hertz theory relates to two elastic bodies coming into contact at a point. Hertz made four basic assumptions when formulating his theory; 1) small strains (surface displacements within the elastic limit) 2) the two contact surfaces are continuous and do not conform 3) the contact region is much smaller than the local radius of curvature 4) there is no friction. Consider two bodies which come into contact at a point. Define a co-ordinate system with origin at the contact point, and z-axis normal to the contact. Each body has a local radius of curvature (e.g., Gaussian curvature) Rix, Riy and R,2X, R-2y for bodies 1 and 2 in the x and y directions respectively. Further, we assume that the principle axes of both bodies are aligned and coincident with the z co-ordinate axis. Hence for x and y close to the contact point, we can express the undeformed 6 surfaces as for k = 1,2. Define the relative radii of curvature by Rx \u00E2\u0080\u0094 and Ry = R R D ly_s_nv \u00E2\u0080\u00A2 Then we can write the separation between the two bodies as rl\y-\-I\2y 2 2 h(x,y) = z1+z2 = ^- + ^--. (2.2) Note that level sets of h are ellipses. We assume further that there is some local deformation around the contact point of maximum magnitude along the z axis. Thus , within an elliptical contact region, we have wi(x,y) + w2(x,y) = \u00E2\u0082\u00ACi + e 2 - h(x,y), (2.3) where Wk(x,y) represents the vertical deformation of a point on body k at position (x,y). B y assuming the contact region is elliptical the displacement can be expressed as a function of the contact pressure, where v is the Poisson ratio and E is Young's modulus of elasticity. Hertz showed that an ellipsoidal pressure distribution could satisfy (2.4) in this case. Combining the above equations, we arrive at the Hertz model of contact, e* f r p&nWdr, y y2 K J Js J(x - a2 + (v - n)2 2i2i 2 V { ' > where e = \u00E2\u0082\u00ACi + e2 and E* = -g^-1- -\\u00E2\u0080\u0094^r2-. This equation relates the pressure to the displacement in the contact zone. R l Figure 2.1: Schematic of Hertz contact scenario. 8 In the case where the contacting bodies are spherical, as shown in 2.1, it is relatively straightforward to show that the total pressure P is related to the displacement e by, P=^E*Rel. (2.6) When the bodies are non-spherical, it becomes difficult or impossible to construct analytic solutions, although perturbations exist for mildly elliptic situations. Most 3 authors simply refer to the relationship / = ket with force / and elastic constant k as the Hertz model. T h e Hertz model is a cornerstone of modern contact theory yet it does not account for many commonly observed collision phenomenon such as energy loss, friction and permanent deformation. We discuss these effects and models for them below. 2.1.2 Rigid Bodies In light of the large number of collision phenomena, most researchers have focused on a few of the key principles. This thesis will be no exception, we will consider physical contact between objects which are accurately modeled as rigid bodies. B y rigid bodies, we mean that collision forces are small enough that significant global deformation does not occur. T h e rigid body assumption restricts our consideration to collisions which are of relatively low velocity and between hard objects. In this narrow band, there are a large number of models and theories which attempt to bridge the gap between the scientific discipline of tribology and practical engineering. For an enlightening discussion of this connection, see [53]. 9 2.2 Contact Model Classification Contact models can be classified into two groups, those which attempt to model contact forces explicitly and those which only attempt to model the effects of a collision. T h e former are called force response rigid models, while the latter are impulse response rigid. This taxonomy was originally introduced by Chatterjee 2.3 Impulse Response Rigid Models Impulse response rigid models only attempt to model the physics immediately before and immediately after a collision. T h e forces generated during contact are not modeled, only their net effect on object velocity. We express this with, where the applied impulse is Ap and the collision starts at time t = 0 and runs until the bodies separate at tj. One of the fundamental assumptions of impulse response rigid models is that tj is very small relative to the rest of the simulated processes. T h e category of impulse response rigid models is further subdivided by three different definitions of a dimensionless coefficient which relates pre-contact states to post-contact states. This quantity is generally known as the coefficient of restitution. We briefly describe the three definitions here, however we must postpone an analysis of their relative merits until after we discuss friction in Chapter 5. 2.3.1 Newton T h e first person to consider this problem in detail was Newton, who proposed that a dimensionless constant e \u00E2\u0082\u00AC [0,1] relate the velocities of bodies before and after a [14]. (2.7) 10 collision (2.8) where is the post collision normal velocity and is the pre-collision normal velocity. A choice of e = 1 corresponds to a perfectly elastic collision whereas setting e = 0 indicates that all kinetic energy is dissipated during a collision. 2.3.2 Poisson Poisson revised Newton's model by dividing an impulsive collision event into two phases, the compression phase and the restitution phase. T h e transition from com-pression to restitution is defined as the point where the normal velocity is zero. T h e resulting coefficient is expressed as a ratio of normal momenta, where pl^ is the integrated normal force of restitution from the time of maximum compression to the final contact time. Also, p~^ is the integrated normal force of compression, from the time of initial contact to the time of maximum compression. Keller presents a detailed analysis of this model for 3D impacts in [39]. 2.3.3 Stronge T h e third widely used definition was proposed in 1991 by Stronge [57] and relates the work done during compression to the work done during restitution. Precisely, (2.10) where Wr, Wc are the work done by the normal force during compression and resti-tution respectively. This model was proposed after it was discovered that the above two models generated energy in certain types of frictional collisions. (2.9) 11 2.3.4 Discussion It is important to note that the above three laws give exactly the same result in contact configurations which do not involve friction during the collision. T h e latter two were constructed to resolve apparent violations of conservation of energy noted by [39], [8] and others. It is also well known that experimentally, the coefficient of restitution is a function of impact velocity (varying approximately as the fourth root, see [25] ), impact geometry [55] as well as a host of other surface and material properties. These functional dependencies lead one to question the merit of just a single restitution coefficient for two colliding bodies. We discuss the impulse response rigid models in more detail in Chapter 5 on frictional contact where the differences are much more clear. 2.4 Force Response Rigid Models Models which attempt to accurately generate the forces that occur during a collision are called force response rigid. These models attempt to be more physically accurate than impulse models but typically come at the price of increased computational cost. Also, when these models are coupled with a friction model, we are ensured that no anomalous energy will be generated by the collision, unlike impulse models. 2.4.1 Penalty Models T h e simplest force response rigid models are penalty models, in which the force is a function of the current penetration depth of two bodies. This corresponds to the term e in the discussion of the Hertz contact law. A simple analysis shows that 12 models which depend only on the local value of this parameter cannot dissipate energy through collision processes. T h e addition of a viscous term depending on e makes energy dissipation possible. Penalty models can thus be expressed in the following functional form, ' /(\u00E2\u0082\u00AC,\u00E2\u0082\u00AC) if \u00E2\u0082\u00AC > 0 (2.11) 0 otherwise Note that this model depends only on the local penetration depth and its derivative (i.e., / only depends on the value at time t). This is an important distinction, because in the next chapter we consider non-local laws of this type. T h e Hertz model presented above can be considered as a nonlinear example of this type of law, with no dependence on k. We consider here further examples, the bilinear model of Stronge [58], the harmonic penalty model and the nonlinear model of Hunt and Crossley [32]. We analyze each of these models in the context of a one dimensional collision. In this case Newton's second law gives, me + /(e,e) = 0. (2.12) 2.4.2 Bilinear T h e bilinear model is a penalty realization of the energetic coefficient of restitution 2.10 defined above. In this model, there are two spring constants, one for compres-sion k\ and one for restitution k2. Given a coefficient of restitution e, we define k2 = fi\" which for e \u00C2\u00A3 [0,1] gives a higher spring constant for unloading, see figure 2.2. Contact is broken when the force returns to zero. T h e energy lost to collision is equal to the area between the two lines which is consistent with 2.10. A s noted by [14], this method suffers from numerical difficulties when e is close to 0 because the unloading slope can be arbitrarily steep. Chatterjee suggests 13 Force e Figure 2.2: Bilinear force model, showing a sample collision path. a possible alternative, where the loading and unloading lines are constrained to pass through the origin, and a 'jump' occurs at the point of maximum compression. In this case, &2 < k\. 2.4.3 Harmonic Penalty It is also possible to achieve energy loss by including a so called dissipation term in the penalty law. This is easily seen in the simple 2nd order linear O D E , me + ce + ke = 0, (2.13) which is a special case of (2.11) with /(e, e) = ^(\u00E2\u0080\u0094ce \u00E2\u0080\u0094 ke). T h e velocity dependent term cx controls the amount of damping in the solution. Adding a damping term to the penalty law removes the need to adjust the spring constant k during the course of the collision. This approach is used by many researchers as it is simple and easy to implement. A complete simulation system based on this type of model, but generalized to 3D, is presented in [26][27] and we examine it in chapter 5. Chatterjee analyzes this model and constructs an approximation to the co-14 efficient of restitution. Define (3 = \u00E2\u0080\u0094| then a good approximation ( maximum relative error less than .06 ) to e is given by, 1 (2.14) e = 1 + 2.5/3 + 4/32 ' One major drawback of this model is that due to the non-zero initial relative velocity, a force discontinuity occurs at the beginning and end of contact. T h e size of this discontinuity is proportional to the impact velocity. This may be a concern for simulation environments in which this type of force model is used to generate haptic feedback, since the discontinuity will certainly be felt by a user. 2.4.4 Hunt and Crossley T h e issue of force discontinuity is one of the motivations behind the seminal work of Hunt and Crossley [32]. They consider a modification to the the linear equation Using this modification, the initial force is zero for zero displacement, independent of the initial velocity. T h e nonlinearities in this equation mean an analytic solution is in general impossible. T h e restitution coefficient, as defined by Newton can be approximated for the case / = n = 1, see [14]. Hunt and Crossley derive a more general expression in the case when A is small (see [32] eqn. 14 ). Relating the energy loss through a collision to the various collision parameters can only be done numerically. We calculate this by integrating eqn. (2.15) for a variety of A and k values in the case / = n = 1, see figure 2.3. T h e initial conditions were e(0) = 0, e = \u00E2\u0080\u00941. Integration was terminated,when e(t) = 0 for t > 0. T h e delta energy was calculated from the ratio of kinetic energies before and after the above TO6+ (Xel)e + ken = 0. (2.15) 15 Hunt Crossley Energy Loss Figure 2.3: Fraction of energy lost during single impact collision. In this case l = n = 1. integration. In the case where m \u00E2\u0080\u0094 n, Hunt and Crossley construct an explicit method for determining A based on a velocity dependent restitution definition. e = 1 - avi, (2.16) where V{ is the impact velocity and a is a small coupling constant. In this case we can write 3 A = -ak. (2.17) 16 For further details, we refer the interested reader to [32]. 2.4.5 Three Dimensions It is important to understand that the conclusions of the one dimensional collision analysis above do not necessarily extend naturally to three dimensions. T h e reason for this is that in three dimensions, we must deal with mass matrices, spring and damping tensors. T h e interaction between these three quantities becomes more subtle and even more dependent on the particular impact configuration. T o reduce this problem somewhat, Goyal , et. al. [26] considers only the spring and damping constants in the normal and tangential directions at the collision point only. In higher dimensions, the number of external model parameters in the force response rigid models presented thus far becomes almost prohibitive. Simplifying assumptions such as those in [26], [27] are required to have a workable simulation. However since the parameters of all the force models have somewhat tenuous relation to the physics of contact, one cannot expect accurate solutions to result from these models. 2.4.6 Discussion In all of the models discussed so far, with the exception of the Hertz model, no mechanism is given to determine the collision parameters. Chatterjee tries to relate all penalty model parameters to the coefficient of restitution. Though this may seem like a reasonable goal, in practice the coefficient of restitution is a difficult quantity to measure, requiring specialized equipment and high speed measuring devices for even simple scenarios (see [55], [20]). Including the important dependencies of e on geometry, velocity and material 17 properties, we quickly end up with a large number of scenarios in which e must be determined. For general simulation environments, a huge database would be needed, which would grow at least as the square of the number of objects. It is highly doubtful that e could even be measured for all possible combinations, even with an advanced facility such as the A C M E testbed [47]. Given the requirement that a simulator run at interactive speeds and pro-duce realistic and accurate data, we can evaluate the models presented thus far. Although impulsive models work very quickly in practice, the restitution coefficient has many subtle dependencies and many measurement related difficulties. T h e force rigid models require that collisions be explicitly integrated, however there is some difficulty in determining the coefficients of the force relations especially in more than one dimension. Accuracy will also have to be traded off to achieve reasonable integration speed on all but the fastest computers. Given these options, we look beyond these simple models to the much more sophisticated models that are currently used in mechanical and civil engineering for performing precise contact analysis. 2.5 Quasi-Static Continuum Models T h e framework for the study of mechanical contact when only small deformations are present is the continuum theory of linear elastic solids. This model consists of a system of P D E ' s together with appropriate boundary and initial conditions and can accurately model many contact and low energy collision scenarios. A s we shall see however, the numerical solution of these problems is both complex and time consuming, effectively removing them from consideration for interactive simulation as well. 18 2.5.1 Linear Elastic Solids in Equilibrium We begin by reviewing the relevant definitions in the subject of continuum mechanics as they relate to elastic solids. Consider a point set Q \u00E2\u0082\u00AC 3t 3 with elements x = (x\, X2, X3). For most applications, we will assume that Q is compact. Now define transformed coordinates of each point by y; = X ; ( x ) a n d define the displacement vector by u = y \u00E2\u0080\u0094 x. Assume further that the transformation x l s invertible and at least C1. F r o m these fundamentals, we define two important quantities, the strain tensor e,-j = | ( u ; j + Ujj) and the stress tensor CT 8J(X) = CT;J(X, et-j). T h e function cTjj encapsulates the material properties and is known as the constitutive relation. When Q is a linear elastic solid, the constitutive relation is where Eujk is the Hooke tensor, completely analogous to the spring constant k in the previous section, but with 81 components. Note that adopt the summation convention in which repeated indices indicate summation. F r o m a force balance we see that, where 6 t(x) are body forces which act over the entire region Q (e.g., gravity). This equation is a force balance and is valid in the absence of time dependent forces which act on time scales smaller than the time for a sound wave to traverse the solid (equilibrium). Full details of the theory of elastic solids can be found in many books, for example [62]. In addition, we consider a dynamic version of this theory in the next chapter. Uij = Eijkieki (2.18) \u00C2\u00B0\"tj,j(x) + 6;(x) = 0, (2.19) 19 Figure 2.4: Schematic of Signorini problem configuration. 2.5.2 Signorini Problem T h e study of elastic solids in contact situations was pioneered by Signorini in 1933 [52] when he formulated a general constrained elastostatic problem. Signorini con-sidered the problem of a linear elastic solid coming into contact with a rigid, semi-infinite foundation, see figure 2.4. Consider an elastic solid 0, which is free to move in the vertical direction in 3-space. Below the solid is a semi-infinite foundation. We let Tc be the region on the surface of Q which will contact the foundation if Q is translated vertically. Assume that a chart (xi, x2)+u3(x1,x2, (xi,x2)) < rp{x1+u1(x1, x2, {xi,x2))). (2.20) Equation 2.20 simply states that the two surfaces cannot penetrate each other. For 20 simplicity we rewrite this using the invertibility of the transformation x, {yx,y2) + u3(y1,y2) < Myuy*)- (2-21) B y considering small increments to the displacement u and linearizing, 2.21 is rewrit-ten as \u00C2\u00AB \u00E2\u0080\u009E ( x ) - g{y) < 0. (2.22) Here \ A + ( i \u00C2\u00A3 ) 2 + ( i \u00C2\u00A3 ) 2 is the normalized separation and w n (x) = u(x) \u00E2\u0080\u00A2 n(x) is the normal displacement. Let r n ( x ) and T T , ( x ) represent the normal and tangential surface stress in the contact region Tc respectively. Then for frictionless contact we require, r T i ( x ) = 0 (2.24) r \u00E2\u0080\u009E ( x ) < 0 . (2.25) Finally, because the contact stresses only occur when there is actual contact, we have the linear complementarity condition, r n ( x ) ( u n ( x ) - 0 ( x ) ) = O. (2.26) T h e Signorini problem is the set of relations 2.19, 2.21, 2.24, 2.26 and any boundary conditions on dQ \u00E2\u0080\u0094 Tc-There are a number of techniques for solving such problems in general. In [40] the problem is formulated as a variational principle which leads to weak solution techniques such as F E M (Finite Element Method) . T h e authors of [2] also use an integral form but apply instead a boundary integral (or B E M , Boundary Element Method) technique. T h e B E M approach has the advantage that only the boundary 21 r needs to be discretized. We discuss this approach in detail in the next section and the following chapter. T h e solution to the Signorini problem and its extensions (two body contact, dynamics) provide high accuracy and depend only on easily known material param-eters such as the shear modulus G and the Poisson ratio v. However the numerics of this approach are formidable and computationally intensive. T h e author has never seen reference to any system which claims to compute such models at or near real time. 2.5.3 Flexibility Matrices Finally we consider an approach to solving the contact problem which is similar to the method presented in the next chapter. T h e boundary element method (see Appendix A ) for the two contacting bodies above yields a dense system of linear equations, Ay = b, (2.27) where A is a discretization of the boundary integrals, y contains the unknown displacement or traction as appropriate and b is a vector of known values of the displacement or traction at each boundary node. T h e flexibility matrix approach [61] to solving 2.27 rapidly, is to pre-compute the elastic response of each node and then solve the boundary integral only in the contact region. T w o matrices are needed to accomplish this, the displacement matrix represents the deformation produced by an incremental external load in the contact region. T h e flexibility matrix gives the displacement when a unit normal or tangential force is applied on the boundary. To calculate the displacement matrix for each body, external forces are ap-22 plied to the appropriate boundary elements, fixed displacement boundary conditions are applied to those elements which are initially in contact and free boundary con-ditions are applied elsewhere. Then using a B E M solution technique, the displace-ments of the elements close to the contact region are calculated. T h e displacement matrix D is iV X 3 for three dimensional contact problems ([61] consider only the 2D case) and is used as follows, S = DP, (2.28) where P is the applied surface traction, and S is the resulting surface displacement. Generation of the flexibility matrix is somewhat more complex. Thi s matrix relates a normal or tangential force at the contact surface to the resulting displacement (of all elements). T h e contact surface is the average of the surface elements weighted by the Young's moduli of the two contacting bodies, 9 a = 9 l . i g i . 0 2 ) - ^ \u00E2\u0080\u0094 . (2.29) This model assumes a reasonable prior knowledge of this contact angle to compute the flexibility matrices for each object and that only small changes will occur during the collision. Using this information, the boundary element method is successively applied with a normal and tangential force for each element and free boundary conditions for all others. T h e result is a rank 4 tensor, F = [fjj], (2.30) where i is the node displaced, j is the node at which a force is applied, k is the direction of the displacement component and / is the direction of the applied traction. Multiplication by a traction matrix T gives a displacement matrix S, FT = S. (2.31) 23 Using an iterative, incremental algorithm the authors [61] show how to use eqns. 2.28, 2.31 to rapidly solve contact problems. We refer to their paper for further details. This approach provides a fast and extremely accurate method for solving quasi-static contact problems. However this algorithm requires a fairly significant amount of knowledge about the type of contact which is occurring to design the deformation and flexibility matrices. T h e authors do not comment on the ease with which the various choices are made. In addition the coupling of the two bodies through the contact angle means that two different flexibility matrices will have to be computed for each object pair. For environments with a large number of objects and/or with complex objects, the storage requirements and the pre-computation times could become prohibitive. 2.6 Impact Loss in Elastic Solids The Signorini and flexibility approach discussed above both assume that the contact forces are applied in a quasi static manner. This assumption is justified in situations where loading and unloading in the contact region take place over a long time period (relative to the wave speed in each body). For collision situations this assumption is often violated. Collision events typically modeled with impulsive laws take place on the scale of tens to hundreds of microseconds. In [55], the authors demonstrate collision events in which as much as 50% of the initial kinetic energy is lost to elastic waves during a collision event. Similar results are reported in [59]. Examining the contact models presented thus far, it is clear that none explic-itly consider this effect. One reason for this omission is that energy loss to elastic 24 waves is very dependent on the contact geometry and the geometry of the solids in contact. Stronge reports that energy loss for a collision between two cylinders is 44% but drops to only 10% for a collision between two spheres. A n y of the models presented which include damping effects could include elastic wave loss implicitly. This approach would necessitate a field of coefficients covering the surface of each body, which could also be a function of the contactors. T h e next chapter presents an alternative to the models presented thus far which also includes effects related to geometry and elastic waves. 2.7 Summary In this chapter, we considered contact models that can be used either for rapid, inaccurate solution of impact configurations or for slow accurate solution of contact configurations. T h e impulsive and penalty models are not suited for long term contact situations because their assumptions lead to either discontinuous dynamics or very stiff O D E ' s . B y comparison, continuum models discussed above are of little use in impact situations in which the internal dynamics of a solid are relevant, due to their quasi-static nature, but provide very accurate solutions for bodies which come to rest against each other. T h e inaccuracies of the impulsive and penalty models is a result of the severe assumptions on the colliding bodies. T h e critical dependence of contact mechanics on geometry is only explicitly included in the Hertz model. In the next chapter we attempt to construct a model that addresses this issue directly. 25 Chapter 3 Contact Response Maps 3.1 Goals This chapter develops a new approach to contact modeling which tries to incorpo-rate geometric and transient effects, without greatly sacrificing performance. T h e approach is to determine the response of a linear elastic solid to a small deformation at each point on the surface and save all these response curves as a contact map. This information can then be used during a simulation to more accurately model the response. One of the primary goals of this approach is to remove arbitrary and/or difficult to measure parameters from the collision model. Nearly all the models discussed in the previous chapter suffer from this parameter indeterminacy problem. In addition, none of the non-continuum models make a direct attempt to handle the dependency of the dynamic collision forces on global object geometry. T h e contact map model attempts to address these deficiencies while remaining rapidly computable. 26 We begin by discussing the materials to which this model is suited, then discuss the dynamic equations which govern them. We take a short detour into boundary elements to describe one possible way to compute a contact map. Fol-lowing this, we discuss a time stepping algorithm for handling collisions that uses contact maps. We conclude with some remarks on the algorithm consistency and energetics. 3.2 Contact Map Model We restrict our attention in this chapter to linear elastic solids which do not deform significantly during a collision process. This category includes most hard objects like steel, aluminum and cast iron but excludes deformable materials like rubber and clay. T h e classification of an object into the hard category is also a function of velocity since, for example most plastics do not deform significantly at low speed. Materials of this type are modeled faithfully [62] by the dynamic Navier equation discussed in the next section. 3.3 Dynamic Linear Elastic Solids In situations where dynamic effects are important, 2.19 is modified to include accel-erations, oihj(t, x) + bi(t, x) = pv.i(t, x ) , (3.1) where U{ = d \" J ^ ' x ^ and p is the mass density. This modification changes the equation into a hyperbolic P D E and greatly complicates the solution process. When the solid 0 is linear, elastic, homogeneous and isotropic the 81 com-ponents of Eijki reduce to just 2 constant values. In general these two parameters 27 are given by the Young's modulus E and the Poisson ratio v. In 3.1, this gives rise to two distinct wave speeds which we denote c\ and c2. In terms of the material constants, they are C l ~ p{l + v){l-2v) [ 6 ' 2 ) and ^ = wnry (3-3) Waves of speed Ci are known as pressure or p waves and those with speed c2 are shear or s waves. Using these constants, we rewrite 3.1 in terms of displacements to obtain the Navier equation for elastic solids, ( C l _ Cf)Uj,ji + clui,jj + = (3-4) We note that other types of waves exist in elastic solids, Raleigh waves (sur-face waves) are one notable type. In our model, these other wave types are either implicit or are not modeled. 3.4 Boundary Elements There has been a great deal of work over the past two decades in the boundary element community on the solution of elastostatic and elastodynamic problems [1]. Boundary element ( B E M ) solution techniques are ideally suited to contact problems because all the unknown quantities of interest are on the boundaries of the colliding bodies. In contrast, F E M (finite element) approaches deal with the entire domain of both bodies and solve for stresses and displacements at every point. T h e basic idea of B E M is to reformulate the partial differential equation as an integral equation. There are several ways of doing this for the Navier field equation. We briefly outline one approach and refer the interested reader to [22] or [1]. 28 Denote the solution to E q . 3.4 by u, with stress field a and body forces b. Consider a second solution u* to the field equation 3.4 with different initial conditions for the same geometry, with stress field a* and body forces b*. T h e weighted residual of the two solutions can be written using the Graffi reciprocal relation [22], / iakj,j * u*k)dQ + / p(bk * u*k)dQ. - / p(uk * u*k)dQ = 0, (3.5) Jn Jn Jn where * denotes a convolution product in time. Using the divergence theorem we can rewrite this as / (pk * u*k)dY + / p(bk *u*k + uQku*k + uoku*k)dQ Jr Jn = (P*k* uk)dT + / p(b*k *uk + u*Qkiik + Uokuk)dQ, Jr Jn (3.6) where uok = Uk(x,0), uok = Wfc(a;,0) are initial displacements and velocities respectively, and pk is the fc'th component of the surface traction. Now, we generate the reciprocal solution u* as a solution of the field equations when an impulsive body load is applied at the point x1 in the direction / , i.e., pbl = 8{t)8{x-xi)8lk. (3.7) Further, assume the initial conditions are identically zero and that the body forces of the primary solution may be neglected. Then we can write E q . 3.6 as ut{x\t) = J^{u*lk*Pk) ~ {uk*p*lk)dT, (3.8) which gives the displacement at point xl in direction / as a function of the boundary. Note that this same expression can be obtained from a variational principle (see [2]). 29 Because the integrals in 3.8 are difficult or impossible to solve analytically, numerical approximation is almost always introduced at this point. Observe that we need only to discretize the boundary, reducing the dimension of the problem by one. In contrast, the F E M approach requires the entire domain Q to be discretized. Typica l B E M algorithms divide up the boundary into piecewise' constant, linear or quadratic elements. Corners require special treatment. Also, it is not uncommon for the resulting integrals to be highly singular, in which case, the Cauchy principal value is used ( see Appendix A ) . A s a result of the discretization of the boundary and approximate evalua-tion of the integrals, we obtain a system of linear algebraic equations for the un-known boundary displacements or tractions. Although the system is much smaller than a similar F E M system, the matrix presents no special symmetry or sparseness properties and is usually solved using dense matrix techniques such as Gaussian elimination. When the integral 3.8 is discretized, the resulting linear equations can be rearranged so that for nodes on which displacements are given, the tractions are solved for and for nodes with given traction, the displacement is solved for. T h e result of this procedure is a linear system of the form, Ax = b, (3.9) where A is a matrix of boundary integrals over each boundary element, b is the known traction or displacement at each boundary element and x is the unknown traction or displacement at each boundary element. This gives kn linear equations for a body with n boundary elements and is a small constant that depends on the order of the integral approximations and the dimensionality of the problem. Mullen and Rencis [45] review a number of iterative schemes for solving 30 boundary element equations in the form of eqn. (3.9). Due to the dense structure of the equations, relaxation schemes such as Gauss-Sidel and Jacobi do not compare favorably with direct Gauss elimination. In fact, the only iterative methods which compared favorably with Gauss elimination were the conjugant gradient (Lanczos) type schemes [4]. We end this section by noting that although boundary element techniques are particularly well suited to contact problems, we are by no means restricted to their use. Highly accurate simulation of elastic bodies can be done with finite elements or, given the proper equipment actual physical experimentation can be done. 3.5 Contact Map Determination We would like to characterize the forces... which occur on the surface of a linear elastic solid when it is deformed in various ways. Surface forces are given by the stress tensor evaluated on the surface, p = c r - n , (3.10) where n is the local surface normal. T h e traction p can be thought of as a surface pressure. In the following, we assume that bodies in contact have the same sized boundary elements or are integral multiples of each other. In addition, we assume that the boundary of the body Q has been discretized into a finite set of nodes {x i , . . . ,x n } . To determine the forces which will occur at a particular boundary node X j , (see fig. 3.1) we consider the traction generated by a normal step displacement at that point (we denote this boundary segment by T i ) on the surface, dj(0\"= -DH{t)hj, (3.11) 31 32 where D is a constant scale factor, H (t) is the Heaviside function and rij is the unit normal at point X j . A l l other points on the surface (denoted by T2) are free, with Neumann condition, where pk is the surface traction at point x^. Because of the dual nature of boundary conditions in B E M it is an ideal candidate for determining the response with mixed boundary conditions such as these. Once the step response has been determined at various points on the bound-ary T, we compute the impulse response by numerically differentiating the response function. We discuss how this is done in some detail in the next section. T h e re-sult is essentially the numerical determination of the Green's function for the IB V P consisting of eq. (3.4), (3.11), (3.12) and zero initial condition. We denote the computed Green's function for body A by where x \u00C2\u00A3 {x i ,X2 , . . . , x \u00E2\u0080\u009E } , is the set of surface nodes in the discretized boundary. For simplicity, we drop the explicit x reference in what follows by assuming that the contact point does not change significantly during a collision event. For collision type impacts with a time scale in the tens to hundreds of micro-seconds this assumption is reasonable. For sustained contact events, some motion may occur and tangential elastic forces will need to be taken into account. We note that given the well posed nature of the above I B V P , we are assured of the existence of the Green's function. Also, since the problem is linear uniqueness follows directly. Pk(t) = 0 V x f c ^ X j (3.12) 9A{x,t), (3.13) 33 2.5e+08 -3e+08 1 1 1 1 1 1 1 1 1 1 1 0 10 20 30 40 50 60 70 80 90 100 Time step (miorosec) Figure 3.2: Traction response curves (step and impulsive) of a steel bar to a step displacement of 1 0 - 6 m normal to one end. T h e impulse response was computed using eq. (3.16) but not scaled by 1 / T . A t this point we are ready to define a contact map. Definition 1 The contact map Ca for a linear elastic solid A is a set of functions {gA(x.i,t), ...,gA(xn,t)} defined at a discrete set of nodal points x \u00E2\u0082\u00AC {xi, x 2 , x \u00E2\u0080\u009E } located on the surface of A, and on the time interval t G \ ( 0 , T ] . The functions \ gA(x.j,t) correspond to the Green's function solution of eq. (3.4) for A due to a normal displacement 8(t, x \u00E2\u0080\u0094 X j ) . A typically contact response curve (for the metal bar) is shown in fig. 3.2. A s is well known from the study of differential equations, the Green's function allows us to solve the P D E for any boundary function by a convolution [23]. We pursue this technique below. 34 3.5.1 Linear Systems Theory Note that the equation for elastic bodies (3.4) is a linear P D E . This linearity allows us to consider the response of an elastic body in terms of linear systems theory [51], which will lead to a technique for using contact maps. In terms of systems theory, contact maps characterize the input-output re-sponse of a solid with no explicit reference to the internal state. Each boundary node on the body corresponds to an input 'terminal' and the corresponding tractions are output 'terminals'. In this sense, the elastic dynamics of the body constitute a black box system which we can only investigate and characterize using these 'terminals'. This approach is known as lumped systems theory and is discussed in detail in [37]. Linear systems theory uses three basic functions to analyze input-output systems, the delta function 8(t), the Heaviside step function H(t) and the one pa-rameter family of sinusoids elwt. For our present purpose, we consider only the first two. We use the T subscript to denote the discrete versions of these functions: 6T(t) = i 1 / M \" \u00E2\u0082\u00AC [ \u00C2\u00B0 ' T 1 , (3.14) 0 otherwise HT(t) = { 0 for t < 0 (3.15) 1 for t>0 T h e subscript is necessary because the discrete delta function contains an explicit time scale, in contrast to the actual delta function. T h e discrete response h(t) of a linear system to a step function also gives the the impulse or delta function response g(t) by backward differencing: 35 for integral k. Note that this relationship is exact in the sense that sampling the continuous impulse response will produce the same result( see [21] for a proof) as (3.16). Given the impulse response, we can use the properties of the delta function and the system linearity to construct the response to arbitrary inputs. Thi s is done with the discrete analogue of the convolution product, ra-1 F(nT) = J2 9((n - k)T)y(kT)T, (3.17) k=0 where y{kT) is the actual system input. We incorporate this very useful result into our time stepping algorithm in the next section. 3.6 Simulation Using Contact Maps Once the contact maps have been determined for a number of bodies, we naturally would like to do simulations in which objects may collide with each other. A full featured simulation environment is discussed in Appendix B , but we need only con-sider an isolated collision event here. For our purposes, the simulation engine needs to provide us with accurate inter-object distances in the sense of minimum transla-tion distances (MTDs see [13]). We define the MTD in terms of the translational configuration space obstacle ( T C S O ) . Definition 2 TCSO{A, B) = {b - a\a eAandb\u00C2\u00A3 B}. This is also known as the Minkowski difference of the two sets A and B. Next we define the distance between to sets as, Definition 3 d(A, B) = ajtleB\a - b\, 36 for some norm | \u00E2\u0080\u00A2 |. T h e MTD can now be defined, Definition 4 MTD(A,B) = d(0,TCSO(A, B)), where 0 is the set containing only the origin. T h e MTD between two objects is the minimum distance that one object must be translated for the two objects to just touch [13]. We discuss algorithms for extremely fast calculation of this distance in Appendix B . Once it has been determined that two objects are colliding, we must find the closest nodes at which we know the contact map on each object. This can be quickly done if the map is stored in a topologically consistent data structure. After the relevant maps have been found, we integrate the collision using a convolution technique described next. 3.6.1 Convolution Integration Given the relevant step response curves, we develop a procedure for resolving colli-sions in a rapid manner. T h e idea of this approach is to consider the penetrations as occurring in discrete steps and use the impulse response to calculate the behav-ior of the solids at each step. In addition, we must ensure that the principle of action-reaction (Newton's third law) is satisfied at each step. In this algorithm, we consider the collision of two bodies A and B. Define the penetration depth e of A as, eA = (uA - uA) \u00E2\u0080\u00A2 n, (3.18) where u is the co-ordinate of a point on the undeformed boundary, u is the co-ordinate of the same point on the deformed boundary and n is the contact normal. T h e penetration depth of B is define similarly. We can think of e as the surface 37 strain of body A. Assume that contact begins when time t = 0. We compute the force at time r > 0 using a convolution product with the penetration trajectory, / ( r ) = f g(r - t)e(t)dt. (3.19) Jo In this section, we assume that the specific impulse response map g(t) has been determined by the methods described above, experimentation or by another method. Further, we assume that the simulator has selected the correct contact map based on the contact location. Now we approximate the penetration depths by a series of impulses and thus compute the force up to time tk as tt+1=J2gA(jAt)etv (3.20) 3=0 where we define fA+l \u00E2\u0080\u0094 fA(tk+i): Equation (3.20) is essentially the endpoint ap-proximation to (3.19) and has a constant error order independent of At. We use this quadrature formula only to simplify the algebra. In practice, a better approximation would be used such as midpoint or Simpson's rule. Note that if the time step At is the same as the constant T used in (3.16) we do not need to multiply (3.20) by At in the discretized equation. A similar calculation is performed for for fB(tk+i). Given the forces fA+1 and fB+1 we compute the updated positions of the two bodies using a standard O D E integrator, such as RK4 (Runge-Kutta-4, see [10]). This will give us the updated values for the boundary coordinates uA and uB. W i t h the updated boundary coordinates, we compute the total penetration depth along the contact normal as Cit+i = (\"*+! -uf + 1 ) i -n . (3.21) To determine the updated positions of the deformed boundary points, we use 38 the principle of action reaction, f*+i + fj?+ 1 = 0. (3-22) We have in the normal direction fk+i=fk+i+9A(te)4+i, (3-23) fk+i=fk+i+9B(At)*k+i- (3-24) for both A and B. A n d also, Cfc+ i -c f+ i = Cfc+i- (3-25) Equations 3.23 and 3.25 form a system in two unknowns eA+l and ek+1. Solving the system yields, A 9B(*t)6k+1-ft+1-fS+1 \u00E2\u0082\u00AC k + 1 ~ gA(At) + gB(At) ^ B _ 9A(^t)ek+1+fA+1 + fkB+1 \u00E2\u0082\u00AC f c + 1 - gA(At) + gB(At) ^ which gives the co-ordinates of the deformed boundary for the next time step. T h e time stepping algorithm is repeated until the value for the contact force has dropped to zero, at which point it is assumed that the bodies have separated. Note that the actual separation distance is not likely to be zero at the same time as the force during object restitution. In this respect, this model is similar to the bilinear force model [58] and to the traditional spring damper model [26]. Further, note that in the case where the impulse response map has no memory (i.e. the step response is constant), the above algorithm is identical to a linear spring penalty model. 39 3.6.2 Complexity T h e complexity of this contact algorithm is certainly greater than that of both the impulsive and penalty approaches (0(1)). This is due to the convolution product in eq. (3.20) whose per step cost grows linearly with the number of time steps. Although we do not pursue it in this thesis, we note that this complexity can be reduced to a constant order by replacing the convolution with the finite difference equation corresponding to the Green's function (see [51]). 3.7 Discussion First note that the time step At in the above algorithm should match the time scale introduced in the calculation of the contact response map. Since the computed contact response map will have a large peak at the first sample point, using a different time scale will introduce extra error into the algorithm. A s presented, the integration algorithm is a zero order technique. Thi s is due to the endpoint approximation to the convolution integral in (3.20). Higher order could be attained if the displacement were a continuous function. Consider the traction response to a set of basis functions in an appropriate function space. One promising possibility would be to use a compactly supported wavelet basis. Compact support would reduce the computational cost of the convo-lution product in equation (3.19). Another possibility it to use a standard Fourier basis. However, for harmonic bases it is simpler to Fourier transform the original equation (3.4) into: {c\ - cffujji + cjuijj + ~= (3.28) where u is the Fourier transformed displacement. Equation (3.28) is much easier 40 to solve than its dynamic counterpart. Numerical techniques that speed up the solution of the corresponding B E M equation, such as multigrid [31] can be applied to (3.28) in a straightforward way. T h e generalization of the contact map algorithm to arbitrary bases is be-yond the scope of this thesis. We note however that for tangential motion or three dimensional impact, a more general approach is needed than the simple scalar step function approach. T h e reason for this is that we have more than just a one di-mensional space to span. To accommodate tangential impact in 2D for example, we would need a matrix of step functions to span the space of allowable impacts. 3.7.1 Relation to Viscoelasticity T h e convolution integral (3.19) used in simulation has a direct correspondence to the constitutive relation in the theory of viscoelasticity. In that theory, stress and strain are related by a Stieltjes integral expressed in convolution form as (see [15]). This integral expresses the memory property of viscoelastic materials. A discretized version of (3.29) would yield an algorithm similar to the one in the previous section, but using the step response map see ([63]). In this theory, elastic materials can be identified as Gijki(t \u00E2\u0080\u0094 r) = EijkiS(t) which reduces (3.29) to the standard elastic constitutive equation. For a given material, the tensor Giju is determined by applying a step strain to a half space of a given material and determining the resulting time dependent stress. T h e relationship between (3.29) and (3.19) is readily apparent. However the contact map is calculated from the equations for a non-dissipative elastic solid, so the memory property can only be related to the boundary conditions. Viscoelastic (3.29) 41 theories also make the additional assumption that the magnitude of the tensor Gijki is a monotone decreasing function of time. In our case, the response can only be said to be bounded in time due to oscillations from reflected strain waves. 3.7.2 Work It is of special importance that the total work done during a collision event be positive. Negative work would result in the objects separating with higher velocity than they initially collided with. We can write the total work done as, W= (tf [* g(t-T)e{T)dTde(t). (3.30) Jo Jo To conserve energy, we must require that W > 0. Without specific assumptions about the structure of g, it is difficult to show this rigorously. In the theory of viscoelasticity, it is relatively straightforward to demonstrate positive work, because the tensor G%jki is monotone decreasing. In our case very little is known about g, primarily because of the indeterminacy of the boundary conditions in general. We have observed cases in which the total work was negative, but only for artificially constructed g's, not for contact maps resulting from B E M analysis. 3.8 Summary This chapter presented a complete algorithm for determining and using contact response maps. Contact response maps attempt to capture more of the physics of the objects using as few input parameters as possible. It is important to note that there is very little ambiguity in the determination of a contact map, since the density and Young's modulus are relatively easy to obtain for a wide variety of materials. T h e object geometry can be determined with a great deal of accuracy as well, using 42 for example the A C M E [47] apparatus. 43 Chapter 4 Numerical Experiments 4.1 Introduction In this section, we exhibit some examples computed using the contact map algorithm introduced in the last chapter. We begin by showing the consistency of the algorithm with the energy conserving harmonic penalty model, then we give examples which have a damping component. We also discuss the possibility of fitting both impulsive and force response rigid models from Chapter 2. 4.2 Consistency We begin by considering a contact map generated from a step response which is flat. A flat step response indicates that the body has no memory and behaves in a perfectly lumped fashion. This can be interpreted physically as an elastic solid with no transient behaviour and no wave properties. A l l the models of Chapter 2 belong to this category. We model the collision of two unit spheres both with identical constant step 44 Step Size Llnf Force Error 2 x l O \" 4 2.45 1 X l O - 4 1.21 5 X 10~ 5 .456 2 x l O \" 5 .240 Table 4.1: M a x i m u m force errors between contact map and harmonic penalty model. response. These step response profiles were chosen arbitrarily and do not correspond to the solution of any elastic solid models. T h e constant was AOOON/m, correspond-ing to a spring penalty model /(e, e) = ke with spring constant k = 4000A r /m. Table 4.1 clearly shows that the contact map algorithm is only of first order. In addition, our implementation uses a fixed step size 4th order Runge-Kutta [10] integrator for the O D E . Since this integrator is explicit, very small steps (on the order of must be taken to resolve the high frequency oscillations. 4.3 Gauging Energy Loss In this section, we try to get a feel for the amount of energy loss associated with a contact map and in particular how the energy loss relates to the transient portion of the response. We collide two spheres, as. before, but add successively larger transients to the step response. One of the spheres is modeled with a non-transient map as in the previous section. T h e step responses are shown in figure 4.1 and their corresponding impulse responses in 4.2. O u r principal measure of energy loss is the relationship of inward to outward velocity for the spheres as in the Newtonian definition of restitution. Since this model is frictionless, the three definitions of the restitution coefficient are equivalent. T h e computed restitution is shown in figure 4.3. Collision times for these 45 7000 6000 5000 4000 3000 2000 1000 \"i 7. r \"responseMapA2000.dat\" \"responseMapA3000.dat\" \"responseMapA4000.dat\" \"responseMapA5000.dat\" \"responseMapA6000.dat\" \"responseMapA7000.dat\" 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 Time (s) Figure 4.1: T h e step response curves used to gauge energy loss. T h e steady state behaviour is constant at 1500 N / m in all cases. 46 Derived Impulse Response 1 1 Peak 2000 Peak 3000 Peak 4000 Peak 5000 Peak 6000 Peak 7000 0.002 0.004 0.006 Time (s) 0.008 0.01 Figure 4.2: Derived impulse response curves used to gauge energy loss. 47 1000 2000 3000 4000 5000 6000 7000 Transient Peak (N/m) Figure 4.3: T h e computed restitution coefficient for each step response, as a function of the transient peak size. cases was around .080s, so most of the force was the result of the steady state, and the transient acted as a damping mechanism. We note that in these simulations, the step size for correct integration was set to 1 X 1 0 - 4 . This is quite small since it takes nearly 800 integration steps to resolve the collision in this case. One possible way to speed this up would be to reformulate the contact algorithm to handle variable step size integration. Another approach is to try and extend the algorithm itself to be 2nd or higher order accurate as discussed in the previous chapter. A third technique for speeding up this approach is to use the results of these simulations to determine coefficients for a model from chapter 2. A s an example, we consider the simple harmonic force with damping given by /(e, e) = ce + ke. Using 48 least squares on the computed trajectories from the responseMapA6000 in figure 4.1, we get k = l59l.3N/m and c = 9.0728./Vs/m, a modestly damped system. Note however, that there is no viscous damping in this system. T h e non-zero value for c is due entirely to internal elastic energy storage. Thi s approach might also benefit existing simulation systems that rely on explicit force models such as [26],[19]. Using this approach, more complex force models could be fit, such as [50]) the model of Hunt and Crossley eqn. (2.15) with a non-linear fit (e.g., Levenberg-Marqd. ) 4.4 Real Contact Maps T h e contact maps discussed in the last section were syntheticaly generated, in that they don't necessarily correspond to any real solid. In this section, we examine some maps generated using the boundary element technique discussed in the last chapter. O u r geometry will be the often studied uniform isotropic bar with rounded ends being dropped at various angles. T h e boundary of the bar is shown in figure 4.4 along with the nodal points which denote the boundary elements. T h e B E M algorithm used quadratic support functions to discretize the boundary [22]. We simulate dropping this bar by setting the normal displacement at successive boundary nodes at one end of the bar. A l l other boundary nodes are given Neumann conditions as discussed above. T h e code developed in [22] is known as QUADPLET and is generally used for B E M calculations. During our use of this code, it was observed that the algorithm was generally unstable, unless the ratio was close to 1 (where L is the average node length). B y subdividing the nodes in the original model, contact maps at the appropriate time scale were obtained. 49 0.15 0.05 -0.05 -0.1 -0.15 Figure 4.4: T h e bar used to determine contact maps 50 E 2.1 X 10uN/m P 7876 .74%/m 3 C\ 9660.m/s C2 5163. m/s Table 4.2: Bar material properties. 2e+11 Contact Response for 18 deg. Impact 1e+11 1 <2 5e+10 I--5e+10 h -1e+11 5e-05 0.0001 0.00015 0.0002 Time (s) 0.00025 Figure 4.5: Derived impulse maps resulting from 18 deg. impact. T h e material properties of this bar, and its physical dimensions are as in [55]. We summarize them here for convenience in table 4.2. T h e resulting derived impulse maps are shown in figures 4.5 through 4.12. Note that the point at time t = 0 is undefined and can be chosen arbitrarily as 0. M a p s for other angles can be determined by rounding to those given angles or by some type of interpolation. We now make a few observations on the contact response maps. T h e initial transient peak is on the same order as the Young's modulus of elasticity for the 51 2e+11 1.5e+11 Contact Response lor 30 deg. Impact 0.0001 0.00015 Time (s) Figure 4.6: Derived impulse maps resulting from 30 deg. impact. Contact Response for 42 deg. Impact 5e-05 0.0001 0.00015 Time (s) 0.0002 0.00025 Figure 4.7: Derived impulse maps resulting from 42 deg. impact. 52 Contact Response for 54 deg. Impact 0.0001 0.00015 Time (s) 0.0002 0.00025 Figure 4.8: Derived impulse maps resulting from 54 deg. impact. 1.46+11 1.2e+11 1e+11 8e+10 S 6e+10 2e+10 -2e+10 -4e+10 Contact Response for 66 deg. Impact n 1 1 r - -| r -66 deg. 2e-05 4e-05 6e-05 8e-05 0.0001 0.00012 0.00014 0.00016 0.00018 0.0002 Time (s) Figure 4.9: Derived impulse maps resulting from 66 deg. impact. 53 Contact Response for 78 deg. Impact 1e+11 8e+10 6e+10 4e+10 -2e+10 -4e+10 5e-05 0.0001 0.00015 0.0002 0.00025 Time (s) Figure 4.10: Derived impulse maps resulting from 78 deg. impact. 1,2e+11 1e+11 Contact Response for 90 deg. Impact 5e-05 0.0001 0.00015 Time (s) 0.0002 0.00025 Figure 4.11: Derived impulse maps resulting from 90 deg. impact. 54 o CO Q. E CD CO c o Q. CO CD DC o S o O 0 ) 0 0)0)0)0)0) CD 0 CD CD CD CD CD TJ T> XJ TJ T3 TJ TJ O 00 CD 1^\" CJ O 00 cn s co in co T -CM O O O CM O O O m o o o d E o o o in o i o in + CD CM + CD in + CD + CD in + CD in + CD ( L U / N ) esuodsey ;ure 4.12: Derived impulse maps resulting from step displacement at each an*. 55 material (2.1 X 1011 N/m). We expect this since in the early stages of compression, the solid should exhibit perfect Hookean response. T h e return wave time of approximately 40/is is consistent with the wave speed in table 4.2 and a bar length of 0.2m. A s the impact angle decreases from normal, the pressure wave will be reflected more times in the interior before returning to the source point. This behaviour is manifested in 4.12 by multiple peaked response for 6 = 18\u00C2\u00B0 (fig. 4.5 )and only a single wave pulse for 6 = 9 0 \u00C2\u00B0 (fig. 4.11 ). Further for angles below about 6 6 \u00C2\u00B0 (fig. 4.9 ), there is an nearly immediate reflection from the opposite side of the bar. Simulation using the maps shown in fig. 4.12 will be discussed in chapter 6, after we have reviewed friction models. 4.5 Summary In this chapter, we showed the consistency of the convolution algorithm presented in the last chapter in the case of a linear penalty model. T h e relationship of the memory length of the contact response map to the loss of energy to internal modes was exhibited for a set of artificially chosen step, response curves. T h e non-local nature of the contact algorithm makes an analytic prediction of the energy loss difficult in all but the simplest maps. It was also shown how to use this algorithm as an off-line analysis tool to calculate the coefficients of impulse and force response rigid models. This technique offers to increase the accuracy of such algorithms with minimal extra computation during a simulation. Finally we generated some actual contact maps for the rigid bar geometry. 56 Chapter 5 Friction 5.1 Introduction U p to this point in the thesis, we have assumed collisions did not involve any velocity tangential to the point of contact. However, in nearly all contact situations there is some tangential motion, which can have a significant effect on the collision resolu-tion. Forces which occur during tangential motion are largely the result of friction, however there has been some discussion in recent years of tangential restitution [8]. We begin by discussing the physics of friction and the commonly used models. Following this, we describe some impulse rigid models which try to capture frictional effects and some of their problems. Finally, we describe models which try to predict frictional forces. 5.2 Historical Friction T h e most common friction law, still widely used, was proposed by Coulomb nearly 2 centuries ago. This law is taught in most first year physics courses and serves as 57 a benchmark against which new models must be measured. Coulomb's most remarkable discovery was that friction is independent of the area of contact and only slightly dependent on velocity. T h e largest component of friction is due to the load or normal force of contact. Also if two objects are pressed together with zero relative velocity, they tend to stick even for moderate non-zero tangential force. Coulomb's law of friction is stated mathematically as, ,here \i \u00E2\u0082\u00AC [0,1] is known as the coefficient of friction. T h e coefficient n is sometimes given by two values, one for static and one for kinetic motion. Strict equality of (5.1) is required for non-zero tangential velocity (kinetic). For zero tangential velocity, the tangential friction force is constrained only to lie within the force cone defined by the inequality. T h e indeterminacy of (5.1) during sticking (zero tangential velocity) leads to non-uniqueness of solutions. This is discussed by several authors [44], [6], [7]. We also discuss this dynamic effect below. 5.2.1 Rotational Friction One immediate extension of the simple Coulomb model is for rotational friction. W h e n two bodies in contact have nonzero relative angular motion, friction will occur over the contact area A. We can write, where fj(x, y) is the coulomb friction at point (x, y) in direction j as given by (5.1). Further extensions of this are discussed in [28], [29]. I / T | < M I / A T (5.1) (5.2) 58 5.3 Tribology Modern theories of surface contact forces and effects are the subject of the field of tribology. From an engineering perspective the Coulomb model is usually sufficient, however it provides little insight into the physics of friction nor do we have a good understanding of the limits of the Coulomb model. 5.3.1 Physics of Contact Nearly all surfaces, no matter how smooth they appear to the eye contain a rough terrain at micro ( 1 0 _ 6 m ) t o nano ( 1 0 - 9 m ) scales. T h e peaks of this terrain are called asperities. Asperities can be observed readily with a profilometer and are usually characterized by their statistical properties [30]. Asperities are the arbiters of contact force between two bodies. In addition to microscale roughness, all solids have chemical properties and some combinations can significantly affect contact forces. For example, coefficient of friction is known depend strongly on the level of oxidation in the contact area. Tabor [60] describes placing two indium metal spheres into contact, when highly cleaned they adhere, but when exposed to air before contact, they slide easily. More complex models of friction can be constructed, based on modern knowl-edge of tribology. One example is a proposal by Pollak [49] which includes adhesive forces, \fr\ = fi(P + 2nRw + y/6nRwP+9n2R2w2), (5.3) where P is the surface pressure, R is the radius of curvature near the contact zone, w is an adhesion coefficient. T h e first term of (5.3) represents the load, the second is the adhesion and the third expresses the variance of adhesion with load. 59 One of the chief problems with models such as (5.3) is that for real bodies, it is difficult to obtain the correct coefficients and once obtained, they may easily change due to environmental effects. In light of this, we are forced to agree with Chatterjee [14] who states that collision models are currently so poor that we can hardly err by using the Coulomb model. 5.4 Impulsive Collisions with Friction When using an impulse response rigid contact model, the addition of friction dra-matically effects model consistency. T h e three definitions of restitution given above are all equivalent in non-frictional collisions but predict dramatically different out-comes when friction is included. Following an insightful paper by Routh [35] on the two dimensional case, almost no research was done in this area until the work of Keller [39], who considered the full 3D contact problem. Since [39], a large number of papers have been published on this topic, see for example [8], [56], [57], [64]. A n excellent review can be found in [14]. Stronge noted in [56] that the standard definitions of restitution, either (2.8) or (2.9) fail to conserve energy when the tangential velocity is zero at some point dur-ing the collision (i.e. contact sticking). This problem is resolved by using Stronge's work related restitution definition. 2 Wc J * \" p \u00E2\u0080\u00A2 ndt wr~ /// p . n d t ' (5-4j where tmc is the time of maximum compression, and tj is the collision termination. Chatterjee [14] analyzes in detail the problem of energy conservation for three dimensional collisions. T h e chief tool in his analysis is the energy ellipse, an 60 ellipsoidal region in impulse space that bounds all energically consistent collision impulses. T h e one parameter impulse models in chapter 2 can be thought of as defining curves within this ellipsoid. T h e addition of tangential forces (and tan-gential coefficients) transform this line into a surface. T h e difficulty encountered by many authors constructing impulse response rigid models is to constrain the parameters of their models to the interior of the energy ellipsoid. In addition to the energetic consistency problem, the so called stick-slip phe-nomenon when the tangential velocity vanishes during collision, introduces non-uniqueness. T h e direction of the tangential force when objects are sticking is not predicted by the Coulomb model. We discuss a possible solution to this problem in the next section. Finally, we note that a complete simulation environment has been written which uses the Stronge definition of restitution along with Coulomb friction. T h e system, called Impulse is the subject of Mirtich's P h D thesis [44]. 5.5 Solution Existence and Uniqueness In this section, we begin by formulating the impulsive contact problem as a linear complementarity problem or L C P . L C P ' s are widely known in the mathematical community and this standard formulation will expose the core issues for consistent solution of frictional collisions. We consider a situation in which several bodies are colliding at n points. B y consideration of the non-penetration constraints at each contact point we can write Newton's law, a n = Af n - b, (5.5) 61 where A describes the relative inverse inertias of the bodies, f\u00E2\u0080\u009E are the contact forces, b are the external accelerations and a n are the resulting outward normal accelerations. For impulsive collisions, we use the same equation but define a n as the resulting velocity, b as the initial velocity and fn as the impulse. For non-frictional collisions, it can be shown that A is positive semi-definite (see [5], [44]). Kinematic consistency is imposed on the problem by requiring nonnegative normal accelerations and forces, a\u00E2\u0080\u009E > 0 fn > 0. (5.6) In addition we require the complementarity condition, a\u00C2\u00A3f n = 0, (5.7) which states that the force is zero in non-constrained configurations. T h e equations 5.6 and 5.7 together define a L C P . Solutions to L C P ' s can be obtained in expected polynomial time if the matrix A is positive semi-definite (PSD) by reformulating the problem into a quadratic program ( Q P ) , see [5], [48]. T h e existence of solutions to the L C P is guaranteed in the case when A is P S D and the forces b lie in its column space, since the corresponding Q P is convex. Uniqueness of solutions can be shown only in the restricted case when A is positive definite (see [18], and [11]). A simple example of this situation is a four legged chair falling onto a flat surface. T h e distribution of constraint forces among the four legs is non-unique (see [46]). T h e addition of Coulomb forces further reduces the exact solubility of this problem. Since frictional forces occur in the tangent plane to the point of contact, we must add a tangential term Qf4 in eqn. (5.5. We then construct an L C P , a = Bf + b, (5.8) 62 where a is the resulting accelerations, f are the applied forces (normal and tan-gential), b are the body forces and B contains the columns of both A and Q. This resulting matrix B is no longer P S D . T h e consequence of this fact is the non-convexity of the corresponding Q P . T h e solution of non-convex Q P ' s is well known to be NP-hard , if a solution even exists. BarafF recommends using the Q P algorithm of Lemke [41] with modifications to ensure termination in the non-convex case. We note that the non-convex problem expressed in Q P form is a problem of constrained global optimization. This is currently an active area of research and new techniques and methods which perform better than the modified Lemke algorithm may soon exist. 5.6 Contact Sticking In this section we consider the case in which sticking occurs during collision inte-gration. This introduces additional non-uniqueness into the problem because static Coulomb friction does not define a direction for the resulting forces. Thus when sticking occurs during a collision and the Coulomb model indicates a resumption of sliding, it is unclear how this transition can be accomplished consistently. In a series of papers, Bhatt and Koechling [6], [7] consider this problem in detail. It is also discussed at length in [44]. Their starting point is the impulsive version of (5.5), v = Ap + v 0. (5.9) Because they are concerned primarily with impulsive collisions, the authors express the collision dynamics in terms of the normal impulse pz instead of the collision time t. This change of variables gives additional insight into the meaning of the coefficient 63 of restitution e. Since collision integration termination is given directly by e it can be thought of as defining the flow length in state space of the corresponding dynamical system. In the tangent plane to the collision point, the Coulomb model can be written, dPx dFx /ivx dPz dFz ^jv2 + V2 dPy dFy /IV y dPz dFz (5.10) (5.11) \Jvl + vy where // is the coefficient of friction, vx and vy are the tangential velocity compo-nents. We insert (5.10) into the derivative of (5.9) to generate a first order dynamical system in the tangent plane, dVX fiVX JlVy -jw = -an . - ai2 y + \u00C2\u00ABi3 (5.12) dVy flVX jlVy dP = ^ I 2 \u00E2\u0080\u00A2 2 \" ^ / 2 , 2 + ( 5 J 3 ) z \JVx + Vy \JVx + Vy B y carefully considering the possible values of the parameters in (5.12), the authors show that the parameter space is a compact subset of -ft3, * = ^ f A = ^ , * = ^ . (5.14) \u00C2\u00AB13 \u00C2\u00AB13 .013 We can express the stick-slip condition in terms of the flow parameters as A < 7 ^ \" 6>- (5-15) Note that in (5.12), the flow is undefined in the case of sticking vx = vy = 0. The singularity can be moved to infinity with the singular transformation, dr = ^ dPz. (5.16) Iv2 + v\ 64 Finally Bhatt and Koechling transform into polar co-ordinates, \u00E2\u0080\u0094\u00E2\u0080\u0094- = vr(\u00E2\u0080\u0094A cos 2 vg \u00E2\u0080\u0094 (j>sm2 vg + cos vg + Ssin v$) (5-17) dr \u00E2\u0080\u0094\u00E2\u0080\u0094 = (A \u00E2\u0080\u0094 ) cos vg sin vg + 6 cos vg \u00E2\u0080\u0094 sin vg. (5.18) dr Obviously this system has a fixed point at vr = 0. Invariant flow directions are obtained from the real roots of ^ = 0. T h e authors show that the conditions in (5.15) exactly correspond to the situation in which there is a unique invariant flow direction on which > 0 (i.e. outflow). This effectively resolves the stick-slip indeterminacy problem for single contact configurations. We note that although this approach is concerned with impulsive collision models, a similar analysis can be carried out for force response rigid models, using eqn. (5.5). In this way, the results of this section could be incorporated into a simulator using either type of collision resolution model. 5.7 Force Response Rigid Collisions with Friction Although much of the preceding discussion applies in principle to force response rigid models, it is insightful to consider the simplification provided by assuming a particular model. In this section, we consider using the 3D harmonic penalty model to resolve a frictional collision. T h e algorithm developed here is presented as part of a simulation system in [26], [27]. We begin with the penalty equation in 3D, fi(t) = K,-(q,- - P i ) + C;(q; - P i ) , (5.19) where p represents the rigid body motion, q is the deformed motion so that their difference can be thought of as the surface strain. Subscripts refer to the two bodies 65 in collision. T h e vectors p and q are in Sft3 and we assume that the third component points in the direction normal to the contact point. T h e penalty parameters are given by two positive definite rank 2 tensors K and C representing spring force and viscous damping respectively. We assume that the third component of f is normal to the collision point. Coulomb friction is incorporated with the quadratic form, f r Q f < 0, (5.20) where Q = diag(l, 1, -fi2). T h e Coulomb model also gives a direction to f in the tangent plane, q 2 - q i = A P f . (5.21) For some positive real value of A. Note that P = diag(l, 1, 0). Goyal uses a first order Taylor expansion to generate the next value of q,-, q,-(i + A t ) = qi + Atq; ( t + A t ) . (5.22) Note that this formula corresponds to a backward Euler numerical integration scheme. Inserting (5.19) into (5.22) and doing a little algebra, we arrive at the familiar equation (A + A P ) f (t + A t ) = a, (5.23) where A and a are functions of known quantities. Further, A is positive definite because of the positive definiteness of K and C. A t this point the only remaining unknowns are A and f. We proceed by considering the case A = 0 (sticking), giving the force as f* = A _ 1 a . If the Coulomb condition f * T Q f * < 0 is satisfied, we set f(t + A t ) = f * . (5.24) 66 If the Coulomb condition is not satisfied then the objects must be slipping. In this case, we can use the quadratic form (5.20) to determine A, Which we solve to get the positive root A*, that gives the next force value as, This algorithm produces unique contact forces while incorporating the Coulomb friction model. Thi s is achieved primarily because of the positive definite assump-tion on the contact parameter tensors K and C which forces the problem to be convex. A s mentioned in chapter 2, the authors of this algorithm give no straight-forward method for determining the values of these tensors for real bodies. Final ly we note that because of the tensor nature of the contact parameters, it is possible to incorporate tangential spring and viscous forces in addition to friction. 5.8 Discussion In this chapter, we have considered the complexity added to contact models by the addition of frictional forces. Because of the nature of the Coulomb model, both existence and uniqueness can be lost. Impulsive models must be especially careful when adding friction since conservation of energy is not guaranteed in this case. Force response rigid models, in particular the harmonic model can be used to generate closed form algorithms with the proper conditions on the contact model. ( ( A + A P ) \" 1 a ) T Q ( ( A + A P ) - x a ) = 0. (5.25) f(t + A t ) = ( A + A * P ) _ 1 a . (5.26) 67 Chapter 6 Contact Map Simulations with Friction 6.1 Introduction In this chapter, we hope to establish the physical consistency of the contact map algorithm by comparison with experimental data. We will focus primarily on a two dimensional collision event for which we have good observational data. Unfortu-nately, it is very difficult to measure directly the physical effects which occur during a contact event. In addition, as discussed in chapter 2, contact mechanics is affected by a large number of variables many of which are difficult to isolate and quantify. For these reasons, there is a shortage of applicable experimental data with which to compare our algorithm. 68 Figure 6.1: Bar experiment configuration. 6.2 Previous Work O u r primary source of data will be a paper by Stoianovici and Hurmuzlu [55]. T h e authors considered the impact of a rigid steel rod onto a massive foundation. T h e rod was dropped at various orientations and with differing impact velocities. D a t a collected during the experiments consisted of pre and post impact tip velocities in the normal and tangential direction. From this data, the authors in-ferred the friction model using the impulse ratio fi, M = ? ^ ^ , (6.1) y - y 69 where x is the tangential velocity, y is the normal velocity and superscripts indicate pre and post collision. T h e experimental data indicated that the impulse ratio was zero at normal impact (90'), increased linearly with decreasing angle and then became constant below some slip angle (roughly 81) . T h e linear portion of the impulse ratio is consistent with static friction and the constant portion with kinetic friction. In all experimental setups below the slip angle, sticking did not occur at any point during the collision. We know that in this situation, the three definitions of the coefficient of restitution are equivalent, hence the authors calculate restitution using the Newtonian definition (2.8). T h e coefficient of restitution was.observed to decrease with decreasing angle to a critical angle (roughly 63'. for a 200mm bar). Below the critical angle the coefficient became constant at a higher value. T h e authors proceeded to construct a numerical model of the experiment by decomposing the rigid bar into n rigid segments each connected by two sets of springs. T h e rigid foundation was modeled using a Hunt and Crossley type penalty model. We present a similar model below. Mode l parameters were calculated by fitting the solution curves with a single set of experimental data. 6.3 Rigid Body Simulations Before we directly compare simulation with contact maps to the results reported in [55], it was considered important to first model the collision as faithfully as possible without including elastic effects. This was done for two reasons, first we wanted to get a good idea of how sensitive the physics was to elastic waves and secondly we wanted a check for the portions of the data in [55] which were indicated to have 70 minimal dependence on the presence of elastic waves. O u r model is based on the setup depicted in figure 6.1. O u r equations of motion, written in first order form are, = x2{t) (6.2) i 2 ( t ) = ^fx(t) (6.3) yi(<) = 2/2 (0 (6.4) Mt) = \u00C2\u00B1fy(t) (6.5) 0i(t) = 62(t) (6.6) Ut) = ^ - ( / ^ O c o s ^ W J - ^ ^ s i n ^ i C t ) ) , (6.7) where m is the bar mass, Ixx is the moment of inertia of a uniform bar of length L rotating about it's geometric center. We also define the tip velocities as xtip = \u00C2\u00B1 2 - ^ s i n ( 0 i ) 0 2 ( 6 - 8 ) Vtip = 2/2 + ^cos(t91)t92. (6.9) T h e normal contact force is given by fy(t) = (K - Bytip)ytip. (6.10) Note that this is a Hunt and Crossley type penalty force and as such incorporates a velocity dependence. In one dimension, this dependence is expressed in [32] with relation to the kinematic coefficient of restitution as e = 1 + avi, (6-11) where a = | j | and the initial normal velocity V{ is assumed negative. Finally, the tangential contact force fx is given by standard Coulomb friction with static coefficient fis and kinetic coefficient [ik-71 L .2 m K 5.5e7 N / m B 2.2e7 N s / m Ms 0.1 0.075 Table 6.1: Experimentally determined collision parameters. 6.4 Model Parameters and Observations T h e parameters K, B, fis and /J^ were determined experimentally in the paper. Values are listed in table 6.1. F r o m the given values of K and B we compute a as .24s /m. In [32], the range of validity of their model is estimated as .08 \u00E2\u0080\u0094 .32. We may thus feel secure in using this contact model in this situation. We present the simulation results in figures 6.2 to 6.5. We first observe that there appears to be negligible angular dependence for the post impact tip velocity. This is due to the local nature of the collision in this model. We expect this because there is no global coupling in the absence of elastic waves. In this experimental setup, we define the critical angle as the impact con-figuration in which the final center of mass normal velocity remains negative. In the context of [55], this is the angle beyond which multiple bounces will be needed to resolve the collision. From figure 6.5, we can clearly see this is in the 60' range. This corresponds to the figure of 64.3' given in [55], so we conclude that this effect is mainly a function of the impact geometry and not of elastic waves. In addition we note that as expected from the penalty model, the energy loss and normal tip restitution are functions of impact velocity. This is not what 72 T i i i : i j o> d co d d co d in d d co d CM d UOJIOBJJ A6jau3 lejoi UOJSJHOQ isod Figure 6.2: Post impact total energy. 73 Figure 6.3: Impact durat ion. 74 co co co co E E E E i - CM CO Tt O CD O CO O 1^ o CD <-. \u00C2\u00AE L O S < _ o o CO Q. E o co o CM -J o LO CM CM LO LO d (S/LU) A l p O j O A d j l I B L U J O N UO|S!||00 } S O d Figure 6.4: Post impact normal t ip velocity. 75 o . . 1 \ \u00E2\u0080\u00A2 l\ \u00E2\u0080\u00A2in to. w co \ T - CM'CO \"* N CM lO m CM m CM (S/LU) AipopA SSBIAI jo JOJUOO JBLUJON UOJSHIOO isod Figure 6.5: Post impact normal center of mass velocity. 76 [55] state nor appear to observe experimentally. It seems unlikely that this velocity dependence is affected by the presence or absence of elastic waves. B y removing the damping from the foundation and since the only other energy loss is via friction, it is possible to use this model to estimate the energy lost solely to friction. This was done for the bar for a variety of angles, the results are plotted in figure 6.6. T h e fractional energy loss to friction also exhibits strong angular dependence, with the distinct peak at the critical angle near 65'. We discuss the importance of this type of energy loss further in the next section. 6.5 Contact Map For contact map simulations, we implemented the algorithm discussed in chapter 3. T h e contact maps used for the bar were generated using the 2D B E M software discussed in [22]. For the rigid foundation, we used the correspondence between the penalty spring constant K above and the step response discussed in chapter 4. Unfortunately this approach makes it difficult to incorporate damping into the foundation so we simply used a flat step response. We only calculated collisions at the nodal points used in the B E M calculation. This eliminated the need for collision detection and element following during contact. A s a result the following graphs have data points at impact angles of 90, 78, 66, 54, 42, 30 and 18 degrees. In figure 6.7, we present a plot of the kinetic energy during the collision as a function of time for each impact angle. T h e final states are summarized in figure 6.8 where we see the kinematic coefficient of restitution, presented as an energy fraction. Clearly there is a peak in the energy loss near the critical angle of 9C = 66'. Thi s 77 J I L 00 CD t CM CD 00 CD O) O) O) OJ Q 00 00 d o d o d d UOjPBJj ABJ9U3 Figure 6.6: Post impact energy fraction due to friction only. 78 peak is clearly absent from figure 6.2 so we may conclude that this effect is related to the presence of elastic waves. Impact durations are shown in 6.9. Comparing these to the curves shown in 6.3 we see very similar results. It appears that as with the linear harmonic penalty model, the impact duration is affected most strongly by the stiffness parameter and very little by the damping. This effect is also noted by [55]. T h e energy lost only to elastic waves was obtained by subtracting the energy lost to friction from fig. 6.6. This data is shown in fig. 6.10. We observe a monotonic decrease in the elastic energy with decreasing impact angle. Intuitively, a normal impact should retain more wave energy because the wave front and the contact region are coincident less often than for a lesser angle. For small angles, the wave energy is more uniformly distributed in the bar and thus has more opportunity to affect and be affected by the contact region. In addition, the collision duration curve is very similar to the elastic wave energy loss, indicating that the total energy transfer is related to the time spent in contact. In [55], the authors further discussed the effect of the multiple collisions which will occur for impact angles below 6C (due to the still negative center of mass velocity). This effect requires that the bar is already vibrating on impact for the second and subsequent collisions. Unfortunately a straightforward implementation of the Green's function model requires that the initial conditions of impact be as-sumed zero. This makes it difficult to compare the contact map approach to the remaining data in [55]. 79 UOIJOBJJ A6J9U3 Figure 6.7: Tota l kinetic energy during collision of bar for various impact angles. 80 (s) eoiji Figure 6.9: Impact duration for various impact angles. 82 (%) A6J9U3 0 | i S B | 3 Figure 6.10: Percentage of energy lost to internal vibrations. 83 6.6 Conclusion In this chapter, we have shown the effectiveness of the contact map approach to the simulation of rigid bodies. The predicted energy loss due to elastic waves qual-itatively matches the observations of [55] with values which are also in the same range. In addition, the data generated using the contact map approach took only a fraction of a second to compute on a Pentium microprocessor using an unoptimized implementation of the algorithm in chapter 3. We believe that with a proper implementation, the contact map algorithm could be fast enough for use in real time simulation. This would allow more accurate simulation with a very minor penalty in computational overhead. 84 Chapter 7 Conclusions In this thesis, we have reviewed the state of the art in rigid body contact models. Both local and global models were presented, but it was shown that local models suffer from a very narrow range of applicability. However global contact models present the greater difficulty of large computational complexity. This thesis proposed a global contact model for point contacts which are of finite duration. T h e contact map model is this approach, which achieves accuracy by solving the wave equation subject to the actual geometry, and achieves speed by caching the solutions. Comparison to a similar model and experimental data showed that this approach behaves correctly, although it predicts different fractions for energy lost to elastic waves. In addition, we considered the problems of existence and uniqueness pre-sented by the inclusion of friction in any contact model. This is still an active area of research. T h e contact map algorithm presented in this paper can already be used in simple applications requiring only normal contact. For example in the design of 85 industrial part feeding machines, simulations requiring hundreds or thousands of contacts per second are necessary. While this algorithm may not give the minimum 30 frames/sec. it will certainly operate at interactive speeds. Other examples include collisions of objects which are decomposable into spherical sub-objects. In this case simple point friction could be used as in chapter 6. We note finally that a full implementation of contact response maps for three dimensions, including friction would be of great value both to researchers and to applied engineering industries. Extremely rapid and accurate contact calculations along with fast collision detection would allow system designers to use automatic optimization tools and have significant feedback starting from the earliest design steps. T h e result of all this would be improved efficiency and reliability in many C A D / C A E designed machines. 86 Bibliography [1] Brebbia C . A . and Dominguez J . Boundary elements : an introductory course. Computational Mechanics, New York, 1989. [2] H . Antes and P. D . Panagiotopoulos. The Boundary Integral Approach to Static and Dynamic Contact Problems. Birkhauser, Basel - Boston - Berlin. , 1992. [3] U . Ascher, D . K . Pai , and B . Cloutier. Forward dynamics, elimination methods, and formulation stiffness in robot simulation. International Journal of Robotics Research, 16:6, December 1997. [4] Owe Axelson. Iterative Solution Methods. Oxford University Press, Oxford U K , 1994. [5] D . Baraff. Dynamic Simulation of Non-Penetrating Rigid Bodies. P h D thesis, Cornell University, 1992. [6] V . Bhatt and J . Koechling. Partitioning the parameter space according to different behaviors during 3d impacts. ASME Journal of Applied Mechanics, 62:740-746, 1995. [7] V . Bhatt and J . Koechling. Three-dimensional frictional rigid-body impact. ASME Journal of Applied Mechanics, 62:893-898, 1995. 87 [8] R. M . Brach. Rigid body collisions. ASME Journal of Appied Mechanics, 56:133-138, 1989. [9] G . Burdea. Force and touch feedback for virtual reality. John Wiley, New York, 1996. [10] R. L . Burden and J . D . Faires. Numerical Analysis. P W S , Boston, 1993. [11] Ferris M . C . and Pang J . S. Engineering and economic applications of comple-mentarity problems. SIAM Review, 39(4):669-714, 1997. [12] C A D S I . Dads - dynamic analysis design system. Technical report, 1998. http:/ /www.cadsi .com/dads.html. [13] S. Cameron. Enhancing gjk: Comput ing minimum and penetration distances between convex polyhedra. In IEEE International Conference on Robotics and Automation, pages 3112-3117, 1997. [14] A . Chatterjee. Rigid Body Collisions: Some General Considerations, New Col-lision Laws, and Some Experimental Data. P h D thesis, Cornell University, 1997. [15] R. M . Chritensen. Theory of Viscoelasticity An Introduction. Academic Press, New York, 1982. [16] V R M L Consort ium. T h e virtual reality modeling language. Technical Report 14772-1:1997, I S O / I E C , 1997. h t tp : / /www.vrml .org /Spec i f l ca t ions /VRML97/ . [17] Sense8 C o r p . World toolkit. Technical report, Engineering Automation Inc., 1998. http://www.sense8.com. 88 [18] R. W . Cottle, J . S. Pang, and R. E . Stone. The Linear Complementarity Problem. Academic Press, Boston, 1992. [19] J . F . Cremer. An Architecture for General Purpose Physical System Simulation - Integrating Geometry, Dynamics and Control. P h D thesis, Cornell University, 1989. [20] Lewis A . D . and Rogers R. J . Experimental and numerical study of forces during oblique impact. Journal of Sound and Vibration, 125:403-412, 1988. [21] D . Delchamps. State Space and Input-Output Linear Systems. Springer-Verlag, New York, 1988. [22] J . Dominguez. Boundary Elements in Dynamics. Elsevier Science, Essex, 1993. [23] P. DuChateau and Zachmann D . Theory and Problems of Partial Differential Equations. M c G r a w Hil l , New York, 1986. [24] D . Flanagan. Java in a Nutshell. O'Reilly, Sebastapol, 1997. [25] W . Goldsmith. Impact: The theory and physical behavior of colliding solids. Edward Arno ld , London, 1960. [26] S. Goyal , E . Pinson, and F . Sinden. Simulation of dynamics of interacting rigid bodies including friction i: General problem and contact model. Engineering with Computers, 10:162-174, 1994. [27] S. Goyal , E . Pinson, and F . Sinden. Simulation of dynamics of interacting rigid bodies including friction ii: Software system design and implementation. Engineering with Computers, 10:175-195, 1994. 89 [28] S. Goyal , A . Ruina, and J . Papadopoulos. Planar sliding with dry friction part 1. limit surface and moment function. Wear, 143:307-330, 1991. [29] S. Goyal , A . Ruina, and J . Papadopoulos. Planar sliding with dry friction part 2. dynamics of motion. Wear, 143:331-352, 1991. [30] J . A . Greenwood. Problems with surface roughness. In Fundamentals of Fric-tion: Macroscopic and Microscopic Processes, pages 57-76, 1991. [31] W . Hackbush. Integral Equations: theory and numerical treatment. Birkhauser Verlag, Basel, 1995. [32] K . H . Hunt and Crossley F . R . E . Coefficient of restitution interpreted as damp-ing in vibroimpact. ASME Journal of Applied Mechanics, 42:440-445, 1975. [33] Mechanical Dynamics Inc. A d a m s simulation package. Technical report, 1998. http: / /www.adams.com. [34] Canny J . Collision detection between moving polyhedra. IEEE Trans. PAMI, 8:200-209,1986. [35] Routh E . J . Dynamics of a System of Rigid Bodies. Dover, New York, 1960. [36] K . L . Johnson. Contact Mechanics. Cambridge University Press, New York, 1985. [37] T . Kai lath . Linear Systems. Prentice-Hall, New Jersey, 1980. [38] H . Keller, H . Stolz, and T . Brunl . Aero - a physically based simulation and animation system. Technical report, T h e University of Western Austral ia , 1997. http:/ /www.ee.uwa.edu.au/ braunl /aero/ . 90 [39] J . B . Keller. Impact with friction. ASME Journal of Applied Mechanics, 53:1-4, 1986. [40] N . Kikuchi and J . T . Oden. Contact Problems in Elasticity: A Study of Vari-ational Inequalities and Finite Element Methods. S I A M , Philadelphia, 1988. [41] C . E . Lemke. Bimatrix equilibrium points and mathematical programming. Management Science, 11:681-689, 1965. [42] M . C . L i n . Efficient collision detection for animation and robotics. P h D thesis, U C Berkeley, 1993. [43] C h . Lubich, Nowak U . , Pohle U . , and Engstler C h . M e x x - numerical software for the integration of constrained mechanical multibody systems. ZIB Preprint SC92-12, 1992. [44] B . Mirt ich . Impulse-based Dynamic Simulation of Rigid Body Systems. P h D thesis, U C Berkeley, 1996. [45] R. L . Mullen and J . J . Rencis. Iterative methods for solving boundary element equations. Computers and Structures, 25(5):713-723, 1988. [46] Ivanov A . P. O n multiple impact. J. Appl. Maths. Mechs., 59(6):887-902,1995. [47] D . Pai . Acme testbed. NSERC Research Grant Proposal, 1997. [48] F Pfeiffer and C Glocker. Multibody Dynamics with Unilateral Contacts. Wiley, New York, 1996. [49] H . M . Pollak. Surface forces and adhesion. In Fundamentals of Friction: Macro-scopic and Microscopic Processes, pages 77-94, 1991. 91 [50] W . Press, Teukolsky S., and Vetterling W . Numerical Recipies in C. Cambridge University Press, New York, 1993. [51] J . G . Reid. Linear System Fundamentals. M c G r a w Hil l , New York, 1983. [52] A . Signorini. Sopra akune questioni di elastostatica. Atti della Societa Italiana per is Progresso delle Scienze, 1933. [53] I. L . Singer and H . M . Pollock. Epilogue. In Fundamentals of Friction: Macro-scopic and Microscopic Processes, pages 569-588, 1991. [54] H Sowizral, K Rushforth, and M Deering. Java 3d api specification. Technical report, Sun Microsystems, 1998. http:/ /java.sun.com. [55] D . Stoianovici and Y . Hurmuzlu. A critical study of the applicability of rigid-body collision theory. ASME Journal of Applied Mechanics, 63:307-316, 1996. [56] W . J . Stronge. Rigid body collisions with friction. Proc. R. Soc. Lond. A, 431:169-181, 1990. [57] W . J . Stronge. Unraveling paradoxical theories for rigid body collisions. ASME Journal of Applied Mechanics, 48:1049-1055, 1991. [58] W . J . Stronge. Planar impact of rough compliant bodies. Int. J. Impact Engng, 15(4):435-450, 1994. [59] W . J . Stronge. Theoretical coefficient of restitution fori planar impact of rough elasto-plastic bodies. Impact, Waves, and Fracture, 205:351-362, 1995. [60] D . Tabor . Friction as adissipative process. In Fundamentals of Friction: Macro-scopic and Microscopic Processes, pages 3-24, 1991. 92 [61] S. Takahashi and Brebbia C . A . Elastic contact analysis with friction using boundary elements flexibility approach. In M . H . Al iabadi and Brebbia C . A . , editors, Computational Methods in Contact Mechanics. Elsevier Science, Essex, 1993. [62] S. P. Timoshenko and Goodier J . N . Theory of Elasticity. M c G r a w - H i l l , New York, 1970. [63] C . Ullrich and D . K . Pai . Contact response maps for real time dynamic simu-lation. In IEEE International Conference on Robotics and Automation, 1998. [64] Y . Wang and M . Mason. Modeling impact dynamics for robotic operations. In IEEE International Conference on Robotics and Automation, pages 678-685, 1987. [65] J Wernecke. The Inventor Mentor. Addison-Wesley, 1995. 93 Appendix A Boundary Element Review This appendix discusses the boundary element method and it's application to the solution of partial differential equations. For simplicity, we use the Laplace equation Au = 0 on some fixed domain Q as the reference example for the technique. A . l Boundary Elements and Finite Elements Before we discuss the details of the boundary element method ( B E M ) , we com-pare some of it's properties with the popular finite element method ( F E M ) solution technique. Both B E M and F E M are integral equation solution procedures for partial differential equations ( P D E ' s ) . Both techniques rely on Green's integral theorem and both produce weak solutions. T h e F E M method solves a P D E by discretizing .the domain over which the solution is defined whereas the B E M approach discretizes only the boundary of the domain. In this sense, the B E M technique always has a lower dimensionality than the F E M procedure for the same problem. In the solution of an F E M problem, a 94 great deal of symmetry is always present, and all equations are of a local nature. T h e B E M solution procedure presents no special symmetry and involves global equations. A . 2 B o u n d a r y Integral E q u a t i o n s A.2.1 Fundamental Solutions To construct a B E M solution for a given P D E , we first find the so called fundamental solution (or Green's function). For the sample equation, V u = 0, the corresponding fundamental equation is, Ato(x) - f5(x!' - x) = 0. ( A . l ) For a given fixed point x \ In this case, the fundamental solutions are just the harmonic functions, to(x) = ^-log(r) intk2 (A.2) t\u00C2\u00BB(x) = m & 3 , (A.3) where r = |x* \u00E2\u0080\u0094 x|. A . 2.2 Residual Now consider the weighted residual, wAu = 0 over the problem domain Q, f wAudQ = 0. (A.4) Jo. Given an approximate solution u \u00E2\u0080\u0094 u', (A.4) gives a measure of the approximation error. T h e choice of w as a weight allows us to make direct use of Green's integral theorem, / wAu + Aw -VudQ= f w^dT. (A.5) Ja Jr on 95 Substituting (A.4) into (A.5) gives, du 0 = / w^dT- f Vw-Vudtt. (A.6) J r on Jn A n d applying Green's integral theorem again to the second integral leads to, du dw 0 = I w^ - u^-dT + / uAwdQ. (A.7) Jr on on Jn Now by ( A . l ) for points in the interior of Q, this reduces to, f du dw . . u = / w- u\u00E2\u0080\u0094dT. (A.8) Jr on on Which contains only boundary integrals. For points on the boundary, slightly more analysis is required (refer to [1] for details). T h e net result is du dw , \ , \ f du d m ., . ^ U ^ = lrWd-n-U^ ^ where in 3?2, c(x) = for x E f i \ for x e r (A.10) for corner of angle a. Equation (A.9) represents the solution to the original Laplace equation in terms of integrals over the boundary of the region Q. A t this point, given a set of boundary conditions the problem can in theory be solved by integration. \u00E2\u0080\u00A2 A.3 Numerical Considerations In practice, the integrals in (A.9) are not soluble by any analytic method. This occurs for even moderately complex boundary geometries or for non-trivial boundary conditions. Also, in engineering, most of the boundary information will be available only numerically. Thus in this section, we summarize the numerical solution of the integral equation (A.9). 96 A.3.1 Boundary Elements We begin by dividing the boundary V into disjoint subsets where, N T=\jTj. ( A . l l ) N ' du 8w This reduces (A.9) to, \u00C2\u00AB M . M = E j ! I ' M - ' K * - (A'12) Now we choose a basis set ,\u00E2\u0080\u00A2 (eg- the standard linear elements of F E M ) such that, Uj = ^fauji (A.13) i fc = ^ f = !>\u00C2\u00AB\u00E2\u0080\u00A2? ;\u00C2\u00AB\u00E2\u0080\u00A2\u00E2\u0080\u00A2 ( A - 1 4 ) i T h e choice of basis functions is an arbitrary component of this method. Typical ly the same bases are used for the geometry, Uj and qj. However this choice can be optimized for specific problems (see [2] or [1]). W i t h a basis, the B I E becomes, N N r dw c(x)w(x) = J2J2qji / i wdT - Yl uJi / (t>i-^-dT. (A.15) .7=1 t J l \" j i = l i J V J 0 7 1 There are several ways to generate a set of linear equations from (A. 15), but we discuss only the collocation method. For a detailed discussion of the other approaches, see [31]. In the collocation method, we choose a set of L boundary nodes at which we require (A.15) to be satisfied exactly. This gives, N Q N cm + ^ Y uoi I fa-p-d? = Y Y / tw^i (A. 16) j=i i Jrj o n j = 1 i JTj 97 where the subscript I indicates evaluation at x/. T h e result of this is a system of L linear equations, A u = B q . (A.17) Given a well posed problem, the system (A.17) can be written in standard form as A x = b . Where A is a rank L dense matrix with no special symmetry. T h e primary difficulty in arriving at the matrix A is the quadrature over the basis functions 4>i- Since these integrals involve the fundamental solution, they are nearly always singular and require special numerical treatment. See [1] and [10] for a detailed discussion of this problem. A.4 Discussion T h e boundary element method is a useful and efficient way to solve complex P D E problems by reducing them to boundary integrals. T h e B E M approach is partic-ularly well suited to problems in which the boundary is the primary location of interest (ie. contact problems) and when the domain Q is large compared with the boundary. In other situations, it is often easier to use F E M since a great deal of good software already exists. We note that both integral techniques are specific to a particular P D E and must be re-derived for each new differential operator. For the solution of time dependent problems such as the Navier equations, or the wave equation more theory must be constructed. However the basic approach is the same as presented here. Details of the mathematics can be found in [2] and of the numerics in [22]. 98 Appendix B Simulation Framework B . l Introduction In this appendix, we outline the design for an interactive simulation framework. T h e framework is known as the IS Simulator and represents the combined effort of a small research group over an 18 month period. T h e IS Simulator is a software framework in which interactive physical simulations can be constructed and experienced. T h e motivations of this design are platform independence, extensibility and physical accuracy. One of the primary goals of this simulation environment is to allow plug and play experimentation for new algorithms, such as the contact maps discussed previously. Thi s simulation framework differs significantly from previous general purpose simulation engines like Newton [19] and Impulse [44] in that most user interaction does not require the simulator to be stopped. B y constructing the simulator with multiple execution threads, it is straightforward to implement this feature. In this way, user interaction becomes another component of the simulation which in turn 99 permits a much wider class of applications. B.2 Design Goals Before the simulator project could begin, a number of goals were selected to drive the design. These goals provide the design constraints to which all other aspects of the simulator must conform to. T h e goals can be summarized as: \u00E2\u0080\u00A2 extensibility \u00E2\u0080\u00A2 powerful numerical engine/core \u00E2\u0080\u00A2 responsive/real time performance \u00E2\u0080\u00A2 ease of use Extensibility implies the ability to easily construct and add extensions to the basic functionality. M a n y previous simulators provide no explicit method for improving or extending them. This goal is accomplished by choosing an object oriented design paradigm and implementing it with an object paradigm, in our case Java B E A N S [24]. T h e numerics of the simulator are critical for both accuracy of results and user interactivity. Thi s design goal was accommodated by the M E X X numerical in-tegration package. M E X X provides robust integration of non-stiff D A E (differential algebraic equations) systems. T h e two main bottlenecks in an interactive simulation are the numerics and collision detection. Most naive collision detection algorithms are computationally very expensive. A s this is still a very active area of research, we chose to encapsulate this functionality to allow for upgraded algorithms when they become available. 100 Previous simulation environments often suffered from difficult or obfuscated interfaces. In this simulator, the user is presented with a number of mouse driven views of the simulation in progress. These views facilitate interaction, design and control of the virtual space. B.3 Design Overview We begin by defining an overview of the functional levels of simulator framework. Then we discuss the key components and how they interact. T h e simulator system is broken down into six levels: 1. Simulator 2. Experimental Design 3 . Universe Design 4. Model Edit ing 5. Model Design 6. Numerical Editor We briefly outline each of these levels below, followed by a detailed discussion of the objects that comprise them. B.3.1 Simulator T h e top level component of the IS simulator is the simulator itself. T h e simulator provides the user(s) with a real time view of the current simulation state, including audio, haptic and any other output mechanisms. In addition, the interface allows 101 the user to directly manipulate the simulation while it is running, using methods defined by the experiment design module (see below). T h e display should ideally be a stand-alone client which can be run from within an internet browser such as Netscape\u00E2\u0084\u00A2. From a design perspective, the simulator encapsulates the majority of the functions associated with traditional virtual reality applications. T h e simulator itself is simply a framework in which a number of components interact. In the most basic situation this will include only a numerical engine and a display. These two elements interact in a variety of ways, but primarily through a third entity called a state registry which acts as a kind of central clearing house for simulation data. B.3.2 Experiment Design T h e initial physical setup of the simulated reality is done with the tools of the experiment design level. T h i n g s can be added to a scene, and their initial values defined. T h e experiment designer will be presented with palettes of physical compo-nents. These palettes consist of collections of constructed t h i n g s , such as metal spheres, robot arm assemblies, microphones, or cameras. Amalgamation of t h i n g s (i.e., building a robot arm from parts) is not done at this level but rather in the model editing level. It should be possible for the experiment designer to co-exist with the exper-iment that is currently running. In this case, it would be possible for a designer to add or remove t h i n g s from the simulation while it runs, or to manipulate ob-ject properties. Systematic design manipulations (i.e., design macros) are to be supported with a scripting mechanism. 102 B.3.3 Universe Design In any simulation, there are various constants, such as the gravitational constant, that are unlikely to change during the run. T h e universe design level will be used for setting these properties, which we refer to as the laws of physics. Properties set in this level cannot be changed while the simulator is executing. Although there are certainly times and places where such restriction is important, we achieve a similar result by adding a t h i n g property type which is immutable during simulation execution. B.3.4 Model Editing T h e fundamental components of simulations are t h i n g s , which are created and manipulated using a collection of model designers that form the model editing level. A n incomplete list of t h i n g properties that can be edited includes: \u00E2\u0080\u00A2 meta - th ing generator \u00E2\u0080\u00A2 geometry \u00E2\u0080\u00A2 acoustic properties \u00E2\u0080\u00A2 haptic properties \u00E2\u0080\u00A2 contact force models \u00E2\u0080\u00A2 control properties There are two fundamentally different states that each of these editors will operate in. First they will be used to create new properties for t h i n g s (including other t h i n g s ) . In this mode, the user generates a t h i n g which is saved to the 103 appropriate t h i n g palette for later use. T h e second mode of operation is to ma-nipulate the properties of a t h i n g which has been instantiated in a simulation. In this mode, the editors are hot in that they are part of the simulation loop and can change t h i n g s that are being simulated. For example, a user could create a robot arm assembly, using the various editors to make a fully functional robot arm which is then saved to the robot-arm palette for later use. Once this arm has been incorporated into a simulation, the control editor could be instantiated as part of the robot arm and the control algorithm could then be tuned, etc. to determine it's behaviour. In addition to the standard editors, the t h i n g class will specialize to external data sources, and for this purpose utilities such as parameter estimation and live data filtering are also included in this level. B.3.5 Model Design Fundamentally, all the properties of a t h i n g must reduce to a set of equations (i.e. a system of D A E ' s ) and algorithms, to be integrated numerically. T h e model design level takes care of this detailed work. Various interfaces to the equations should be provided, such as a bond graph editor, block diagram (such as Simulink) editor, etc. Digital filters and control algorithms should also be laid out using the editors of this level. Given a rich set of such equations, it is unlikely the average user would have much use for the editors at this level. However, many users have needs that cannot be foreseen or may want to extend the simulator to novel situations, such as fluid dynamics, or novel types of materials, etc. 104 B.3.6 Numerical Editor A s stated in the last section, underneath it all is a (possibly large) set of equations. T h e equations are accessed though the numerical editor. T h e intent is to provide equation classes that may be specialized for particular problems. Various classes of D A E ' s (differential algebraic equations), P D E ' s (partial differential equations) and O D E ' s (ordinary differential equations) will be represented as well as more esoteric objects, such as D D E ' s (delay differential equations). T h e interface to the user at this level will be through the Java language. B.3.7 Summary T h e six layers of the simulator are much like an onion, they form a nested hierarchy. It is possible to use the simulator at any level, but each level requires knowledge of the levels above it. Casual users will most likely restrict their interaction to the simulator and perhaps the experiment design level. More regular users will use all the levels down to the model editors whereas the serious user and researcher will use all levels. B.4 Object Components T h e simulator levels are realized in the following objects. Functional relationships among these objects are discussed in the next section. B.4.1 Display T h e display object encapsulates all the input/output streams for the user. Each display will run as a separate thread, perhaps remotely and render state informa-105 tion received from the state registry. The display uses information from the thing database to generate appropriate output, such as a 3D display, haptic force, etc. Inputs to the display object are set by default for a standard 'user' thing, but displays may be constructed for other things, such as scene components. Outputs from the display include pick events, mouse motion, speech, etc. The mapping of these values into state variables again will have defaults for the user thing, but are mutable. The concept of the current time is local to each display object, though mul-tiple displays may be configured to remain consistent. The display time is separate from the current state time (discussed below). B.4.2 State Registry Underlying any simulation are a set of equations which relate variables to one an-other. While a simulation is running, the variables evolve according to these equa-tions but at a given fixed moment in time, the numeric value of all the variables is collectively known as the state. The purpose of the state registry is to maintain a record of the system state in time. In our simulator, we extend the state concept slightly to include the user state, which could consist of such things as mouse location, velocity, pick location, etc. and also arbitrary values obtained from external sources of state. The state registry is a passive object which accepts state variable updates, and generates state variable lists in response to queries. The state registry also keeps track of the state time, which is defined as the greatest time at which all state variables are known. The registry may be queried to provide all or some defined subset of the 106 current state. In addition, trajectories are available for all state variables, stored to a pre-defined accuracy. T h e registry may also do internal curve fitting or state history compression in an arbitrary manner within the accuracy constraints. Future states are added for all or some of the variables from a select set of classes. Note that future states are not limited to the next time step, and may be arbitrarily ahead of the state time. B.4.3 Relationship Managers Relationship managers are auxiliary state registry entities which store second order state information. One example of this is inter-object distance which is derived from the current state. Relationship managers may generate state information or serve as repositories for other higher order objects. Relationship managers encapsulate algorithms that operate on the state but are not equations or external inputs. A t this point, the primary relationship manager is the collision manager. This object will monitor the state registry along with the geometries defined in the t h i n g database. Using this information, it calculates inter-object minimum translation distances ( M T D ) and generates collision reports to the integration engine. T h e advantage of a collision manager as a separate object is that it encap-sulates collision detection functionality. Collision detection can be divided into two steps, wide mode (determining which objects to check) and narrow mode (checking distances of a list of object pairs). There are very efficient algorithms for both modes [34], [42], [13]. However a clearly optimal approach has yet to be demonstrated. For this reason the current simulator design only defines an interface to the collision manager. 107 B.4.4 Engine T h e engine is a numerical solver which collects together all the equations governing the simulation variables from the t h i n g database and the managers, and generates the next time ordered set of states. T h e simulator will contain a D A E solver by default, which should be sufficient for most users requirements. T h e current prototype simulator uses the D A E solver M E X X [43] to perform numerical integrations. This solver integrates equations which can be formulated as, P M(t,p)v 0 u = T(T,p)v f(t,p,v,\,u)-G(t,pf\ G(t,p)-v + gI(t,p) d(t,p, v, A , u). (B. l ) (B.2) (B.3) (B.4) (B.5) A great many dynamical mechanical systems can be formulated in this way. In the solution of such systems, matrix equations of the form, (B.6) ( M G T \ ( G 0 V ) must be solved at each time step. B y using recursive elimination algorithms [43], [3] the computational work to solve such systems grows only linearly with the number of bodies being simulated. One drawback of the M E X X integrator is it's inability to deal with stiff equations, for example when there is significant friction or collision damping. This is a result of the extrapolation time stepping technique used. In practice, the engine should run far ahead of the users, predicting the states for a significant future time. If the users do nothing, this simply continues, 108 but when users perform some action or there is an unforeseen event, the engine must reset to that event and continue. There is deliberate latitude for tricks to optimize the engine, for example speculative execution using multiple solvers, or partial re-integration. B.4.5 Thing Database T h e thing database is the collection of entities in the simulation. Thi s database is similar to a scene graph, as found in Java3D [54], Openlnventor [65], V R M L [16], etc. In those systems, the database consists of objects and their properties in a hierarchy. In the IS system however, the database will be a collection of things, which contain their own properties. T h e thing database members are constructed using the experiment design editor and the root node is constructed using the universe design editor. B.4.6 Thing A thing is a simulation entity. Examples of things are spheres, robot arms and user interfaces. Things may be hierarchical in that they may contain things as children. A l l things have a set of properties or parameters which are determined by using a property editor, such as the geometry editor. User related properties might consist of focal length, hearing sensitivity and handedness. Things contain methods which the simulation engine uses to advance the state. For example, a robot arm would contain a control law dictating how it behaves and how it would respond to an external force. Things are passive in that they do not directly generate events or calls to other threads. Things may contain clocks or refer to the simulation or user clock 109 however. In addition, t h i n g s should be implemented using the bean component model and may contain arbitrary code segments (even running in their own threads). B.4.7 Clocks Each separate execution thread must contain a clock object to drive it. T h e exe-cutable code in a clock can be the display update, a simulation engine, an editor, etc. Clock threads usually contain some notion of time, such as user time, state time or other. B.4.8 Data Streams D a t a streams are sources or sinks of state modifying information. D a t a streams may have arbitrary filters and may be connected to external sources, even over the internet. They form the socket abstraction within the simulation environment. B.5 Object Interaction Figure B . l shows the operational relationship among the various components. B.6 IS Simulator Implementation T h e implementation of this simulator is currently in the prototype phase. We seek to characterize the key problems with the realization of this design. B y virtue of it's Java implementation, the simulator runs on a number of platforms including IRIX, Linux, H P - U X and Solaris. O n the Sun UltraSparc-2 computer, the simulator runs in sub-realtime for small scenes with a few constraints. 110 IS Simulator Framework Design Relationship Managers distance manager geometry editor sound editor haptic editor model editor force model editor parameter estimator ( clock \ display user input live data thing thing Figure B . l : IS simulator design overview. I l l Figure B .2: Experiment, design window. Currently the user interface to the simulator consists of three basic windows. T h e experiment design window (see fig. B.2 is where the simulation is setup and all initial values are set. T h e running simulation is viewed through the current attached user view. A n example view is shown in fig. B .3 . T h e user view can be manipulated either directly with a pointing device such as a mouse or with the navigation panel (fig. B .4). Finally, while a simulation is running, t h i n g s may be selected and interacted with dynamically or through a property sheet (see fig. B .5). T h e contact map algorithm discussed in this thesis has not yet been imple-mented in the simulator. Further work to consider the cases of tangential force and multiple contact regions, as well as the extension to 3D need to be completed before contact maps can be used. However prototype simulator has demonstrated that a major bottleneck in dynamic simulation is due to collision detection. Although this was implemented as a native code block, it is clear that collision detection must be optimized in any way possible. 112 Applet started. Figure B.3: A running simulation shown from the user view. 113 Figure B.4: T h e scene navigation bar. B.7 IS Simulator Improvements/Future Work A great deal of work still needs to be done on the IS Simulator project before any firm conclusions can be drawn with regard to it's usefulness. Certainly the ability to interact and simulate in real time has been shown. However the numerical core is still a native component of the simulator. There are no existing numerical integrators written in Java which provide the power of the M E X X engine. One critical future improvement is thus to make the simulator pure Java. T h e current prototype presents only limited user interaction. Upgrades to the prototype could include support for haptic (force feedback) rendering, stereo display and 6 degree of freedom pointer devices. T h e last two are supported in the current specification of Java3D [54]. Haptic rendering presents a special challenge because it must run at a high rate (300-1000 Hz) for objects to feel firm [9]. Finally, a critically important feature is the support of multiple views and/or users. A s an example, consider using the IS simulator to teach an introductory physics course. T h e professor could setup an experiment demonstrating some im-portant dynamic principle. T h e simulator would run on a server and students would connect and be permitted to interact with and observe the apparatus. This func-tionality is not discussed in the design explicitly, but it will be an important feature to provide an implementation for. 114 Figure B.5: A n example of a live t h i n g property editor. While the simulation running, the selected t h i n g s properties may be edited. 115 "@en . "Thesis/Dissertation"@en . "1998-11"@en . "10.14288/1.0080034"@en . "eng"@en . "Mathematics"@en . "Vancouver : University of British Columbia Library"@en . "University of British Columbia"@en . "For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use."@en . "Graduate"@en . "Contact mechanics using Green's functions for interactive simulated environments"@en . "Text"@en . "http://hdl.handle.net/2429/8295"@en .