Modia3D: Modeling and Simulation of 3D-Systems in Julia

Modia3D is an experimental Julia package to model and simulate 3D mechanical systems. Ideas from modern game engines are used to achieve a highly ﬂexible setup and features of multi-body algorithms are used to get a rigid mathematical formulation and support, for example, of closed kinematic loops. Collision handling is performed on convex geometries with elastic response calculation. A Modia3D model is solved with a variable-step solver. This is important to combine Modia3D with the equation-based modeling system Modia in the future.


Introduction
The open source Modia prototype platform, see table 1, consists of various packages based on the Julia programming language [3] to model and simulate physical systems described by differential and algebraic equations. The goal is to simulate, e.g. robots, vehicles, aircraft, buildings or power plants, and to experiment with novel features for the next Modelica © language generation 1 . Modia3D 2 is an experimental modeling and simulation environment to provide 3D geometry and 3D mechanical systems for the Modia platform. The goal is to fully integrate Modia3D with Modia's equation-based modeling and provide a common graphical user interface with the web app Modiator. This will allow to model, for example, the 3D-mechanics of a robot with Modia3D and the electrical motors, gear boxes, sensors and controllers with Modia. One enhancement with respect to the widely-used Modelica language will be that specialized algorithms for 3D kinematics and mechanics, combined with equation-based modeling, allows more robust and efficient simulations of complex systems. Modia3D uses ideas of modern computer game engines 3 , to achieve a highly flexible setup of mechanical systems including collision handling. Contrary to game engines, numerical integration is performed with a variable-step solver (IDA via the Sundials.jl Julia package [9,16]) as needed for applications where system simulations have to match reality with a certain precision. In Modia3D, collision handling is performed with elastic response calculation and not with impulses, as it is common for game engines. The reason is that simulation results are closer to reality and it is easier to On the other hand, multi-body programs are usually far from the flexible setup of games (see also section 2). Modia3D provides a generic interface to visualize simulation results with different 3D renderers. Currently, the free community edition as well as the professional edition 5 of the DLR Visualization library 6 [1,8] are supported. Another team is developing Modiator, a free 2D/3D web-based authoring and rendering tool. The user's view of Modia3D was introduced in [12] showing the very flexible definition of 3D systems. Some key algorithms are discussed in [11,13]. Collision handling with elastic response calculation and error controlled integration is challenging and this article discusses some of the difficulties and how they are solved. Further, it gives an overview of Modia3D from a user's perspective, and in particular how collisions between objects are defined.

Flexible Definition of 3D Systems
Modia3D follows the approach of modern game engines to provide a coordinate system as a primitive that is located in 3D and has a container with optional components (such as geometry, visualization, dynamics, collision properties, light, camera, sound, etc.), see for example [14] 7 . 8 Such types of objects are called GameObject 9 in Unity, Actor 10 in Unreal Engine, and Object3D 11 in Three.js. In Modia3D the name Object3D is used. This very flexible approach allows to define many optional components and variants and treat them in a modular way. The Julia programming language is particularly suited for this component-oriented programming pattern and therefore key-concepts of Julia, such as multiple dispatch, are heavily used in Modia3D. Hierarchical structuring for grouping and aggregation is performed with the Modia3D macro @assembly. Julia macros are metaprogramming 12 language elements and a macro name starts with @. It generates an abstract syntax tree (AST) of Julia code which is automatically compiled and executed at the line where the macro is called. Fig. 1 shows a bar that is constructed from several Object3D elements and is defined with the following Julia code: § ¤  The first Object3D is a SolidBeam geometry that is made of aluminium and is visualized in blue color. Mass, center of mass and inertia matrix are computed internally from this information. The next two objects place coordinate systems on the first object. Three instances of the Bar assembly together with a ground object (= the fourth bar) are used to build up the four-bar mechanism 13 shown in Fig

¦ ¥
The angle of revolute joint rev1 is moved kinematically via a signal. The resulting system is simulated by calling function simulate! (for more details of the Modia3D elements, see [12])

