Master Thesis
Validation toolbox for a Physics Engine
Author:
Emma Sundling
emsu0033@student.umu.se
Advisers:
Kenneth Bodin kenneth.bodin@algoryx.se Claude Lacoursi`ere claude.lacoursiere@algoryx.se
June 2016
Physics Department
Examiner Eddie Wadbro eddiew@cs.umu.se
W. Clement Stone
Master thesis 30hp
Validation toolbox for a Physics Engine
Physics engines become more and more common due to the rapid development and increasing demand of simulations. With this comes a need of testing the engine, a way to measure its performance, not only its speed but also its accuracy and stability. The purpose of this thesis has been to create a set of benchmark tests. They aim to check the physical aspects, especially mechanics, of the engine. A strategy and export functions for the test results in order to automate the testing have also been developed.
The resulting tests became a beam on piles which analyses constraint stability, an overdetermined system consisting of a static door on multiple hinges, a falling object investigating the accuracy of the integrator, a box on an inclined plane for testing the friction model, a single pendulum as well as a multibody pendulum checking constraint accuracy and energy conservation, the Earth orbiting around the Sun which tests the stability of the integrator and finally a cantilever beam that is a static test of a real scenario. After the tests are performed the results are presented on an HTML-page. A prototype of a Web application is also established as well as a set of scalar tests that can be performed continuously, in order to follow trends or compare the engine’s performance from time to time.
This thesis was initialized by Algoryx Simulation AB which also provided the engine, AgX Dynamics, with the numerical method called Spook. It mainly performed well on all tests. In order to build a fully general toolbox more tests need to be added such as material interactions, scalable test with thousands of bodies, torque tests as well as more complex scenarios, for example a scissor lift and robots. The work can also be extended with more developed export functions, both to the Web and to documents. Hopefully this thesis can be seen as a complement to the earlier efforts done in creating a general set of benchmarks and automation framework for continuous integration and testing.
Examensarbete 30hp
Valideringsverktyg f¨or en fysikmotor
Fysikmotorer blir mer och mer vanliga p˚a grund av den snabba utvecklingen och efterfr˚agan p˚a simuleringar. I och med detta ¨okar ocks˚a behovet av att testa motorerna och ett s¨att att m¨ata prestandan, inte bara snabbheten utan ocks˚a noggrannheten och sta- biliteten. Syftet med detta examensarbete har varit att skapa ett set av prestandatester.
De syftar till att testa de fysikaliska aspekterna av fysikmotorn, s¨arskilt inom mekanik.
En strategi och exportfunktioner f¨or testresultaten f¨or att automatisera testningen har ocks˚a utvecklats.
De resulterande testerna blev en balk p˚a p˚alar som analyserar stabiliteten hos villkoren, ett ¨overbest¨amt system best˚aende av en statisk d¨orr p˚a flera g˚angj¨arn, ett fallande objekt som granskar precisionen hos integratorn, en l˚ada p˚a ett lutande plan som testar friktionsmodellen, en enkel pendel samt en flerkropppspendel som kontrollerar villkorsprecisionen och energikonservering, jordens bana runt solen som testar integratorns stabilitet och slutligen en utskjutande balk som ¨ar ett statiskt test av ett verkligt fall.
N¨ar testerna ¨ar genomf¨orda presenteras resultaten p˚a en HTML-sida. En prototyp av en webb-applikation har ocks˚a utvecklats samt ett set med skal¨ara tester som kan utf¨oras kontinuerligt f¨or att f¨olja upp trender och j¨amf¨ora motorns prestanda ¨over tid.
Det h¨ar examensarbetet initierades av Algoryx Simulation AB som ¨aven tillhandah˚allit fysikmotorn, AgX Dynamics, med den numeriska metoden Spook. Motorn presterade
¨
overlag bra p˚a testerna. F¨or att bygga en allm¨an verktygsl˚ada beh¨ovs fler tester s˚a som interaktion mellan material, skalbara tester med tusentals kroppar samt mer komplexa simuleringar, t.ex. en saxlyft och robotar. Arbetet kan ocks˚a ut¨okas med mer utvecklade exportfunktioner, b˚ade mot webben och som dokument. F¨orhoppningsvis kan detta ses som ett komplement till de tidigare anstr¨agningar som gjorts f¨or att skapa ett generellt set av prestandatester och ett automatiskt ramverk f¨or kontinuerlig testning.
Firstly I want to thank my advisers Kenneth Bodin and Claude Lacoursi`ere for all the help and support throughout this work. The road has certainly not been straight and I am very grateful for all your patience, knowledge and encouragement. Besides my advisers, I would like to sincerely thank Eddie Wadbro who always been ready with answers and good advices, no matter when or where. I also want to express my gratitude for all the help and support provided by the people at Algoryx. They have not only contributed with knowledge and guidance but likewise in practical ways. Finally I would like to thank my family and partner for their support both during this work and during life in general.
iv
Abstract ii
Sammanfattning iii
Acknowledgements iv
Contents v
List of Figures vii
Physical Constants and Notation ix
1 Introduction 1
1.1 Background . . . 2
1.2 Related Work . . . 2
1.2.1 Robotics . . . 3
1.2.2 Offshore . . . 4
1.2.3 Web and Games . . . 4
1.2.4 Conclusion . . . 5
1.3 Purpose and Goals . . . 6
1.4 Delimitations . . . 6
2 Preliminaries 7 2.1 Physics engine . . . 8
2.1.1 Collision detection . . . 8
2.1.2 Simulation . . . 10
2.1.2.1 SPOOK - Pendulum . . . 10
2.1.2.2 SPOOK - Harmonic Oscillator . . . 14
2.1.3 Key features . . . 20
2.2 Stable simulations . . . 21
3 Test cases 24 3.1 Work flow . . . 24
3.1.1 Pre-study . . . 24
3.1.2 Implementation . . . 25
3.2 Test cases . . . 25
v
3.2.1 Beam — Constraint stability . . . 25
3.2.1.1 Results . . . 26
3.2.1.2 Discussion . . . 29
3.2.2 Door — Overdetermined system . . . 29
3.2.2.1 Results . . . 31
3.2.2.2 Discussion . . . 33
3.2.3 Falling Object — Accuracy . . . 33
3.2.3.1 Results . . . 34
3.2.3.2 Discussion . . . 36
3.2.4 Friction . . . 37
3.2.4.1 Results . . . 39
3.2.4.2 Discussion . . . 41
3.2.5 Pendulum — Constraint accuracy . . . 42
3.2.5.1 Results . . . 43
3.2.5.2 Discussion . . . 46
3.2.6 N-pendulum . . . 46
3.2.6.1 Results . . . 47
3.2.6.2 Discussion . . . 50
3.2.7 Planetarium — Integrator stability . . . 51
3.2.7.1 Results . . . 52
3.2.7.2 Discussion . . . 54
3.2.8 Cantilever beam . . . 55
3.2.8.1 Results . . . 56
3.2.8.2 Discussion . . . 57
4 Additional work 58 4.1 Result presentation . . . 58
4.2 Continuous testing . . . 59
4.3 Interactive Web based Simulation . . . 60
5 Discussion 62 5.1 General . . . 62
5.2 Continuous testing . . . 63
5.3 Interactive Web based simulation . . . 63
5.4 Conclusion . . . 63
A Strategy and Work flow 65
Bibliography 66
2.1 Different bounding volumes . . . 9
2.2 The theoretical set-up of the pendulum. . . 11
2.3 The steps of shake . . . 12
2.4 Spring - position vs time, both simulation and manual stepping . . . 15
2.5 Spring - position vs time using a distance constraint . . . 15
2.6 Spring - position vs time using distance constraint with high frequency and no damping . . . 16
2.7 Spring - position vs time using distance constraint and damping . . . 17
2.8 Spring - amplitude versus number of time steps for different frequencies . 19 3.1 Beam set-up . . . 25
3.2 Beam — maximum deviation vs time step . . . 27
3.3 Beam — max deviation vs time for time step of 9.8 · 10−5 s . . . 27
3.4 Beam — deviation for each box . . . 28
3.5 Beam — maximum deviation vs number of boxes . . . 28
3.6 Door set-up . . . 30
3.7 Door — deviation in total Fz compared to weight vs number of hinges . . 31
3.8 Door — Fx deviation for two hinges vs mass . . . 32
3.9 Door — deviation in total Fz compared to weight vs mass . . . 32
3.10 Falling Object — position vs time . . . 35
3.11 Falling Object — error vs time . . . 35
3.12 Falling Object — position vs time with altered time step . . . 36
3.13 Falling Object — error vs time with altered time step . . . 36
3.14 Friction set-up number 1 . . . 37
3.15 Friction set-up number 2 . . . 39
3.16 Friction — θ vs µs . . . 40
3.17 Friction — the force ratio vs tan(θ) . . . 40
3.18 Friction — anisotropy plot . . . 41
3.19 Pendulum set-up . . . 42
3.20 Pendulum — length deviation vs time step . . . 44
3.21 Pendulum — energy deviation vs time step . . . 44
3.22 Pendulum — energy vs time . . . 45
3.23 Pendulum — energy deviation vs time step, stiff system . . . 45
3.24 Pendulum — length deviation vs mass . . . 46
3.25 N-pendulum set-up . . . 47
3.26 N-pendulum — max deviation for each constraint . . . 48
3.27 N-pendulum — max deviation vs time step . . . 48
3.28 N-pendulum — energy deviation vs time step . . . 49 vii
3.29 N-pendulum — max deviation vs number . . . 49
3.30 N-pendulum — energy deviation vs number . . . 50
3.31 Planetarium set-up . . . 51
3.32 Planetarium — the orbit . . . 52
3.33 Planetarium — energy vs time, longer simulation . . . 53
3.34 Planetarium — distance vs angle . . . 54
3.35 Planetarium — energy deviation versus time step . . . 54
3.36 Cantilever beam set-up . . . 55
3.37 Cantilever beam — deflection vs position on beam . . . 57
4.1 Screen shot of the presentation of the results . . . 59
4.2 The result table of the continuous test . . . 60
4.3 Screen shot of the friction simulation on the web . . . 61
Acceleration a (ax, ay, az) m/s2
Angle θ rad or degrees,o
Compliance
Damping τ
Force F (Fx, Fy, Fz) N
Friction coefficient µ
Gaussian gravitational constant k = 2π/365 Gravitational acceleration g m/s2
Mass m kg
Momentum p N s
Number N
Pi π = 3.14159..
Position x (x, y, z) m
Spring constant κ N/m
Static Friction Coefficient µs
Time t s
Time step h s
Velocity v (vx, vy, vz) m/s
ix
Introduction
Simulations are an important part of today’s society. Not only for performance optimiza- tion, but also for training, education, research and games. It is a way to study systems that for some reasons are hard to examine in real life, for example scenarios in outer space or in construction purposes since it is often cheaper to simulate a system than to build a proper model. Simulations are also a way to construct virtual environment training systems used for both airline and military pilots, car drivers and do not forget its important role in the entertainment industry, especially for computer games [24]. It is also a way to gain access to deep information about systems, such as all positions and velocities. Many different frameworks for simulations have been developed, both commercial and open source. They are often referred to as physics engines.
A physics engine is a software program that keeps track of the laws from the physical world as well as handles forces and interactions between objects [5]. The engines are often divided into game physics engines and scientific physics engines where the first one focuses mainly on real-time simulations and the experience. For the scientific case it is more important that the calculations are precise for the accuracy of the simulations. This is the balance all engines face, the compromise between speed, stability and accuracy.
So far there are no standard for measuring the accuracy of a physics engine. Some efforts have been made to compare engines, often in order to choose for a particular task. Additionally benchmark tests have been developed but at present these concentrate mainly on game engines and examine speed and memory usage rather than accuracy and stability.
In this thesis a validation toolbox for a physics engine is developed. It examines the accuracy of the simulations through a collection of benchmark tests. These simulations check the engine’s performance in different areas such as constraint accuracy and stability, friction, stability of the integrator and overdetermined systems. Guidelines for choosing
1
benchmark tests are also produced. The following part of Chapter 1 describes the background of the thesis and a study of related work as well as goals, limitations of the work and methods. In Chapter 2 the theory of the engine is discussed and thereafter, in Chapter 3, each test case is presented with both theory, results and a short discussion.
Later, in Chapter 4, the additional work is displayed followed by a discussion and conclusions in Chapter 5.
1.1 Background
This thesis is initialized by Algoryx Simulation AB which develops software for physics based simulations. It was founded in 2007 as a spin-out from Ume˚a University and they are the producers of the AgX Dynamics. AgX Dynamics is a physics engine for simulating a wide range of scenarios such as wires and cables, rigid multibody system dynamics, vehicles and granular material for example [1]. Algoryx wishes a toolbox for validation of the accuracy of their product where focus is on non-smooth phenomena like impacts, contacts and dry friction for example. The aim is not only to validate what works but also to detect flaws in order to improve the product and additionally to be able to follow the performance of the simulations over time.
1.2 Related Work
There have been several studies performed with the aim to compare physics engines, e.g.
for choosing a suitable engine for a certain case. Seugling and R¨olin [38] published a master thesis where they perform nine run time tests in order to select a suitable physics engine to use in a physics module. They are testing friction, gyroscopic forces, energy conservation, constraint stability, accuracy, constraint scalability, scalability of contacts, piling and complex contacts. The benchmarks cover the most common scenarios in games and simple interactive simulations and most of the results are compared with analytical solutions. Still some of the more interesting cases are only judged using visual inspection and they do not put the engines to the test with a large amount of objects. In many complex simulation scenarios there are thousands of objects rather than 20. A typical training scenario for offshore simulations could be 2 to 3 ships, each with cranes on them with hundreds of bodies and cables that are kilometres long created by hundreds of elements. Also, due to the nature of the thesis, a part of the evaluation considers other aspects such as included features and cost rather than the accuracy of the engine.
An often referred work is the evaluation performed by Boeing and Br¨aunl [6] that is inspired from the study of Seugling and R¨olin. It is a qualitative study of free physics
engines using the Physics Abstraction Layer (PAL). The tests cover a free falling sphere, a bouncing sphere, friction with a box on an inclined plane, chain of spheres, a collision system and a stacking test. All tests are performed and chosen carefully but could be developed with larger mass ratios which can cause ill-conditioning, more objects and more complex scenarios to make things harder. Another similar study is the master thesis performed by Hansson [20]. It covers the basic benchmarks with a free falling box, collisions, a friction test and pendulum in order to test the constraint stability. All tests are varied with different friction coefficients, masses etc. Especially the mass is substantially increased in some of the tests. Even though this study is not that wide each test has a purpose and is thoroughly tested. Similar to the work of Boeing it would be useful to develop the tests further, towards more complex and realistic scenarios.
Mirtich and Canny [32] have published a work on an impulsed-based method for simula- tions. Despite that the article is not about evaluation of physics engines they perform several interesting benchmarks in order to test their method, especially for testing colli- sion detection. Some of the test results are compared to analytical solutions and some are visually inspected. These benchmarks may not alone found a decent validation of a physics engine but they provide a good starting point. This is the direction I followed.
1.2.1 Robotics
Robotics is another area where there have been different comparisons of physics engines.
As Roennau et al. [35] claim in their work does detailed simulations of robotic systems become more and more important in the development of the robotics area and its applications. They focus their tests on requirements for robotics which includes collision detection, complex collision detection, different stacking cases, constraint stability and two friction tests. In some tests the calculation time is investigated, in others the deviation from the analytical solution. Additionally they examine cases where the engines either fail or pass, for example if a boxes can be stacked into a pile without it collapsing. This is an interesting approach but it may contain weaknesses. For example, if it is only examined if the pile is collapsing or not the engine may pass the test even though the boxes are starting to penetrate each other or perform other unrealistic behaviour. In general the work is interesting and the cases are very well described which makes them easy to replicate. A similar comparison of physics engines for robotics is made by Erez et al. [10]. In this work an attempt is made to develop tests that are universal and differ from the regular benchmarks with rigid bodies in order to be more suited for evaluating engines for robotics applications. They test grasping of a capsule with a hand, a falling humanoid, a planar chain and 27 capsules falling on the floor. Over all the tests investigate which time steps the engines can use before the simulations fail and the CPU
time but also energy and momentum are examined when possible. The conclusion from the work is that each engine usually perform well on systems which it was designed for.
Even though the benchmarks are well performed and complex the work did not test common features such as friction and scalability. The test with the dropped capsules uses up to 250 capsules as an attempt to scale up the system but the simulation is only run for 10 seconds which may not be the case in real scenarios. Scalability is more deeply examined in the work published by Hummel et al. [21] where collision tests with 1 to 2500 objects are performed. This study also focuses on robotics and the additional tests are a bouncing tennis ball, a chain of spheres, interpenetration tests and a scenario where a screw is screwed into a nut. Apart from scaling up the scenarios this work also tries large masses. All tests are of importance when testing a physics engine but as they conclude in their work the benchmarks can be extended in order to check more features such as collisions with non convex geometries and soft body dynamics.
1.2.2 Offshore
The offshore exploitation industry is another field of study where simulations come in handy due the complex nature of the systems because the important aspect is wire and anchor handling, involving very long lines under very high tension in addition to hydrodynamics. One study performed by Metrikin et al. [30] compares different physics engines, mostly with respect to simulations of the broken-ice fields. The tests performed in this study are a free falling sphere, an ice floe colliding with a floater and energy dissipation when two ice floes are colliding. Additional tests such as friction tests and collision tests were suggested but not performed. All executed tests compared the error towards the time step to examine which physics engine was most suited for simulations of a broken-ice field. The study uses rather simple benchmarks which could be further developed for a more general evaluation, as suggested in the conclusion. Despite this the study is interesting since it deals with a specific and important area.
1.2.3 Web and Games
Apart from articles there are also web pages that present a set of test cases in order to evaluate physics engines. One of them is a product review by Lander and Hecker [27].
It lists 12 benchmarks which test collision detection, large size ratios, constraints, large mass ratios etc. The range of the tests are broad and especially the collision detection is well tested. Furthermore are each test focused on testing one thing which makes it easy to grasp the results. Unfortunately the results are not compared to analytical solutions and the grading of the engines is subjective. This could be a result of the work being mainly
focused on game engines. A similar set of tests is published at the web pages of MBS Benchmark, MultiBody Systems Benchmark [33], which is built on the work of Gonzalez et al. [14] [15]. It aims towards producing a standardized set of problems and procedures to measure computational efficiency and provide a platform where information can be collected, organized and shared. The main focus is multibody dynamics simulations and perhaps not physics engine principally. The simple pendulum, the N-bar mechanism, Andrew’s mechanism, Bricard’s mechanism and fly ball governor are the tests presented on the web page and all results are compared to analytical solutions. These cases are in the group Basic Problems but there is also a group called Industrial Applications.
Some of the tests listed above are very simple but there are also a stiff system and an over-constrained system. As many other studies mentioned the set of scenarios is a good starting point but should be developed further to really test the simulations.
One of the more extensive work to produce benchmark tests is called PEEL, Physics Engine Evaluation Lab, constructed by Pierre Terdiman. It is a collection of more than 250 tests that aims to evaluate and compare engines with focus on speed and memory usage. Some of the tests are described in his blog Coder Corner and they cover rigid bodies, stacking of different objects, joints, Raycasts, sweep tests among others [41]. As mentioned this collection of tests evaluates the speed and the usage of memory and it is mainly aimed at games applications. It does not consider the validation of the physics in the simulations.
1.2.4 Conclusion
To summarize, there are several studies performed where some kind of tests are executed.
The most common benchmarks are the simple pendulum, a free falling body, collision of two bodies and a chain of sphere, probably since they are simple to perform and the results can be compared to the analytical solution. These cases are of course important to test but it is crucial to keep in mind that if a simple scenario is accurately simulated it does not imply that also complex cases are correct. This is why the tests need to cover a wide range of cases, from simple ones to more complex ones. It is of course impossible to enclose everything but the studies discussed earlier may provide a starting point and a direction towards which tests that should be included. In most of the studies the aim is to choose a physics engine for a certain case or to determine which engine is the best. The aim of this thesis is not a comparison between different engines but rather a validation of the accuracy of the physics in the simulations today and a possibility to follow the change in performance during time. An additional goal is to develop a way to automate the tests and result production.
1.3 Purpose and Goals
The overall purpose of this thesis is to develop a product for validating the accuracy and stability of simulations performed by a physics engine and to create a set of benchmarks.
Since there are no similar work performed with the aim to create a method for validation that is universal the description of the assignment and the goals are flexible and may be changed during the project. The goals of the thesis consists of four points
• Create a set of benchmarks (10–20) for testing the accuracy of the simulations
• Design a strategy for creating new test cases which can be used with continuous integration and testing
• Develop a product that performs the validation of the physics engine
• Develop an export functionality for visualization of the results, either to the World Wide Web or as document, in collaboration with the staff on Algoryx
1.4 Delimitations
Since the number of possible things to test, both physical laws and numerical problems, are infinite the extent of the work needs to be limited in order to fit the scope of a master thesis (30 hp). This thesis is therefore restricted to 10–20 benchmark tests depending on the proportions of the tests. Additionally the export functionality and the possibility to use the Web is not developed from scratch but the existing tools is used together with newly developed tools produced by Algoryx.
Preliminaries
In this chapter there is an introduction to what a physics engine is and how it works, both regarding collision detection and the integrator. The numerical method Spook [24]
is briefly explained using the pendulum and the harmonic oscillator. In the beginning there is a reminder of the basics of physics and the chapter ends with some common features and problems when simulating.
Newton’s laws of motion
Newton’s three laws of motion are considered to be the fundamentals in classical mechanics.
The Newtonian laws describe how objects behave when exposed to forces, i.e. the response to the forces. Almost all textbook in mechanics states these laws, see for example Physics for Game Developers by Bourg [7, p.1]. Below are the laws stated in vectorial form.
Remember that the mass is assumed to be constant.
1. An object will continue to be at rest as long as the net forces acting on it is 0. If the object is in motion its velocity will not be changed. Mathematically it can be written as
XF = 0 ⇔ dv
dt = 0. (2.1)
2. The change of momentum (p) is equal in both magnitude and direction as the applied force if the mass of the body is constant, or
F = d(mv) dt = dp
dt. (2.2)
3. For a force FAthere is always a counterforce that is equal in magnitude but opposite in direction, i.e.
FA= −FB. (2.3)
7
The second law can be rewritten using the definition of momentum, p = mv. Since it is assumed that the mass is constant this gives
F = mdv
dt = ma, (2.4)
where F is the resultant force [7, p.15]. These equations are ODE’s of motion and need numerical integration. From here on, the time derivative will be denoted with a dot.
That is, the velocity will be written as ˙x and the acceleration as ¨x.
2.1 Physics engine
In real life bodies behave according to the Newtonian laws. In a simulation it is the physics engine’s responsibility to make sure that the laws are obeyed to make the simulation realistic. It is the engine that simulates these laws in the computer, it works as the heart of the simulation [9, p.295]. There are several ways to describe the structure of the engine, see both [9], [11], [20] and [22]. According to Erleben [11] the engine is described to consists of a collision detection part and a simulation part. Both sections will be described further in this chapter but shortly one can say that the collision detection handles the geometry of the state, i.e. analyses if there are any collisions or contacts while the simulation takes care of the physics related part with numerical integration, forces, constraints and so on.
2.1.1 Collision detection
Contacts and impacts are an important part of a realistic simulation. It is one of the ways that the objects can interact with each other. In real life the phenomenon is common, it happens every time you grasp something, stand or walk and it is natural that two solid objects do not occupy the same point in space. This is not the case in dynamics simulations.
First, before the program can handle a collision, the intersection between two objects needs to be detected, which requires a collision detection system. It is a pure geometrical problem of determining which objects are in contact with each other or will be in the short future time. This is a computational hard task. If a simulation has n bodies the computer needs to perform, in principle, n(n − 1)/2 checks in worst case. In order to reduce the work and speed up the process several methods have been developed. Usually the process starts with a broad phase that filters out objects into a list of potential collisions. This list then goes into the narrow phase that determine if there is a collision and also extracts additional information.
The Broad Phase
As mentioned above the broad phase aims to find objects that are close enough to each other for further testing [11]. A basic method for this is to use bounding-volume tests.
The idea is to use a simple geometries to cover the more complex objects, called bounding volumes [42, p.193], and the benefit is that it is easier to perform collision tests between the simpler primitives. Three different bounding volumes are shown in Fig. 2.1 inspired from Kimmerle [23]. It is the bounding sphere, the AABB (Axis Aligned Bounding Box) and the OBB (Oriented Bounding Box). The figure illustrates that the OBB is a better approximation than the sphere but in return it is computationally more costly. Examples of methods that can be used to test if pairs intersect are Exhaustive search, Sweep n’
prune and (multilevel) grids [11] as well as Space Partitioning and Model Partitioning which are described in for example Bergen [42].
Bounding sphere AABB OBB
Better approximations, higher build and update cost
Smaller computational costs for overlap test
Figure 2.1: Different bounding volumes, the bounding sphere, the AABB (Axis Aligned Bounding Box) and the OBB (Oriented Bounding Box) and their different properties regarding approximations and computational cost. Figure inspired from Kimmerle [23].
There are different ways to perform intersection tests. Firstly, to check if two objects intersect using only their instant position is called a static intersection test [17]. Secondly there is a dynamic intersection test which uses information about the motion of the object and constructs a sweeping volume instead of a bounding volume [11].
When the broad phase is performed the list of possible colliding pairs are reduced compared to the start. This list is then passed on to the next step in the process, the narrow phase.
The Narrow Phase
The goal with the narrow phase is to determine whether a collision is occurring or not.
Additionally other information such as contact points, penetrating features and proximity information can be gathered [11]. When using discrete time penetrations are unavoidable.
This phase also approximates the depth of these penetrations. The algorithms used for the collision tests in this phase are usually more refined but also slower compared to the ones in the broad phase. One algorithm that can be used is the Lin-Canny closest features algorithm which is described by Mirtich [31]. The advantage of the method is that it is fast, but on the other hand it has a problem terminating when two objects are penetrating and its convergence may be poor in degenerate cases. Another approach to the problem is to treat the object as a convex hull of a point set rather than focusing on polyhedral features. Gilbert, Johnson and Keerthi (GJK) have designed an algorithm with that approach [31]. The GJK algorithm computes the distance between convex objects iteratively. It is a fast and functional method that can be used on general convex objects. The drawbacks are that it is hard to understand and since it is iterative it is also sensitive to numerical errors [42, p.121]. There are many alternative methods such as the method of rotating calipers [9, p.464], the I-Collide method and the V-Clip (which are developed from both the Lin-Canny and the GJK methods) [9, p.464] [31].
2.1.2 Simulation
As mentioned earlier the simulation module is responsible for the physics in the simulations.
In AgX the numerical method Spook is used for solving the differential algebraic equations regarding the motion of the objects. For a complete description of Spook, see the PhD thesis by Lacouris`ere [24].
2.1.2.1 SPOOK - Pendulum
Spook is a discrete time variational method derived from the least action principle [24].
The action is defined as the time integral of the Lagrangian which is the difference between the kinetic energy and the potential energy. Shortly the principle can be explained as the fact that a body always want as small action as possible or the least action possible.
A thrown ball, for example, will always have a parabolic path in the air since it gives the least action. If the path would be irregular or bumpy this principle would be violated [4].
To simplify the description of Spook, let’s consider the pendulum. It is a bob of mass m moving at a distance l from the origin under constant gravity g (the bob is affected by the gravity force Fg). Often the motion is considered in 2D, Fig 2.2 shows an illustration of the scenario. The governing equations for the pendulum can be obtained from Newton’s second law, Eq. (2.4),
m¨x = Fc+ Fg (2.5)
where x is the position vector of the bob and Fc is the radial force keeping the bob on the path. Due to the attachment of the motion of the bob is restricted, it can only move on a trajectory a specific distance from the attachment point, l0. This constraint can be mathematically expressed as
g(x) = ||x|| − l0 = 0. (2.6)
This is a holonomic constraint since it is a strict equality and only depending on the coordinates and time [24, p.395]. It is the force Fc that constrains the motion to the allowed trajectory and it therefore needs to be directed radially towards the attachment point, i.e.
Fc= λˆx. (2.7)
In order to get something that is directed radially inwards the constraint is used. The steepest direction from the allowed trajectory g(x) = 0 is the gradient of the constraint,
∇g, denoted G. It is always perpendicular to the trajectory by nature. Therefore Fc = GTλ. From here on let Fg = Fe, for external forces. If all equations are put together it yields
m¨x = Fc+ Fg = GTλ + Fe
g(x) = 0. (2.8)
In 2D this is a system of 3 unknowns and 3 equations, i.e. a closed system, also called a differential algebraic equation of index 3.
g(x)
l
m Fc v
Fg θ
Figure 2.2: The theoretical set-up of the pendulum.
The system can be solved by a computer using discretization, both in space and in time.
Let xk be the position at time k and use the time step h which gives mvk+1= mvk+ hFk
xk+1= xk+ hvk+1.
(2.9)
Let Fk= GTkλh + Fe,k. Fe,k is the external force at time k and note the new definition of λ. Putting everything together yields
mvk+1 = mvk+ GTkλ + hFe,k
g(xk+1) = 0
xk+1 = xk+ hvk+1.
(2.10)
This is called Shake which is a method that is almost energy conserving. It preserves energy like quantity E = E0+ O(h) for exponentially long time [19]. In order to solve g(xk+1) = 0 the Newton-Raphson method can be used. Then a non-symmetric matrix appears which can make the problem very hard. Additionally it can fail to find a solution.
In Fig. 2.3 the different steps of Shake are illustrated, first mvk is added to xk. Then hFg is added and finally GTkλ. These three additions determine where the position xk+1 is situated and whether it fulfils the constraint or not. In the figure it does not cross the constraint circle which is a drawback of this method, that it sometimes may not find a solution. It happens when velocities are too high compared to constraint curvature.
There are projection methods that project back on the constraint surface but these requires some computational work and are still not guaranteed to find a solution [25].
g(x) xk
mvk
hFg
GTkλ
Figure 2.3: The steps of shake.
Spook is a way to handle this since it rewrites the constraint to 1
4(gk+1+ 2gk+ gk−1) + λ −τ
hGk+1(xk+1− xk) = 0 (2.11) where is the relaxation of the constraint and τ is the damping. That is, the constraint is not solved exactly but a small deviation is allowed, if the constraint is degenerate. is also called the compliance. Recall that xk+1h−xk = vk+1 and put all equations together
and it will give
mvk+1 = mvk+ GTkλ + hFe,k 1
4(gk+1+ 2gk+ gk−1) + λ + τ Gk+1vk+1 = 0
xk+1 = xk+ hvk+1.
(2.12)
Eq. (2.12) is Spook with holonomic constraints and damping in velocity format [24, p.98]. Also Gk+1 is approximated with Gk and gk+1+ 2gk+ gk−1≡ 4hgki is linearized around k. By doing this linearization and stabilization it is made sure that the constraint surface is always crossed, i.e. it is always possible to find a solution.
If the damping is 0, then λ can be calculated from Eq. (2.11) as λ = −1
4(gk+1+ 2gk+ gk−1) ≡ −1
hgki. (2.13)
Substituting Eq. (2.13) for the constraint force gives Fc= GTkλ = −1
GTkhgki. (2.14)
For the pendulum GTk is linear, i.e.
GTk = 1
||xk||xTk. (2.15)
Eq. (2.15) can be compared to the constraint force for a spring, Fspring= −κ∆x. It is proportional to the spring constant κ, a positive real number characteristic for the spring, times the deviation from the equilibrium length. Since this force is also linear and by comparing the forces it can be seen that the compliance (the relaxation parameter) is related to the spring constant as = κ1
Fc= −1
GTkhgki Fspring= −κ∆x.
(2.16)
To prevent division the solver is designed in a way that the division with is avoid but it illustrates the aspect of the theory. This relaxation corresponds to arbitrary stiff forces and the discretization makes them almost energy preserving, if there is no damping. It is a crucial part of the method since they prevent the simulation from exploding. This oscillations or crossing of the surface is similar to an idea by Hairer [18]. A method like this is often called a penalty method since it penalizes violations of the constraint, to keep it short.
In summary, Spook is a method that relaxes the constraint in order to be able to find solutions in all cases, including constraint degeneracy and ill-conditioning. This decreases the accuracy but makes the method stable. The simulation will oscillate around the
”true” constraint but with limited oscillations. Next an harmonic oscillator is considered in order to further investigate the damping parameter.
2.1.2.2 SPOOK - Harmonic Oscillator
Let’s consider a simple spring with an attached ball with the force Fspring = −κ∆x. The analytical solution for the spring, or harmonic oscillator, is
xz(t) = Ae−ζω0tcos(p
1 − ζ2ω0t) (2.17)
where xz is the position in z-direction, A is the amplitude, ζ = τ κ/(2√
mκ) is the damping ratio and ω0 =pκ/m is the undamped angular frequency [12], [43].
In order to get smooth and nice oscillations the following needs to be satisfied
hr κ
m < 2. (2.18)
The mass of the ball, m, is 1 kg. To satisfy the condition in Eq. (2.18) the spring constant is chosen to
κ = 0.01 4
h2. (2.19)
Additionally to simplify the scenario the gravity is put to 0.
First look at the case when the spring is only simulated using a force in every time step.
In this case there are no constraints, i.e. neither elasticity nor damping is present. Due to this Spook works just like the Verlet stepper
vk+1 = vk+ Fh xk+1 = xk+ vk+1h.
(2.20)
Fig. 2.4 shows the similarities between Spook and Verlet. It illustrates the position of the ball versus time for the two methods and it can be seen that the graphs follow each other completely. Actually the difference is always 0.
0 1 2 3 4 5
−2 0 2
Time [s]
Position,z
Position vs time, spring with only forces Simulation
Verlet
Figure 2.4: The position of the ball versus time both for the simulation and the Verlet stepper in the spring simulation.
Next the spring is simulated by a distance constraint with an compliance, = 1κ, also called constraint relaxation. The result is shown in Fig. 2.5. There it can be seen that the distance constraint still simulates the oscillating spring but with a phase shift. This case is similar to the constraint explained for the pendulum since it has no damping.
0 1 2 3 4 5
−2 0 2
Time [s]
Position,z
Position vs time, spring with distance constraint
Simulation Verlet
0 1 2 3 4 5
−1 0 1
Time [s]
Simulation minusVerlet
Difference between the spring with distance constraint and the Verlet stepper
Figure 2.5: The position of the ball versus time both for the Spook simulation using a distance constraint and the Verlet stepper in the spring simulation. Also the difference
between the methods is plotted.
It can also be interesting to see what happens when the spring oscillates faster than the condition stated in Eq. (2.18) using the distance constraint. Fig. 2.6 illustrates the position of the ball versus time and there it can be seen that the oscillations are uneven and not as smooth as the previous cases. This examples clarifies the importance of an appropriate time step in comparison to the velocity scale of the scene. If the bodies move too fast compared to the time step their true movement cannot be simulated. This is called aliasing, low frequency sampling of high frequency dynamics which gives a result that makes no sense. Also worth noticing is that even though the time step is too large for the velocities in the simulation the stability is still maintained.
0 1 2 3 4 5
−10 0 10 20
Time [s]
Position,z
Position vs time, spring with distance constraint, high frequencies
Figure 2.6: The position of the ball versus time in the spring simulation simulated with a distance constraint for the case of high frequency without damping. Stability is
maintained.
Finally numerical damping is imposed to the constraint to see its impact. The damping parameter, τ is chosen to τ = 2h which provides critical damping for undesired frequencies.
In Spook the system is linearized which can cause constraint violations. The damping is imposed to avoid these violations with the cost of energy conservation though not as much as for projection methods, for example [24]. A damped system, even if it is conservative, will not have conserved energy but in return the simulation is stable. A thing with Spook is that the decay of its damping can be predicted as √
αk when θ < 2 where
α = 1 − 4θλS
(1 + 4θ)(λS+ ). (2.21)
k is integer numbers from 0 and up, θ = τ /h and λS = 1 is an eigenvalue of a matrix S considering the constraint matrices and λs is also called the natural frequency [24, p.101]. Fig. 2.7 shows the result of the damped spring, both from the simulation and the analytical expression presented in Eq. (2.17). In the figure it can be seen that the
simulation follows the analytical expression but more importantly that it also follows the damping rate √
αk.
0 0.5 1 1.5 2 2.5 3
−1 0 1 2
Time [s]
Position,z
Pos vs time spring with distance constraint and damping Simulation
Analytical Damping
Figure 2.7: The position of the ball simulated with a distance constraint for the case with damping in the spring simulation.
It is possible to theoretically determine if a numerical method is stable by looking at the eigenvalues of the matrices involved. For the harmonic oscillator, consider the rewriting of Eq. (2.20) which become
"
I −hI
0 I
# "
xk+1 vk+1
#
=
"
I 0 0 I
# "
xk vk
# +
"
0 hF
#
. (2.22)
In this case the force expression is known which makes it possible to write the equations on the form Azk+1 = Bzk where
zk=
"
xk vk
#
(2.23) and A, B are matrices. Azk+1 = Bzk can be solved for zk+1 by multiplying with the inverse of A, i.e
zk+1 = A−1Bzk≡ Hzk. (2.24)
The stability of this system is totally determined by the eigenvalues of H. From now on the eigenvalues will be denoted as λ, not to confuse with the factor introduced to the radial force for the pendulum and in the constraint. In order for the system to be stable the eigenvalues need to be less than unity, i.e. λ < 1. Under certain conditions, see the thesis by Lacoursi`ere, [24, p.102], it can be seen that the equation for the eigenvalues reduces to
λ2− 2λβ + α = 0 (2.25)
where
β = 1 − 2(θ + 1)λS (1 + 4θ)(λS+ ) α = 1 − 4θλS
(1 + 4θ)(λS+ ).
(2.26)
For β and α the following holds:
0 < β < 1 when θ > 1/2, > 0 0 < β < 1/3 when θ > 2, > 0
(2.27)
0 < α < 1 when θ > 0, > 0. (2.28) Eq. (2.25) has the roots
λ±= β ±p
β2− α. (2.29)
The modulus of this is |λ±| = β ±pβ2− α. This can be examined in two different cases, β2− α < 0 for an oscillatory system and the case when β2− α > 0.
First, consider the case when β2− α < 0 which, as mentioned above, corresponds to an oscillatory system and occurs when θ < 2 and > 0. Then the modulus can be calculated as
|λ±| = r
β2+p
α − β22
=√
α. (2.30)
As stated earlier 0 < α < 1 when θ, > 0 and since it is the case in this regime
|λ±| =√
α < 1, i.e. the system is stable.
Next β2− α > 0 for which |λ±| = |α|. By examining different cases it is revealed that also this type of system will have a modulus less than unity and therefore is stable. For the complete analysis see the thesis by Lacoursi`ere, [24, p.103].
Finally, to further show some properties of Spook, let’s study the cases when is large and also when it is small or 0. The case with a large > λS corresponds to a low frequency system with the spring constant κ = 1/ ≤ 1. α can then be approximated as
α = 1 − 4θλS
(1 + 4θ)(λS+ ) ≈ 1 − 4θλS
(1 + 4θ)h2(1+4θ)4
= 1 −4h2θλ 4 =
= − hτ
= 1 −hτ
= 1 − hτ K = 1 − θ(hω)2
(2.31)
where it is used that the perturbation can be written as h2(1+4θ)4 , see [24, p.100]. In this case, a low frequency system, the aim is λ close to 1 in order to conserve energy and for stable simulations. As can be seen in Eq. (2.31) α is one minus something small if
the parameters are chosen correctly. That is, λ =√
α will be close to 1 in most cases for large and the simulation will conserve energy and be stable.
The opposite case, ≈ 0, is a high frequency system with h2ω2 1. If is put to 0, α can be written as
α = 1 − 4θλS (1 + 4θ)λS
λS (1 + 4θ)λS
=
= 1
1 + 4θ = 1 1 + 4τ /h.
(2.32)
That is, for = 0 the damping will be eigenvalue independent but the system is still stable. If a non-zero is introduced it will add a small filtering to the system. It is a low pass filter which removes high frequencies and represents the cut-off frequency. Fig.
2.8 shows this filtering property of Spook. The graph illustrates the amplitude of the oscillations of the spring for different frequencies versus number of time steps. It can be seen that higher frequencies are filtered more effectively than lower frequencies and the damping can even be calculated in number of time steps. For example, for τ = 2h, oscillations with a frequency mode higher than 50 steps per period will be all gone within 10 time steps.
0 10 20 30 40 50
−2
−1 0 1 2
Number of time steps
Amplitude
Amplitude for different frequencies vs number of time steps
Period = 100.0 steps Period = 50.0 steps Period = 10.0 steps Period = 5.0 steps Period = 1.0 steps Period = 0.1 steps
Figure 2.8: The amplitude of the oscillations of the spring for different frequencies versus number of time steps.
This simple spring example does not only explain how the stepper and the damping works but also shows the filtering property of Spook. Note that for low damping and low frequency, there is no numerical dissipation, as shown in Fig. 2.7. But choosing τ = 2h for constraint stabilization guarantees stability, without any tuning requirements and it is independent of masses. This example additionally provides a start for a stability
analysis. The stability of a method is extremely important but it is also a rather massive topic, huge enough to itself cover a thesis. It is therefore just briefly discussed in this chapter, for a deeper analysis see the thesis and the manuscript draft by Lacoursi`ere and Linde [24], [25].
In summary Spook is method that filters high frequencies. This high energy dissipation is similar to the property of the Generalized-α methods. They are stable implicit methods with second order accuracy but that requires solving nonlinear system of equations [37].
Also, as long as the ratio between the damping parameter and the time step is small enough, τ /h < 2, the system is not overdamped and the constraint will be satisfied eventually. These properties makes Spook a stable tool even for complex cases but with a cost of precision and some energy loss.
2.1.3 Key features
The everyday problem for a physics engine is the compromise between the following
• speed
• stability
• ”faithfulness”
• accuracy
• efficiency (accuracy versus work)
Speed
In order to create a realistic simulation the frequency of the images needs to be such that the human apprehend it as realistic. To achieve that the refresh rate needs to exceed 20 Hz. The standard for games is near 60 Hz which also is the frequency for visual comfort. In AgX the default time step is 1/60 s. This frequency gives little time for the computational work required of the engine which is why it needs to be fast.
Stability
An important part of simulations is the stability of the integrator. If the simulation is run for a long time time it is not eligible if the scenario explodes due to a constant increase in energy. If, of course, no energy is added on purpose. For a conservative system the energy should be constant or it could be oscillatory but with global bounds.
”Faithfulness”
”Faithfulness” means, in this case, that fundamental aspects of the laws of physics are
maintained as well as the damping introduced corresponds to physical damping. Shortly one can say that the simulation should aim to replicate the real world. The solutions provided by the physics engine is in fact the exact solution but of a slightly different physical problem of the same nature [19].
Accuracy
Especially for simulations used in investigation purposes are forced to be accurate in order to give the correct information. But accuracy is also important for other cases, for example in training simulations. The human mind have an incredible understanding of what is natural and reasonable even without having to measure things.
Efficiency
In many cases there is often possible to choose a numerical method that is very accurate but the cost is often the speed. Therefore it is wanted to have as high accuracy as possible with as little work as possible.
As mentioned the challenge is to compromise between these aspects. For example does a stable method lacks in speed in many cases. Is is therefore important to choose the correct method for the purpose depending on what is most important. Additionally it is of desirable to avoid tuning as much as possible, i.e. to be able to simulate a wide variety of scenarios without having to search for the perfect parameters [24][p.6-9].
2.2 Stable simulations
There are a large variety of things that can cause problems when simulating. In many cases the troubles can be problem specific and requires knowledge about the physics in the particular situation to solve. Some problems, though, are more general than others and in this section some of them are mentioned and what to do about it.
Choice of time step
There are many problems that can occur when the time step is not chosen properly.
Consider for example the spring mentioned above where the frequency is too high. This corresponds to choosing a too large time step in comparison to the velocity scale of the system. The engine then fails to correctly simulate the oscillations of the spring. Also, if a too large time step is chosen, the engine may not detect a collision between two objects which can result in a tunneling, i.e. two solid objects going through each other.
For example a ball falling towards a floor. If the time step is too large the ball may travel through the floor without the engine noticing it and preventing it from happening.
Some engine may have special methods for preventing this kind of events but not all.
The ball can also end up in the middle of the floor in one time step and this creates
large stabilization forces. These stabilization forces try to counter effect the constraint violation and since they are large the simulation can be unstable. Additionally the choice of time step can impose errors in the simulation. So, if the time step is too large it may cause
• Large constraint errors, i.e. inaccurate simulations
• Tunneling
• Large interpenetrations which causes large stabilization forces and therefore unstable simulations
The advantage of a large time step is that the simulation is fast. If the time step is too small the simulation instead becomes slow [39].
Initial state
A properly chosen initial state provides a good start for the simulation. If two solid bodies are initialized penetrating each other it will result in large forces and therefore also high velocities. Other invalid initial states are unrealistic forces on motors for example or unrealistic velocities. Additionally is it not good to initialize a system with large constraint violations. All of these scenarios can result in unstable simulations [29].
Mass ratio
Another important thing to consider in simulations is the mass ratio. Very large mass ratios can result in ill-conditioned matrices for the method to solve which can lead to undesirable simulations. An example of a bad mass ratio is to attach an object of 109 kg to a thin wire which weigh around 10−5 kg. This gives a mass ratio of 10−14 which can be compared to machine precision that is 10−16. That means only 2 digits of accuracy which will give all wrong results. It is therefore impossible for the computer so solve a problem of this kind. [34].
Solver
In previous sections the numerical method Spook and its key features have been described.
It ended up with a system of equations that needs to be solved which can be done in different ways. First there is the direct solver which solves the system ”exact” except for rounding errors. The Gaussian elimination, for example, can be used for solving linear system of equations. In AgX the direct solver is called Saber [26]. Next is the iterative solver which generates an approximation of the solution and then improves the approximation until a termination criteria or the maximum number of iterations are reached. An example of this method is the Gauss Seidel method. The iterative solver can have problems with stiff systems and large mass ratios where the simulation may end up
with rubber-like behaviour. But for non-linear cases the direct solver may have a hard time solving the problem and then it can be a good idea to use the iterative solver [13].
For AgX dynamics the non-linearity is not a problem, since everything is linearized. Its direct solver can have a problem, though, when the simulation contains many bodies that are tightly connected.
The default solver in AgX is a hybrid solver which solves all constraint with the direct solver, except the frictional contacts. It uses a direct method for joints and normals and is coupled with iterative methods for friction, see Lacoursi`ere and Linde for details [25].
Test cases
This chapter first describes the work flow of the thesis. Then the different test cases are presented with both model set-up, a short theory introduction and the specifics of the simulation. It is followed by the results of the test and all cases are summed up with a short discussion.
3.1 Work flow
3.1.1 Pre-study
In the beginning of the work a pre-study was performed with information gathering from articles, master’s dissertations, theses and web pages. The goal was to understand how a physics engine works in general and which scenarios that can cause problems.
Furthermore it was investigated what has been done earlier regarding evaluation of physics engines and which benchmarks that usually are used. The aim of the thesis is to cover as wide range of test cases as possible but also to sort out the most useful ones for the company. To be able to prioritize the test cases a number of interviews with stakeholders at Algoryx were performed. Then it was discussed which cases that are suitable and interesting for the company.
To be able to implement the platform an introductory course in AgX Dynamics was taken since it is the engine provided for the work. The version 2.14.0.0 of AgX Dynamics is used.
Additionally, for writing scripts, some effort was put into learning the scripting language Lua [28] which can be used to create simple simulations with AgX. The pre-study ended up as a project plan.
24
3.1.2 Implementation
When the pre-study was finished the implementation of the platform began. The simulations are written as script in Lua and the analyses are performed with Octave 3.8.2 which also produces the plots for the HTML-pages. First a few easy benchmarks were implemented in order to learn how the programs work. Those were then evolved by adding more complex scenarios. Simultaneously an export function was developed. Due to limited time this was restricted to the web as an HTML-page. In parallel with the implementation, the strategy for choosing and creating benchmark tests was developed as well as the export and visualization functionality. Additionally, after the first few test cases were almost done, the work with the report started.
When a couple of tests were done it was investigated how continuous tests could be performed. Scalar trials for some of the test cases were made and put together in a continuous test with additional result presentation. Also a prototype of a Web application was tested in cooperation with the developers on Algoryx.
3.2 Test cases
Below is a presentation of all tests performed with the set-up, the results and a discussion.
All simulations are run with the default solver, i.e. the hybrid solver.
3.2.1 Beam — Constraint stability
The stability of lock constraints is the topic of this test case. It is examined by simulating a beam situated on two piles. The beam is constructed using a number of boxes, i.e. it is discretized as a chain of rigid bodies constrained to each other with lock joints. Fig.
3.1 illustrates the set-up of the model.
x z
Figure 3.1: The set-up of the beam created by boxes situated on piles.
As can be seen in the figure the beam is situated on two piles which exert a lifting force on the beam to counteract the gravity pointing in the negative z-direction. The wish is
as small height deviations as possible. In order to have stable simulations the damping needs to be non zero, otherwise the beam started to oscillate and finally explode. In the simulation both the constraint and the contact between the piles and the beam is controlled regarding compliance and damping.
Two tests are performed, one where the maximum deviation in height is measured as a function of time step and one where the number of boxes in the beam is varied. All simulations are run for 1000 time steps for a fair comparison between the tests and all objects are initially at rest. Due to the oscillations that may occur in the beginning of the simulations the maximum deviation is measured for the last 50 time steps.
The parameters for the first simulations are:
• Mass/box: 2.5 kg
• Size of box (x, y, z): (0.125, 0.1, 0.1) m
• Number of boxes: 80
• Gravity (g): 9.81 m/s2
• Compliance (): 0
• Damping (τ ): 2h
In this simulation the beam reaches a total length of 10 m, which is common length at the steel beam market [40].
The second simulation investigates the effect of the number of boxes, i.e. the length of the beam. The number goes from 10 up to 200 where 200 represents a beam of 25 m and the simulation is run for 20 seconds. All other parameters are the same as in the previous simulation and the time step used is 1/60 s.
3.2.1.1 Results
Fig. 3.2 illustrates the maximum deviation on the beam compared to the total length versus the time step. As can be seen in the graph the deviations are overall small, smaller than 0.1 ppm, except for a time step of about 10−4 s. For that time step the deviation is approximately −0.2 · 10−6%. In the next figure, Fig. 3.3, the maximum deviation for the beam simulated with the time step of 9.8 · 10−5 s is plotted versus time. There it can be seen that the beam never stops oscillating which is the reason for that increased deviation. For the other beams the oscillations eventually stopped and the simulations
stabilized. In the lower plot in Fig. 3.2 the deviation for larger time step is visualized showing that for time step larger than 0.01 s the deviations start to increase. Since the deviation is negative it implies that the beam is too low compared to the ideal height, i.e. it has dropped down.
10−4 10−3 10−2 10−1
−0.2
−0.15
−0.1
−5 · 10−2 0
Time step [s]
Maxheightdeviation[10−6%]
Maximum height deviation compared to total length vs time step
Figure 3.2: The maximum deviation in height of the beam versus the time step.
0 2 · 10−2 4 · 10−2 6 · 10−2 8 · 10−2 0.1
−1
−0.8
−0.6
−0.4
−0.2 0
Time [s]
Maxheightdeviation[10−6 %]
Max height deviation vs time, h = 0.000098 s
Figure 3.3: The maximum deviation in height of the beam versus the time for time step of 9.8 · 10−5 s
In the next figure, Fig. 3.4, the default time step is used and the deviation for each beam element is plotted. It shows that the beam bends like a half moon, as expected, with the largest deviation in the middle. The largest deviation is about −5 · 10−9 m and the bending seems symmetrical.
0 20 40 60 80
−5
−4
−3
−2
−1 0
·10−3
Position on beam Heightdeviation[10−6m]
Height deviation vs position on beam
Figure 3.4: The deviation for each box in the beam, i.e. the position on the beam, with default time step, 1/60 s.
In Fig. 3.5, the last figure, the maximum deviation of the beam height compared to its total length is plotted versus the number of elements in the beam. Since the length of each beam element and their weight are kept constant both the total length and the total mass of the beam are increased as well. In the figure it is shown that the deviation decreases rapidly with increasing number of elements, especially for 100 and more elements which is longer than approximately 12 m.
101 102
−8
−6
−4
−2 0
·10−2
Number of elements Maxheightdeviationcompared tototallength[10−6%]
Maximum height deviation vs number of elements
Figure 3.5: The maximum deviation in height compared to the total length versus the number of boxes constructing the beam. All values are negative.
