Virtual Mechanics Simulation and Animation of Rigid Body Systems Hartmut Keller, Horst Stolz, Andreas Ziegler, Thomas Bräunl Univ. Stuttgart, IPVR, Computer Vision Group, Breitwiesenstr. 20-22, D-70565 Stuttgart, Germany e-mail: braunl@informatik.uni-stuttgart.de Abstract: The animation system AERO is presented in this paper. It allows the creation of virtual environments from scenes with simple three-dimensional geometrical objects (sphere, cylinder, box, point, plane). Objects can be linked to each other by a variety of methods (rod, spring, damper, joint). Realistic object movements are achieved by a simulation procedure, where forces can be applied to the objects in addition to gravity, friction, and air resistance. Collisions between objects are simulated using either the ?penalty method? or the ?analytical method?, described in the text. Wire frame animations can be generated in real time, depending on scene complexity and computer system performance. AERO also generates scene description sequences, which can be used as input for ray tracing programs to generate photorealistic animations as batch jobs, a task that is also called ?physically based modelling?. AERO has a number of special features, such as stereo graphics output, camera control, and even the mounting of the camera to an object in the scene. Keywords: simulation, animation, physically based modelling, rigid body systems, virtual mechanics, ray tracing. 1. Introduction Computer animation has been becoming increasingly important in recent years. Of particular interest is animation that results from simulation. In the area of mechanics, reality is reduced to relatively simple abstractions, such as the velocity and position of a mass at a single point, in order to reduce the computational complexity. However, if one gives up this simple abstraction and considers real objects with actual dimensions, the problem becomes increasingly more complicated. In the new scenario there is rotational velocity and acceleration and multiple objects can collide at different angles. In addition, bodies touching each other no longer exert point forces, instead they exert forces over entire areas. These complex systems are no longer expressible in closed form, instead they must be simulated iteratively. For example, how do several balls behave when being dropped from different heights onto an uneven floor? When and how do the balls collide (see Figure 1) ? AERO ( ?Animation Editor for Realistic Object motion?) is an animation system for the visualization of such problems. In AERO a virtual world is simulated, in which all defined bodies move according to the operative physical laws. This system allows the design and execution of physically based computer animations. While the motion of objects in other animation systems often has to be specified by the user, AERO at- Dieses Dokument wurde erstellt mit FrameMaker 4.0.2. Virtual Mechanics 2 tempts to achieve realistic animation solely through the computational application of the laws of mechanics. A scenario containing objects and applied forces is produced with the help of an editor for 3D-data, which was developed during the project. After starting the simulation, the movement and interaction of objects can be observed. Particularly interesting is the option to make changes in the scenario at any point during execution and either continue or restart the simulation. The laudable goal of generating photo realistic animation in real time is not achievable on current workstations. Therefore, an implementation was chosen which allows the presentation of the animation sequences at three levels: 1. Interactive mode: The scenes are produced, calculated and presented in an interactive user interface. The viewing of a simulation sequence can be stopped, wound forward and backwards and altered under user control just like with a video mixer. It is possible to react to the events of the simulation through this interactive intervention. For example, if an external force is to be applied at the exact moment of a collision, one must first know the exact moment of the collision. In AERO the sequence could be started and simulated until the collision occurred. At this animation frame, the sequence would be stopped and the new force applied. 2. Precomputed mode: It is possible with particularly complex scenes containing many objects and complicated interdependencies that the computation of the individual sequence takes a relatively long time. Therefore, it is possible to calculate a sequence off-line (in batch mode) and have it displayed at a later time. 3. Rendering (full calculation): Due to time constraints, the previous two modes allow only a very simple graphical wire-frame representation. In order to achieve photo realistic results, enormous computational expense and therefore considerable time is required. In this mode, a sequence with full details (light, shadows, surfaces, etc.) can be calculated off-line. AERO generates source files for an external ray tracing program. The result, basically a sequence of rendered images, can be displayed as animation on the computer screen using special presentation programs (e.g. ImageMagic or MPEG compression) or stored on videodisc or videotape. Figure 1: Example of a simulation [[[[ gravity x y z Virtual Mechanics 3 2. Creating a Virtual Environment 2.1 Scene Editor After starting AERO, the user enters the scene editor (see Figure 2). The editor allows to enter three-dimensional objects, forces and connections. However, since only twodimensional media are available for input (such as keyboard and mouse) and output (screen), the input occurs through four views. These are (clockwise from bottom left): a view from above (XZ view), a view from the front (XY view), a view from the right (YZ view), and a 3D view diagonally from above to the front and the right. Objects can be selected and placed in the different views using the menus provided. Various fundamental bodies may be selected: spheres, cylinders, cuboids 1 , planes and fixed points. Object parameters, such as size and color, can be altered using a special object window (see Figure 3). By selecting different materials, the user can influence physical characteristics, such as density (thus determining the mass and weight of an object) or coefficient of friction and elasticity. Objects can be connected to one another. Two objects are selected for every connection and the type of connection is specified. Available are rigid connections, rods, springs, dampers and joints, where additional parameters may be specified (for example minimum active separation and spring constant for the spring connections). Multiple connections between two bodies are also possible. For example, one could construct a cuboid with multiple springs connected to a plane at different locations or a shock ab- 1. Cuboids are right parallelepipeds. Basically, these are three-dimensional boxes in which the height, width and length can differ from one another, but all surfaces have right angles at their corners. Figure 2: Editor window Virtual Mechanics 4 sorber realized through a spring and a damper. Complex objects can be constructed this way. The editor window in Figure 4 displays a scenario with a number of bodies and connections between them. In addition, forces can be specified. The user can accelerate or induce rotation of bodies in a defined direction using forces. Forces are specified in three different ways: Figure 3: Object settings Figure 4: Editor window with objects Virtual Mechanics 5 - local to a single body (for example, a force constantly applied tangentially to a sphere in order to create a rotational moment) - relative to another body (a force that is always exerted in the direction of another body, allowing following) - absolute with respect to the reference frame (for the unidirectional acceleration of an object) For the application of forces, a certain time interval has to be specified. 2.2 Animation In Figure 5 an impression of a wire frame animation sequence is given. An animation can be started using a scenario entered via the methods described above. At that point, a new window will be opened, in which the animation occurs (see Figure 6). At the right border of the animation window are buttons resembling those of a video recorder (see Figure 7). These have a very similar functionality in the AERO system. Figure 5: Images from animation sequence Figure 6: Animation window Virtual Mechanics 6 The start button is used to begin the viewing of an animation. The other buttons can be used to influence the viewing of the animation accordingly. The position of the (virtual) camera in the animation window can be specified, in order to view the animation from a different location. This camera position can be altered at any point during the animation playback. An uncommon concept in the AERO system are synchronization points (sync points). Only the start state and current state of the animation are maintained during its playback. All playback functions use state transition functions that convert the current state into the following state. Unfortunately, the calculations cannot simply be reversed in order to achieve the corresponding reverse playback. However, one should still be able to reverse the animation in small steps. This is where sync points are useful. The state of the animation at any point can be declared a sync point. This means that the playback control stores this state and the user can jump back to it at any time. Since any number of sync points can be declared, the sequence can be divided up as finely as necessary. If one wants to reverse the animation by a small amount, one simply jumps back to the previous sync point and inches forward image by image until the desired point is reached. The scene can be altered at any point in time. This is similar to a sync point, however, the transition to the next state does not occur via the state transition functions. Instead, the new scene information is adopted at this scene sync point. If one returns to the editor by closing the animation window, the scenario can be altered at that point in the animation. Film cuts in the animation can be generated this way. The generation of output for the ray tracer can be engaged during the animation. In this case, a new file is generated for every image frame, which contains the current state in the form of the scene description language of the applicable ray tracing program. 2.3 Presentation of Objects AERO bodies are presented differently in the different modes (see Figure 8). In the edtor, it is necessary to be able to determine accurately the distances between objects. A perspective contortion would be out of place and therefore a parallel projection is used for presentation. In addition, connections are displayed only as symbols for the purpose of clarity. Figure 7: Animation functions Start: run sequence step single frame go to start of sequence go back one sync-point Stop: halt sequence go to end of sequence go forward one sync-point Virtual Mechanics 7 The presentation in the animation window of AERO should occur as realistically as possible. Therefore, a perspective projection was selected. In addition, the connections have their correct shapes. A spring connection is shown as a spiral spring and a damper is cylinder-shaped. If photo realistic output via the ray tracer is selected, the objects are rendered with a material-dependent surface texture. In this way, a body made of wood will have a wood grain texture. Glass bodies are clear and light refraction is also taken into account. 2.4 Stereo Image Display Images are two-dimensional due to the fact that they are displayed on a surface and the third dimension, depth, is missing. Humans, can perceive things three-dimensionally, viewing images with two eyes from slightly different viewpoints. The brain can deduce distance information from the differences, or more generally the brain senses depth. In order to achieve a three-dimensional effect for an animation sequence, each of the eyes must be provided with the two-dimensional image that it would have seen in reality (images which differ slightly in their viewpoint). This is achieved in 3D films by using two cameras mounted side by side with eye separation on a tripod. The second image can be generated by a rendering tool just like the first, except that the second virtual camera is about 6-7cm to the side of the first and has a slight inward convergence in order to compensate for the parallax of human eyes. One common method of providing each eye with separate images is the so called redgreen 3D image representation (see Figure 9). The observer wears a pair of glasses, whose right lens is green and left lens is red. The image for the right eye is presented in red simultaneous to the image for the left eye, which is green. Due to the color filtering of the lenses and the fact that green and red are complementary colors, both eyes view only the image intended for them. However, only a pseudo-monochrome presentation is possible. This 3D presentation method is implemented in the animation window of AERO as wire frame graphics. Another more convenient way of viewing stereo images are 3D displays. A liquid crystal display (LCD) in front of the CRT monitor alternately polarizes the light coming from the tube in two opposing directions. In order to prevent flicker, these displays have to use twice the regular image refresh rate. Images corresponding to the Editor Animation Ray Tracing Figure 8: Representation of various objects and connections Virtual Mechanics 8 right or left eye are displayed in synchronization with the direction of polarization. If the observer wears a special pair of glasses with accordingly polarized lenses, then each eye sees exactly one image in full color. Due to the alternating polarization of the display, the observer sees alternately one image in the left eye and then the other image in the right eye. If the image alternation is fast enough, the switching will not be noticed and a three-dimensional effect is achieved. AERO allows the generation of such images when used together with a ray tracing rendering system. 3. Simulation Model The purpose of a simulation model is a computable representation of mechanical processes, which approximates reality accurately enough under the desired conditions. On the one hand, the response time in the interactive system should be as short as possible, and on the other hand as many physical effects and laws as possible should be taken into consideration in order to insure realistic movement. In order tocompute object movements over time, so-called motion equations have to be established. These provide a mathematical relationship between the motion variables like position, velocity and acceleration, and the applied forces. Another important point is the modelling of physical characteristics. The restriction to rigid spatial bodies comes from the well founded and tested knowledge of those systems which can be applied. There are also models which allow for elastic, deformable bodies. However, since the bases of these models are finite element methods, they are not workable in a time critical system such as AERO. Therefore, liquids and gases are also not modelled. Only a rudimentary air resistance is provided. In order to prevent objects from passing through one another, either the impulse of a collision is applied or the forces of touching objects are applied. Before this collision processing can be undertaken, collision recognition has to be done. This includes computation of the contact section between the objects and further data required for collision processing. Constraints 1 define the coupling between bodies in the form of connections. These are input into the motion equations as forces just like gravity and friction. Finally, user specific forces can be applied in order to set bodies in motion. 1. These are usually equations depending on the position of the body and having the form . Figure 9: Stereoscopic red/green presentation of animation j x1 ? xn , , ( ) 0 = Virtual Mechanics 9 3.1 Position and Orientation of a Body in the Virtual Environment To start with, the user must define the position and orientation of a body in threedimensional space. A vector 1 from the origin O of the fixed space coordinate system K to the center of gravity S is sufficient to describe the position of a body. Euler parameters (Quaternions) are used to describe the rotation between the K' coordinate system, which is fixed with respect to the body, and the fixed space coordinate system K. K' is derived by rotating K clockwise around the rotational axis by the angle . Using this scheme, the four Euler parameters are calculated according to: One of the nice features of quaternions is the simple summation of rotations through . Matrix A provides a translation from vector in coordinate system K' to vector in coordinate system K: Using this matrix, we can derive the relationships and . The inverse matrix is obtained by simply reversing the rotational axis . This is achieved by replacing with in the matrix above. A detailed description of Euler parameters, Euler angles and Bryant angles can be found in [Wittenburg 77]; the relationship between two rotations is described in [Shoemake 85]. 1. In figures, vectors are printed in bold type face. Vectors with an apostrophe are with respect to the K' coordinate system. Figure 10: Position of a body and rotation between two coordinate systems rOS u ? q q0 q1 q2 q3 , , , ( ) = q0 ? 2--- cos = q1 q2 q3 , , ( ) T u ? 2--- sin = q u? q0 q , ( ) u0 u, ( ) ? q0 u0 ? q ? u? q0 u? u0 q ? q u? + + , ( ) = = r' r A 2 q0 2 q1 2 + ( ) 1 ? 2 q1q2 q0q3 + ( ) 2 q1q3 q0q2 ? ( ) 2 q1q2 q0q3 ? ( ) 2 q0 2 q2 2 + ( ) 1 ? 2 q2q3 q0q1 + ( ) 2 q1q3 q0q2 + ( ) 2 q2q3 q0q1 ? ( ) 2 q0 2 q3 2 + ( ) 1 ? = r' A r ? = r A 1? r' ? = A 1? u q0 q1 q2 q3 , , , ( ) q0 q ? 1 q ? 2 q ? 3 , , , ( ) OS r x y z x? y? z? S O K? K rigid body x y z g x? y? z? u Virtual Mechanics 10 The following formula, using the matrix above, defines the location of a point fixed with respect to a body: In addition, the following holds for the velocity of the point P: Here represents the angular velocity vector, whose direction specifies the rotational axis about the center of gravity. represents the angular velocity about this axis. 3.2 Rigid Body Systems A model using six fundamental bodies was implemented in an effort to provide basic building blocks, with which the user can experiment. These bodies have geometric forms and physical characteristics corresponding to real materials. In order to minimize computational complexity, a system of rigid or non-deformable bodies should be used. These systems are known in mechanics as rigid body systems. The use of a few simple but generalized spatial forms should avoid the necessity of approximating complex objects with polygonal surfaces and thereby reduce the difficulty of collision processing. A body has several different attributes: mass 1 , physical size and material characteristics. Each of the attributes has impact on the simulation. For instance, motion due to the application of force is only possible when the mass is known. Only bodies with physical size can collide. Finally, the material characteristics influence collision, contact and friction between bodies. Not all of the bodies in Figure 12 have all of the attributes. Except for the immovable point and point mass, all of the bodies have physical size and therefore material characteristics as well. However, an infinite plane has no mass, since it is only a two-dimensional object. This means that planes cannot move. Since the immovable point also has mass zero, it provides a fixed point from which other bodies can be influenced via connections. The point mass can move, however it is not influenced by physical size, i.e. from collisions or contact with other objects. 1. Specifically, mass m and inertia tensor J, where the value of the latter for common spatial forms can be found in any technical reference book. Figure 11: Position and velocity of a point fixed with respect to a body rOP rOS rSP + rOS A 1? r'SP ? + = = u r? = uP vS w rSP ? + = w w OS r OP r SP r? x y x? y? z? S O P z OS r vP vS x y z O P S w rSP Virtual Mechanics 11 Compound bodies can be composed from point masses, cuboids, spheres, and cylinders. These allow the construction of complex rigid bodies. The body parts, which are fixed with respect to each other, are simulated with their corresponding material characteristics. In addition to the three attributes described above, there are others which control visual presentation such as color, texture, etc. However, these attributes are irrelevant for the calculation of the state transition. 3.3 Motion Equations By choosing the center of gravity as reference point for a fixed body, the resulting equations for momentum and spin take on their simplest form. Together they represent the motion equations for a body. Where are the forces applied to the body of mass m, is the acceleration of the mass at the center of gravity and is the rotational acceleration. The matrix is the so-called inertia tensor and represents the spatial mass distribution about the center of gravity. It is comparable to the inertial characteristics of the mass m in the momentum equation. For instance, the rotational velocity of pirouetting ice skaters is smaller with their arms extended due to a larger rotational moment than with their arms retracted. Figure 12: Body types c b a R R cuboid sphere cylinder h plane point of mass immovable point composed object x y z m r??O S ? Fi i 1 = n ? = LinearMomentum JS w?S ? MS ri Fi ? ( ) i 1 = n ? + = AngularMomentum F i r??O S u?S aS = = w?S JS Virtual Mechanics 12 Since the unknown is represented by its derivative with respect to time, these equations are known as a system of differential equations. This system can be solved using numerical integration. 3.4 Connections A connection links a point A of one body with a point B of another body. This is achieved by applying a force corresponding to the type of the connection 1 . Springs, dampers and rods, which consist of a combined spring and damper members, were selected as connection types due to their universal applicability, wellknown functionality and simple implementation. By setting the rod length to zero, a ball-and-socket joint can be approximated.The force equations are calculated at each point in time and the results are fed into the motion equations. It should be made clear at this point, that a rod consisting of a spring and damper must have its length altered in order to exert a force. This is a type of shock absorber. However, by setting the appropriate parameters, one can simulate a rod with sufficient accuracy. 1. is the normal vector between points A and B with and is the distance between the two points. The velocity is given by: . Figure 13: Forces exerted against a rigid body Figure 14: Connections between two bodies Fn r1 MS F1 F2 Fi ri 2r rn w S v S rOS pA pB S SA r Ar B body B body A B A O B n connection n rA rB ? = n 1 = x rA rB ? = r? r?S w p ? + = Virtual Mechanics 13 A general model for coupling different bodies can be achieved using constraints. [Witkin 90], [Barzel et al. 88] and [Isaacs 87] give analytical methods for the more general connection types, including the implementation of fixed rods and joints. These methods are commonly used for technical and scientific purposes 1 and are based on analytical mechanics. These methods were not selected for two reasons: - If contact or frictional forces must also be analytically calculated, the resulting systems of equations are virtually unsolvable. - The rigid connections would also impart momentum during collisions, which would make the collision processing even more complex. 3.5 Adding External Forces In order to manipulate bodies over time, constant forces can be applied at a fixed point on a body during specified intervals. The direction of the force can be specified either in the coordinate system of the body, such as , in fixed space coordinates, such as , or by giving another point on a different body, such as . Smooth, accelerated translational and rotational motion can be generated using these three types of forces. For example, rotational motion can be achieved, as shown in Figure 16, by the application of two equal forces directed within the coordinate system of the body. The translational acceleration of the two forces annihilate each other, leaving the resulting moment: . 3.6 Collisions A generalized representation of a collision between two bodies is shown in Figure 17. In addition to the collision normal and the penetration depth , the determination of the normal velocity of the two bodies at the collision point P is important. 1. Due to the fact that they allow exact calculation of constraint forces. For more information, see [Haug 84]. Figure 15: Types of force F'1 F2 r'4 F3 , ( ) r? r? r? F F F?1 2 3 S SA B r? 4 3 1 2 guided local inertial M 2 r' F' ? ? = n ? uN n u1 w1 p1 u2 ? w2 p2 ? ? ? + ( ) ? = Virtual Mechanics 14 If , then P is a point of contact (see next section). In the case of , the two bodies are separating. A collision only occurs if . The collision value represents the material behavior of the bodies involved: : elastic collision -- a ball thrown against the wall returns with the same velocity ( and exchange values). : partially elastic collision -- incidental velocities are distributed according to . : non-elastic collision -- incidental velocities of the two bodies are the same after the collision. Two routines for collision handling have been implemented: Collision Processing via Springs When a collision occurs, a very stiff spring is inserted at the point of contact. The strength of the spring is proportional to the penetration depth. The spring constant depends on the materials of the colliding bodies. The collision elasticity of the materials is also taken into account. The disadvantage of this method are the large computational costs arising from the high values of the spring constant increasing the time granularity of the integral so- Figure 16: Generation of rotational moment using two forces Figure 17: Collision between two bodies S F? -F? r? -r? M S p p1 S S body 2 body 1 2 1 2 v v 1 2 w w1 2 P n body 1 body 2 collision plane collision normal P uN 0 = uN 0 < uN 0 > e e 1 = uN1 uN2 0 e 1 < < uN1 uN2 ? ( ) e u'N1 u'N2 ? ( ) = e 0 = c F ? cn for uN 0 ? (approaching) e ? cn ? for uN 0 < (separating) ???? ? = c Virtual Mechanics 15 lution process. If is selected too small, the bodies will penetrate each other depending on their mass and velocity and may even pass through each other. Analytical Collision Handling The analytical method determines the new velocities after the collision from the current translational and angular velocities by applying momentum and spin conservation principles (no energy is lost). In order to ensure that the collision impulse , which is exchanged between the bodies during the collision, is exerted in the direction of the collision, the following restrictions are used in addition: This system of linear equations with 15 unknowns 1 , and can be solved with standard methods such as the Gauss algorithm. Altogether the analytical method is faster than the method using springs because it only needs to be carried out once at the point of collision and not over a long period of time. The path of a sphere falling to the ground from a height of 0.05m can be seen in Figure 19. Based on [MacMillan 60], [Moore et al. 88] indicates that one can take both friction as well as the propagation of collision impulses through connected bodies into account. [Wang et al. 92] describes the case of friction for two-dimensional applications and [Wittenburg 77] dedicates a chapter to collision processing for systems containing compound bodies. However, the critical problem of handling multiple simultaneous collisions at one body remains therein unsolved. 1. This can easily be reduced to 13 unknowns since runs parallel to . Figure 18: Collision processing via springs c body 1 body 2 F -F g P u'1 w'1 u'2 w'2 , , , ( ) u1 w1 u2 w2 , , , ( ) m1u'1 m1u1 R + = m2u'2 m2u2 R ? = J1w'1 J1w1 p1 R ? + = J2w'2 J2w2 p2 R ? ? = R R e1 ? 0 = R e2 ? 0 = n u'2 w'2 p2 u'1 w'1 p1 ? + ( ) ? ? + ( ) ? e n u'1 w'1 p1 u'2 w'2 p2 ? + ( ) ? ? + ( ) ? ? = u'1 u'2 w'1 w'2 , , , R R n Virtual Mechanics 16 Consider a cuboid, such as the one in Figure 20, that experiences collisions at two different places simultaneously. Simplified, the following results from the collisions: This means that the cuboid would not begin to rotate, even though the collisions are asymmetrical. As long as there is no easily applicable physical model for multiple, simultaneous collisions, one has to handle simultaneous collisions serially. The AERO system carries out collision processing on points with until all points have separated from each other with . However, a maximum of steps are carried out, where is the number of collisions occurring. Figure 19: Path of a sphere falling to the ground with = 0.81. Figure 20: Cuboid with two asymmetrical points of collision -0.01 0 0.01 0.02 0.03 0.04 0.05 0.06 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 "analytical solution" e v v R R a b a b m mu' mu Ra Rb + + = J w' Jw Ra ? Rb ? = u' w'a + 0 = u' w'b + 0 = ?? ??ü since a b ? w' fi 0 = uN 0 > uN 0 ? 2k k Virtual Mechanics 17 This form of collision propagation has the disadvantage, that if a cuboid falls flat onto a level surface, it begins to rotate on an axis parallel to the surface as shown in Figure 22. For this reason, the user can select for collision handling either the analytical method or the spring method (which does not have the problem shown above), whichever is best suited for the intended simulation. 3.7 Handling of Physical Contact Theoretically, one could apply the collision handling methods discussed in the previous chapter to points of contact having as well. In fact, this is exactly what happens when the user selects the spring method for collision handling. In contrast, the analytical method of collision handling is too computationally complex to use with every state transition. Instead, the ?penalty method? is used, which has a better transient behavior than the spring method. Every point having less than a specifiable maximum contact velocity will have a force applied opposing the direction of penetration. This force is proportional to the penetration depth and proportional to the velocity along the contact normal . The elasticity factor and the damping factor are available as parameters for every type of material. Extended transient periods as well as strong oscillation upon contact can be avoided during the contact process through a suitable selection of . Figure 23 shows the settling of an iron sphere in meters. The horizontal time axis is in units of seconds. The analytical handling of contact forces, as described in [Baraff 91], was not implemented. However, it certainly represents a worthwhile goal, because it avoids the problem of different penetration depths of bodies as shown in Figure 24. In addition, it would allow the consideration of constraint forces. A subroutine for solving LCP 1 problems using Lemke's Complementary Pivot Algorithm was implemented as well and may be included in an extended version. 1. Linear Complementarity Problem: find the vectors and for , where the matrix and the vector are given. For more details, see [Murty 88]. Figure 21: Propagation of collision impulse with three spheres Figure 22: Cuboid falls flat onto plane v v v 1st collision 2nd collision ? w w?? R R 1 2 v? v?? v uN 0 ? uN uNmax ? uN c d d w z w Mz ? q w ? 0 ? z ? 0 ? wizi ? 0 = = i " M q Virtual Mechanics 18 3.8 Friction When two touching bodies move relative to one another, friction creates a force opposing the direction of movement. This type of friction is known as dynamic friction. The force created by the dynamic friction can easily be calculated from: Figure 23: Comparison of methods for settling a sphere in to a plane Figure 24: Contact handling using the penalty method Figure 25: Dynamic friction induces rotation in a sliding sphere -0.0009 -0.0008 -0.0007 -0.0006 -0.0005 -0.0004 -0.0003 -0.0002 -0.0001 0 0 0.02 0.04 0.06 0.08 0.1 0.12 "penality-method" "spring-impact-method" light sphere heavy cuboid plane F F v N R S G FR0 FR0 uT uT --------- ? ?0 FN ? ? = Virtual Mechanics 19 is the force of contact which was calculated in the previous section. The coefficient of dynamic friction is a value in the range and depends on the pair of materials in contact. As has been seen before, implementing dynamic friction is no problem. Static friction, on the other hand, must be approximated. The force of static friction usually applies only when the tangential velocity at the point of contact and holds the body at a standstill. In addition, always lies within the plane of contact and is never larger than the force of contact , or more exactly: . In AERO, the restriction is lifted and all points of contact are handled with . The dynamic friction formula is used for points of contact having a tangential velocity above the maximum . Instead of having an unknown direction of application, the force of static friction is applied opposing the direction of movement within the plane of contact. This results in the following approximation: It should be apparent, that this static friction will not bring a body to an absolute standstill. One can only select the smallest possible parameter values so that no visual affects are apparent in the animation. The problem here is that this would require the integration method to use very fine steps. Actually, the contact forces and frictional forces should be placed in a large system of equations in order to achieve an analytical solution, which would provide the correct solution. As shown in [Baraff 91], such mixed systems of equations are difficult to solve with current mathematical methods. 3.9 Gravity and Air Resistance In order to prevent bodies from behaving as if they were in zero gravity in outer space, a constant force of gravity is applied everywhere in the negative direction. Naturally, a value for differing from 9,81 can be selected in order to move the simulation to the surface of another planet. Every body also has an extra attribute ?gravitation? that allows one to turn off the gravity for just that particular object. Another option is allowing the force of gravity to apply somewhat randomly with a specified percentage. This has the advantage, that special cases such as stacking spheres on top of one another do end with the expected collapse. On the other hand, it should be pointed out that replacing constant gravity with a random function having an average value also has unwanted side effects. This will lead to uncontrolled energy changes in free falling bodies. The energy of a body is calculated as: The random directly influences . For a pendulum, which changes potential energy into kinetic energy, the energy changes would lead to changes in the height and period of the swing. Another side effect is that stationary objects sitting on the ground are more easily brought to sliding. In addition to the force of gravity, a force of air resistance can be applied opposing the direction of travel. This force increases with the square of the velocity and includes a FN ?0 0?1 [ ] FR uT 0 = FR FN FR ? FN ? ? uT 0 = uT uTmax ? uTmax FR uT uTmax -------------- ? ? FN ? ? = g y g m s2 ? g E Ekin Epot + m 2---- u2 1 2--- J w2 + = = g u? a = Virtual Mechanics 20 body specific air resistance factor . The attack surface is calculated as a constant from the dimensions of the body. The formula used in AERO does not reflect true flow behaviour, but serves quite well as a simple way for generating a resisting force: 3.10 Collision Detection In order to engage the collision handling process, collision detection must recognize the penetration of one body within another and calculate the data required for the collision or contact handling methods. Collisions are recognized by checking each body against every other body for overlapping regions. Those pairs of bodies that are in contact are collected in a list. The assumption is made, that only a few of the points of contact need to be taken into consideration. For instance, when two cuboids are in contact, only the corners of the overlapping region are calculated. These points are then used for the application of contact and frictional forces as well as for the collision impulse. The parameters relevant for a collision are the collision point and the vectors from the bodies' centers to that point: and . Also important are the collision normal and the penetration depth . It is necessary, that a suitably large countering force be applied to a point so that the penetration depth is reduced. Exactly this behavior is required of the contact methods. In spite of the fact that geometrically simple bodies were selected, routines are requires to handle the bodies (cuboid, sphere, cylinder and plane). Unfortunately, even with these simple fundamental bodies, the collision tests proved to be computationally expensive. Therefore, due to run time considerations, the cylinder collision routines have not been imlemented yet. Cylinders are approximated by cuboids and handled using the cuboid collision routines. Figure 26: Collision detection cw A FAirResistance c ? w A u u ? ? ? = rP2 body 1 body 2 P r 0 P1 n collision plane e e 1 2 g Pr normal a rP rP1 rP2 n ? P ? n 2-- n 1 + ( ) ? 10 = n 4 = Virtual Mechanics 21 4. Computational Methods 4.1 Integration of the Motion Equations Since the motion equations form a system of differential equations, solving them via numerical integration is appropriate. Most integration methods are only applicable for systems of the form: The process will calculate the new state from the old state . However, the motion equations represent a system of differential equations of the second order. In order to construct a system of the first order, the following assumptions are made: Here, the main reason for using the Euler parameters to describe the orientation of a body is revealed. In contrast to Euler angles, which use three serially executed axis rotations, Euler parameters are more easily integrable. The following equation holds: The state vector contains for each body the position , the orientation , the linear velocity and the angular velocity . The right side of the equation above corresponds to the problem dependent function vector . The integration method is responsible for accumulating or summing the accelerations and , the velocities and , and the position and orientation . 4.2 Integration Method The selection of an integration method was suggested by [Heinzel 92] and involves an imbedded Runge-Kutta procedure 1 with error control and interpolation of 1. This is known as RK5(4)7FM in [Dormand et al. 80] and [Heinzel 92]. x? f t x , ( ) x t0 ( ) ? x0 = = x t h + ( ) x t() r?? 1 m---- F i i? = w? J 1? M r i F i ? ( ) i? + ? ? ? ö ? = r? u = q? g q w , ( ) = q? 0 q? 1 q? 2 q? 3 g q w , ( ) 1 2--- 0 w ? 1 w ? 2 w ? 3 w1 0 w3 w ? 2 w2 w ? 3 0 w1 w3 w2 w ? 1 0 q0 q1 q2 q 3 ? = = x r q u w f t x ,( ) r?? w? r? w r q Virtual Mechanics 22 intermediate states. The method of the order calculates from the intermediate values the new state as The step size is selected through the error control according to the desired accuracy. By selecting another state value of a lower order using one can get an approximation of the error between the values for and from . If the accuracy is specified, then whenever holds, the state must be discarded and recomputed with a new step size , to get more accurate results. A suitable formula for the new step size is: Here, is the order of the state equation for , which is one less than that of . In an actual implementation, it is useful to constrain the step size to a specific interval in order to avoid unacceptably long response times. By using interpolation based on the previously generated intermediate values , the system can efficiently generate intermediate states within the interval . The state changes within this interval are approximated by using a polynomial of degree . This polynomial's term coefficients are Computing an intermediate state occurs without evaluating the function by using: The constants and for the formulas above can be found in [Heinzel 92] or in the original texts [Dormand et al. 80], [Dormand et al. 81] and [Dormand et al. 86]. s 7 = ko f t x t() , ( ) = k i f t c i h ? + x t() h a i j k j j 0 = i 1 ? ? ? + , ( ) = i 1 ? s 1 ? , , = x t h + ( ) x h b i k i i 0 = s 1 ? ? ? + = h x' x' t h + ( ) x h b' i k i i 0 = s 1 ? ? ? + = t t h + d x x' ? = e d e > x t h + ( ) hnew hnew 0,9h e d-- ? ? ? ö 1 p 1 + ------------ = p 4 = x ' q p 1 + 5 = = k0 ? k s 1 ? , , x t sh + ( ) t sh + t t h + , [ ] ? d 5 = a0 x t() = a j h b j 1 ? i,ki i 0 = s 1 ? ? = j 1 ? d , , = f t x ,( ) x t sh + ( ) a j s j j 0 = d ? = a i j ci bi b'i , , , bi j, Virtual Mechanics 23 4.3 General State Transition Operation The computational structure in Figure 27 shows the interaction of the individual components of the state transition calculations used in the AERO system. From a state containing the location and movement of all bodies as well as the states of the connections and user specifics forces in effect, a new state at the point is calcu- lated. The integration method has the task of summing or accumulating the results of the motion equations from to . However, it became apparent, that the interval should first be divided into smaller pieces known as in order to allow a finer grained collision detection within the interval. In addition, when using the analytical collision handling method, the integration method must be restarted, since the newly calculated velocities are seen as start values for the systems of differential equations. The motion equations are only nested within the integration method in order to achieve the correct processing order. The integration method is implemented fully independent of the problem at hand. The selected modelling of the forces (gravity, contact, friction and air resistance) allows direct calculation from the state variable . Since they only use the position and velocity, they are directly inserted into the linear and angular momentum equations. However, in order to take the contact forces into account, a collision detection must have occurred previously. This is handled by recording the collisions as soon as they are recognized and storing them to be handled at a later point. One should note that since the step size is variable, the integration method occurs asynchronously to the requested states at . When the integration method takes larger steps, the intermediate states are approximated via interpolation. If it takes smaller steps, intermediate states must be inserted. Figure 27: Schematic of the state transition calculations t D t + N FOR t' = t TO t + t STEP t Kol D D D equations of motion Y state t + tD impact impact at t? = t , restart of integration needed impact t within [t?, t?+ t ] D Col determine earliest collisiontime t interpolation if last step t+h > t? + t D Col supply k - k for the computation of f(t+h,x) 0 s-1 Runge-Kutta-method with variable stepsize h asynchronous integration from t? up to t?+ t by Col use penality-method to determinate the contact forces analytical impact computation state t, timestep t D impact collision detection function evaluation of f(t,x) at intermediate points t t Dt + D t Co l x h t' t 'Col + Virtual Mechanics 24 4.4 Simulation Requirements In order to carry out a simulation true to the model, four conditions must be met. Due to their mutual interaction, these conditions cannot be viewed as being independent of each other. Finally, the resulting animation sequence should be analyzed for unexplained motion changes or large position changes, in order to fine tune the simulation parameters. Numerical Stability The integration method works stably if the accumulation of a derived function sufficiently approximates the function itself: In order to achieve this, the step size must be selected appropriately for the problem function . The error control selects a suitable value for independently, to which the accuracy of the intermediate steps is applied. Since only varies within a specified interval, if the lower bounds of the interval, corresponding to the minimal step size, is reached, this indicates that the accuracy specified can no longer be met. In this case, the sum may diverge considerably from the nominal value . In extreme cases, this will lead to full separation, which will be noticeable in the simulation as sudden position changes of the affected object. This could deteriorate to the point that an object totally disappears when its position coordinates diverge rapidly toward infinity. Sufficient Collision Detection A timely recognition of collisions is crucial for processing collisions and physical contact. The collision step size represents the maximum time interval between two collision detection invocations. It should be selected based on the spatial dimensions of the simulated bodies and the maximum velocity. Selecting too large a value for has many negative consequences: - Bodies pass through each other because the time interval of their overlap is smaller than . - The penetration of bodies is recognized too late, which leads to very large penetration depths . This leads to enormous contact forces, which cause the objects to jump away from each other. Furthermore, this breaks up the smooth application of the contact forces, which in turn drives the step size of the integration method down to its minimum value. Moderate Collision Classification Using the analytical collision handling method, the system must first classify the intersection of two bodies as a collision or contact. This is accomplished via a comparison with the maximum contact velocity . This velocity cannot be considered to be independent of the collision step size. f x() I x( ) ? f? x' ( ) x' x0 = x ? = h f x() h e h I x() f x() D t Col D t Col DtCol ? uNmax Virtual Mechanics 25 A body which rests upon an unmoving base should not be supported by collision forces, it should be supported by contact forces. The difference in velocity between two bodies after a collision is approximately . Practical values for should be 10 to 100 times larger, in order to cover higher velocities as well. These could occur when forces in addition to gravity are applied to a body. Having too small a value leads to collisions during the settling process. In addition to extending the settling period, this also consumes more computing time due to the unused collision inertia calculations and the automatic restart of the integration method. Appropriate Parameters for Body and Material The modelling of contact forces together with the spring collision method and the connection types spring, damper, and rod result in a feedback control loop which must be kept stable. The principal parameters of this feedback control loop are the elasticity and damping factors of the materials and connections as well as the gravitational force of the bodies. Since the mass of a body comes from its physical size, a prudent selection of the material type is important. Radically different body masses should be avoided. Elasticity and damping factors should protect the settling process from the elastic behavior of the connections. This property is comparable to the effect of shockabsorbers in an automobile. The goal here is to have as short a settling period as possible and to avoid over reaction upon contact. Previously, Figure 23 showed the settling period of a sphere on a plane using the penalty method and the spring method with the same elasticity factor. Since the spring method only smooths collisions based on the collision value , slowly fading oscillations with high amplitudes can occur. For this reason, simulation of collision problems using high velocities should not use the spring method. Only by selecting larger material factors can the maximum penetration depth be contained. However, such ?rigid? problems have large computational costs. Problems with radically different time factors in the differential equations are known as ?rigid? differential equations. The are only solvable by using very fine time intervals. In contrast, the penalty method takes the difference in velocity at the collision point into account and thereby more quickly achieves a standstill without large oscillations. The default values for the damping and spring constants were generated by analyzing the settling behavior of spheres from 1 to 600kg. A universal selection of the factors, which provides optimal results for every thinkable case is not possible. A simulation should therefore be built from bodies and factors that are suitable for each other. Bodies in the simulation which should not move can be tagged as massless. This saves processing time since these bodies are not added to the motion equations and collision detection between massless bodies is not required. 5. Sample Animations 5.1 Pendulum Five rubber balls are connected to fixed points via rods as shown in Figure 28. One of the outside balls is brought to the horizontal position so that it falls due to its weight in a circular path toward the other still balls. In this way, a system of five mathemati- g D t Col ? uNmax uNmax e Virtual Mechanics 26 cal pendulums is constructed. The energy of the outer pendulums is propagated through collision impulses to the other pendulums. In order to beautify the simulation, a supporting stand can be built. The cylinders comprising the stand must be defined as immobile. In order to avoid changes in the length of the rods, one can increase the elasticity and damping constants by a factor of ten. After all bodies and connections have been entered, the animation can be started. Normally the user first selects a suitable frame frequency. For example, Figure 29 shows an animation sequence that was produced with 15 frames per second. Every image frame is stored, so after recording, the whole animation sequence can be played forward and backwards without a large computational overhead. In this way, one can find the optimal camera position for viewing the sequence. Once the animation is setup perfectly, the ray tracer code for each image frame is generated and the ray-traced images can be computed in batch mode. 5.2 Can Toss The next example shows a can toss simulation well known from county fairs: six cans are set up in a pyramid shape and all of them must be knocked down with a single ball (see Figure 30). It is interesting to see how realistically the cylindrical cans fly in Figure 28: Step-by-step construction of a five sphere pendulum Figure 29: Animation sequence of the five sphere pendulum Virtual Mechanics 27 all directions after being struck by the ball. Or expressed more mathematically: how the impulse of the ball is propagated to the previously still cylinders. 5.3 Vehicle Simulation This is an example of a longer sequence. Due to space restrictions, only a few images are shown here (approximately every 30th image). The original version of this sequence is approximately 500 image frames long and takes up 440MB of disk space in uncompressed format. Each image is pixels large, and when shown with 25 frames a second, the sequence has a duration of about 20 seconds. This high playback speed can not be directly achieved on a workstation class computer. However, using a video recorder with single frame write capability or a laser disc system, one can write each image frame individually. Thereafter, the image can be played back at full speed. In the future, image compression algorithms such as MPEG and faster workstations could allow software playback at full speed and at full size. The sequence (see Figure 31) shows a three wheeled vehicle with shock absorbers and spherical wheels. The vehicle starts out standing on small elevated surface, which is implemented as a cuboid with a simple pendulum in the background. It is accelerated by inducing torque on the front wheel, so it rolls off the stairs (note the shock absorbers). The vehicle passes under a stylized tree, encounters a small hill, causing it to turn to the right onto an ice surface, where it receives an angular momentum and spins to a stop. Figure 30: Can toss 640 480 ? Virtual Mechanics 28 6. Acknowledgements In the beginning, the idea behind the AERO system was to create a simulation of the real world according to elementary physics. Particularly the articles from David Baraff provide convincing evidence of the feasibility of such a project. A collection of various articles came together quickly, in which each article very accurately takes individual physical effects into consideration. However, the sum of the articles did not provide a functioning model. It became apparent, that the well known formulas of mechanics are mainly directed toward special cases and introductory lectures on the subject are normally simplified to demonstrate the special cases. Universal analytical handling, in contrast, is only possible through the principles of mechanics or through approximations such as the penalty method. We would like to thank Brian Blevins for translating the original German manuscript and Prof. Levi for his support. We also acknowledge the work of Drew Wells and colleagues for their public domain POV (?persistance of vision?) ray tracing package, available via ?anonymous ftp? from alfred.ccs.carleton.ca (134.117.1.1), as well as the work from Lawrence Rowe, Kevin Gong, Ketan Patel, Brian Smith, and Dan Wallach from the University of California at Berkeley for their public domain MPEG encoding and decoding package, available via ?anonymous ftp? from toe.cs.berkeley.edu (128.32.149.117), directory pub/multimedia/mpeg. We used their systems for rendering frames and composing animation sequences. The AERO system has been implemented in C/Unix with X-Windows and tested on Sun workstations and IBM-PC/linux systems by Andreas Ziegler (3D-scene editor), Hartmut Keller (3D-graphics, spatial representation, scene generation for ray tracing) and Horst Stolz (simulation computation) under the direction of Thomas Bräunl. The Figure 31: Vehicle on ice Virtual Mechanics 29 AERO system is also available as public domain software via ?anomymous ftp? and may be copied from our server: ftp.informatik.uni-stuttgart.de (currently 129.69.211.2) in subdirectory: pub/AERO . 7. References [Barzel et al. 88] Ronen Barzel, Alan H. Barr: A Modeling System Based on Dynamic Constraints; pp. 179- 187 in Computer Graphics (Proceedings of SIGGRAPH'88), ACM Press, Atlanta, vol. 22, no. 4, August 1988 [Baraff 91] David Baraff: Coping with Friction for Non-penetrating Rigid Body Simulation; pp. 31-40, in Computer Graphics (Proceedings of SIGGRAPH `91), ACM Press, Las Vegas, vol. 25, no. 4, July 1991 [Dormand et al. 80] J. R. Dormand, P. J. Prince: A Family of Imbedded Runge-Kutta Formulae; pp. 19-27 in Journal of Computational and Applied Mathematics, Elsevier Science Publishers B.V., North-Holland, vol. 6, no. 1, 1980 [Dormand et al. 81] J. R. Dormand, P. J. Prince: High Order Embedded Runge-Kutta Formulae; pp. 203-211 in Journal of Computational and Applied Mathematics, Elsevier Science Publishers B.V., North- Holland, vol. 7, no. 1, 1981 [Dormand et al. 86] J. R. Dormand, P. J. Prince: A Reconsideration of some Embedded Runge-Kutta Formulae; pp. 203-211 in Journal of Computational and Applied Mathematics, Elsevier Science Publishers B.V., North-Holland, vol. 15, 1986 [Haug 84] E. J. Haug (Ed.): Computer Aided Analysis and Optimization of Mechanical System Dynamics; NATO ASI Series, Series F - Computer and System Science, vol. 9, Springer Verlag, Berlin, Heidelberg, 1984 [Heinzel 92] Gerhard Heinzel: Beliebig genau -- Moderne Runga-Kutta-Verfahren zur Lösung von Differentialgleichungen; pp. 172-185 in c't, Verlag Heinz Heise GmbH & Co KG, Hannover, no. 8, 1992 [Isaacs 87] Paul M. Isaacs, Michaell F. Cohen: Controlling Dynamic Simulation with Kinematic Constraints, Behavior Functions and Inverse Dynamics; pp.215-224 in Computer Graphics (Proceedings of SIG- GRAPH'87), ACM Press, Anaheim vol. 21, no. 4, July 1987 [MacMillan 60] William D. MacMillan: Dynamics of Rigid Bodies; Dover Publications, New York, 1960 [Moore et al. 88] Matthew Moore, Jane Wilhelms: Collision Detection and Response for Computer Animation; pp. 289-298 in Computer Graphics (Proceedings of SIGGRAPH `88), ACM Press, Atlanta, vol. 22, no. 4, August 1988 [Murty 88] K. G. Murty: Linear Complementary, Linear and Nonlinear Programming; Heldermann-Verlag, Berlin, 1988 [Shoemake 85] Ken Shoemake: Animating Rotation with Quaternion Curves; pp. 245-254 in Proceedings of SIGGRAPH'85, ACM Press, San Francisco, July 1985 [Wang et al. 92] Yu Wang, Matther T. Mason: Two-Dimensional Rigid-Body Collision with Friction; pp. 635- 642 in Journal of Applied Mechanics, vol. 59, September 1992 [Witkin 90] Andrew Witkin, Micheal Gleicher, William Welch: Interactive Dynamics; pp. 11-21 in Computer Graphics, ACM Press, 1990 [Wittenburg 77] Jens Wittenburg: Dynamics of Systems of Rigid Bodies; Teubner, Stuttgart, 1977