Collision Handling
Modeling and simulating the collisions between geometrical objects is difficult and there are many methods dedicated to particular purposes. For example, in a game it is important that collisions of many objects are supported in real-time and that a simulation looks reasonably realistic. Game engines typically make compromises regarding physics (see for example [5], where it is pointed out that NVIDIA PhysX and Havok neglect Coriolis forces). When designing an industrial product, it is important that collision models can be validated by measurements. Therefore, the simulation results must be much closer to reality as in a game.

User's View
Objects only take part in collision handling, if a contactMaterial is associated with them. Modia3D uses two sets of material data: (1) Solid material constants In a dictionary, the name of one material is used as a key, for example Steel. The values are typical material constants of a solid. For elastic response calculation, Young's modules E and Poisson's ratio ν are used. From this data of two contacting objects, a spring constant is computed to describe the elastic deformation in the contact area.

(2) Contact pair material constants
In another dictionary, the names of two materials are used as a key, for example Steel and DryWood. The values are the coefficient of restitution cor, the kinetic/sliding friction force coefficient µ k and the rotational resistance torque coefficient µ r between these two contacting objects of the respective materials, see section 3.3.
A simple example is shown in Fig. 3, where a steel ball is bouncing on a wooden table. This model is defined with the following Julia program (r is the position vector from the parent -defined with the first argument -to the reference frame of the object and fixed defines whether the object is fixed or not fixed in space): Proceedings of JuliaCon

Numerical Solution
A Modia3D model is mathematically defined by the hybrid DAE system (1), where x = x(t) and J (1c) is a regular Jacobian: When differentiating f c once, it is (conceptually) possible to solve (1c) forẋ because J is regular. Therefore, (1a) is an index 1 DAE. (1b) defines zero-crossing functions z(t). Whenever a z i crosses zero, an event is triggered, simulation is halted, functions f d , f c can be changed, and simulation is restarted. (1) is numerically solved with the variable-step DAE integrator IDA of the Sundials suite [9] via the Sundials.jl [16] Julia package. Since contact forces and torque lead to extreme changes of f d , f c , the distances between convex shapes are used as zero-crossing functions z i (t) so that the start and the end of a contact phase triggers an event. It is expected that this approach improves the reliability of a simulation. Furthermore, the elastic response calculation of section 3.3 needs the relative velocity when a contact starts, to compute an appropriate damping factor d. Therefore, an event at the start of a contact is mandatory.
Distances between convex shapes, as well as penetration depths are computed with an improved version of the Minkowski Portal Refinement algorithm (MPR-algorithm) [18]. The MPR-algorithm is much simpler to implement and has less numerical problems than the often used GJK/EPA-standard algorithms [2,7], because it only works with triangles and not with tetrahedrons. In the original version of the MPR-algorithm [18] only penetration depths are determined. In Modia3D improvements of the MPR-algorithm are utilized that have been proposed in [10,11], in particular to compute the distances of shapes that are not in contact and treat special collision situations properly.  Fig. 4. Contact normal force f n , contact tangential force f t (= sliding friction force) and contact torque τ ω between two penetrating objects. e n , e t , e ω are unit vectors in direction of the respective relative movement.
Collisions of n potentially colliding shapes are handled in the following (mostly standard) way: 1. Broad Phase The shapes are approximated by Axis Aligned Bounding Boxes, see e.g. [2], where potential collisions and approximate distances can be determined with low computation cost resulting in O(n 2 ) low-cost tests. When using special data structures (such as octrees or kd-trees), it is even possible to reduce the number of low-cost tests to O(n log(n)).

Narrow Phase
For the potentially colliding shape pairs as identified in the broad phase, the signed distances are computed with the improved MPR-algorithm [11].

Response Calculation
If two shapes are penetrating, a normal and a tangential force, as well as a torque are applied at the contact point. For details, see section 3.3 and Fig. 4.

Elastic Response Calculation
A contact normal force f n , a contact tangential force f t (= sliding friction force), and a contact torque τ ω are calculated when two 3D objects penetrate each other with a penetration depth δ ≥ 0, as shown in Fig. 4. The intuition is that there is a contact area with a certain pressure distribution in normal and a stress distribution in tangential direction and that the response characteristics provides an approximation of the resultant force and torque of these distributions.
The MPR-algorithm computes the contact point, δ, and a unit vector e n that is orthogonal to the contacting surfaces. The novel response characteristic of (2) utilizes ideas from [6,15,17] with some extensions and corrections. Variables with index geo are computed from the contacting geometries and vectors with index reg are regularized with smooth characteristics to avoid divisions by zero (for example, | e t,reg | = 1 with exception of a small region around | v rel,t | < v small , where | e t,reg | = 0 at | v rel,t | = 0 and it smoothly approaches 1 at | v rel,t | = v small ): The symbols of (2) have the following meaning: e n Unit vector normal to the contacting surfaces. e t Unit vector in direction of the relative tangential velocity. e ω Unit vector in direction of the relative angular velocity. c res (E 1 , E 2 , ν 1 , ν 2 ) Resultant spring constant in normal direction. c res = 1/(1/c 1 + 1/c 2 ); d(cor reg ,δ − ) Damping coefficient as function of the regularized coefficient of restitution cor reg andδ when contact starts (δ − ≥ 0). Variable cor reg is computed according to [6] and regularized so that cor reg = 0.001 if the normal relative velocity v rel,n = 0 m s −1 and approaches cor when v rel,n = v small . One reason for this is that otherwise [6] would result in a division by zero if cor = 0. The other reason is that a bouncing object stays at rest if v rel,n becomes small enough and therefore cor reg must be drastically reduced for small v rel,n (this effect is typically not taken into account in other contact laws, such as in [17]). µ k Kinetic/sliding friction force coefficient (≥ 0). µ r Rotational resistance torque coefficient (≥ 0). Its effect is that torque τ ω is computed to reduce the relative angular velocity ω rel between the two objects until ω rel = 0 rad s −1 . For a ball, µ r is the (standard) rolling resistance coefficient and µ geo is the ball radius.
c geo , n geo , µ geo depend on the geometries of the two objects in contact. If not enough information is available, these factors are set to one. c geo , n geo take the contact volume into account, under the assumption of Hertz' pressure (e.g. n geo = 3/2 if at least one of the two contacting objects is a sphere). The response characteristics (2) shall be clarified with a few special experiments: Comparison between elastic and impulsive contact response. In Fig. 5, the height of the bouncing ball (see model of Fig. 3) is shown as a function of time. The red curve is the reaction when using the elastic response calculation of (2) with cor = 0.7 and the solid material constants of steel and of dry wood. The blue curve is the result when the contact force is computed with an impulse for the same cor value. The green curve is cor reg (for small velocities it becomes small). This shows that f n in (2) leads to a similar reaction if compared to an impulsive response. Note, this was the intention for the development of this force law in [6]. A systematic comparison was made in [17].
Sliding and rolling ball. In Fig. 6 an animation of a billiard ball is shown that is sliding and rolling on a billiard table. The billiard ball is a free flying object with 6 degrees-of-freedom where the rotation  and ideal rolling of the ball takes place. Since v rel,t = 0 m s −1 , the sliding friction force f t = 0 N, because e t,reg = 0. Therefore, the ball would roll forever, contrary to reality. On the other hand, the rotational resistance torque τ ω (µ r = 0.02) acts as a rolling resistance that continuously reduces the angular velocity ω ball,2 and the ball comes to rest at some point in the future outside of the plot area (see right upper and lower plot of Fig. 7). The result is that (2) is able to reproduce the effect of a sliding and rolling ball. However, for simplicity of the formulation, in this simplistic elastic response law, velocity dependency of the coefficients cor, µ k , µ r is neglected.
Collision of two balls. Two billiard balls are positioned on a billiard table, see Fig. 8. One of them has an initial velocity and hits the other resting ball after some time. Simulation results are shown in Fig. 9. Before the first ball hits the second one, the effects of sliding and rolling occur as analyzed before. At time = 0.55 s, the first ball hits the second ball. Since the coefficient of restitution between the two balls is one, a fully elastic collision takes place. In this case the first ball transfers most of its kinetic energy to the second ball which starts moving with the velocity of the first ball. This "exchange" of velocity can be observed in the middle plot of Fig. 9 at time = 0.55 s. However, since the first ball was rolling, the angular momentum was greater zero. This momentum is conserved. Therefore, the first ball continues rolling and velocity v ball,1 rises from zero again. Since the relative velocity is no longer zero due to the impact, again a friction force f t is acting that introduces a counter torque at the balls axis which quickly reduces the angular velocity until again the relative velocity is zero around time = 0.65 s. Both balls are again ideally rolling, and due to the rotational resistance torque, the angular velocities are slowly reduced until both balls come to rest which is not shown in the plots.

Regularizing the Contact Area
A DAE solver such as IDA, see Section 3.2, solves a nonlinear algebraic equation system at every time instant. If f d , f c in (1) are not smooth or even not continuous, it is most likely that no solution of this algebraic equation system is found and simulation stops. This situation occurs easily, for example, if the edges or vertices of two boxes collide on each other, because very small changes of the positions of the boxes can change the contact situation drastically. To improve the reliability of the simulation, Modia3D uses the approach proposed in [2], and is smoothing every shape with a small sphere (default radius = 1 mm) that is (conceptually) moved over all surfaces. Practically, this smoothing can be very easily and cheaply incorporated for the convex collision handling. Every method that considers only contact points on two colliding objects and the corresponding penetration depth has the disadvantage that an infinite number of solutions can occur in some cases. For example this happens, when a box falls on a table and the faces  of box and table are parallel. Imagine small changes in the configuration can result in identification of very different contact points in the contact area. It is then highly probable that simulation with a variable-step solver fails. Furthermore, the contact force and torque depend heavily on the size of the contact area. The penetration depth alone does not provide enough information to calculate reasonable values that are in accordance with physics. For both cases, it would be necessary to take the complete contact area and volume into account, as done e.g. in [4]. This method is however much more costly and probably only works reasonably for soft contacts where the contact volume is large enough. To summarize, collision handling with variable-step solvers and penetration depth computation can most likely only work reliably for point contacts that are smoothly changing.

Example: Billiard Table with 16 Billiard Balls
The billiard table in Fig. 10 has 16 billiard balls. The material constants are shown in Table 2. The cue ball has an initial velocity pointing to the right and hits the the center of the rack (15 other balls) exactly after a short time interval. This results in a symmetric evolution of the balls, as one would expect. All previously described effects (sliding, rolling, colliding) act together. The hybrid DAE system has dim(x) = 13 · 16 = 208 and there are about 200 possible collision pairs. Simulation on a standard PC needs about 20 min for 5 s of simulation time. At the moment, the Modia3D code is implemented for functionality and not tuned for efficiency, so a speed-up is expected in the future.

Conclusion
In this article, a short overview about the experimental 3D modeling environment Modia3D is given. In particular, collision handling with a variable-step solver has been sketched and a novel formulation for elastic response calculation is proposed. The Modia3D collision and contact handling is demonstrated with several examples. Modia3D combines ideas from different communities. The architecture with component-oriented modeling is inspired by game engines so that 3D models can be setup in a very flexible way, as well as several elements for collision handling. Other features are from multi-body programs, like hierarchical structuring, support of closed kinematic loops, and algorithms to compute results close to real physics. Modia3D is still a prototype implementation and several important parts are under development. Especially, the integration with Modia is missing at the moment. Furthermore, the code was currently mainly developed for its functionality and is not yet tuned for efficiency. For these reasons, benchmarks and comparisons with other programs with respect to simulation efficiency have not yet been performed.