REBOUND Hard-Sphere Fails for Long Contacts Under Gravitation

[ErickWhiteDev]

Environment
Which version of REBOUND are you using and on what operating system?

• REBOUND Version: 4.4.10
• API interface: Python
• Operating System (including version): Linux 24.04, running in WSL 2 on Windows 10

Describe the bug
The default REBOUND hard-sphere implementation, while effective for brief collisions, breaks down under periods of extended contact. This can be seen in, e.g., integrations of two rigid bodies resting against each other; assuming a coefficient of restitution e < 1, after the relative velocity between the spheres drops low enough, rather than resting against each other (as would be expected for rigid bodies), the two begin to sink into each other.

To Reproduce

import rebound

def coefficient_of_restitution(r, v):
    return 0.5

sim = rebound.Simulation()
sim.G = 1.0
sim.collision = "direct"
sim.collision_resolve = "hardsphere"
sim.coefficient_of_restitution = coefficient_of_restitution

sim.start_server(port=1234)

sim.add(m=1, r=2, x=10)
sim.add(m=1, r=2, x=-10)

sim.usleep = 16667

sim.move_to_com()

sim.integrate(float('inf'))

sim = None

This code snippet shows an example of the sinking behavior in a simple two-particle case with e = 0.5.

Additional context
@dangcpham and I discussed this issue in a meeting today; my current method of circumventing this involves reverse-time integration upon detection of a collision, followed by a bisection to determine the actual time of collision. This conserves energy very well (typically to within around 1E-13 relative error), but is computationally expensive.

Name/Affiliation
Erick White, Celestial and Spaceflight Mechanics Laboratory, Colorado Center for Astrodynamics Research, University of Colorado Boulder

[dangcpham]

I think it might be helpful to implement a REB_COLLISION_EXACT method to find the exact time of collision (following a bisection method as you mentioned). If we can find the time of contact, the velocity at that point would be sufficiently high to prevent the sinking bodies

[hannorein]

This is a well known effect. There are a couple of ways around this. The easiest way is to enforce a minimum collision velocity (see https://rebound.readthedocs.io/en/latest/simulationvariables/#collisions). This works well in dense planetary rings with many collisions per particle/timestep. The minimum collision velocity needs to be fine tuned for the specific setup. And of course this formally violates energy conservation. But if the collisions are dissipative anyway, this should not be a big deal.

Another way is what you try to do, find the exact time of collision. However, this scales badly, O(N^3) in the worst case. If there are many collisions per particle/timestep, then you need to find the first collision, resolve that, then search for the next one in the remainder of the timestep (any previously found collision might have changed), and so on.

Other people implement hard spheres as a repulsive force. This is nice because no actual collision detection needs to be performed in principle. But getting this right is tricky as it makes the equations of motions much more stiff and might lead to unwanted effects, for example vibrational modes.

[dangcpham]

Thanks for the background Hanno! I was thinking to use the line collision algorithm until a collision is detected, then refining the collision time for the 2-body collision with the bisection method. This way during most time steps, collisions are checked faster with the line algorithm (presumably O(n^2)?)

[hannorein]

What system are interested in?

For dense planetary rings, this is not a good solution as it scales at O(N^3). Use the minimum collision velocity instead.

If you have a few planets and want to know the exact time of collision then the sort of approach you're suggesting is good.

If you have some system where two particles are constantly touching each other, then neither approach is good and you'd need to look into the force model.

[ErickWhiteDev]

At the moment, the system consists of three equal-mass, equal-radius bodies which often do end up in contact over long periods of time. I was worried about the effectiveness of a hard-sphere model in simulating this, but using a force based or soft-sphere model led, as you suggested, to stiffness (and typically to non-physical increases in overall system energy and angular momentum).

The current approach uses a velocity cutoff in addition to that bisection method; this is undesirable as it is a bit of an arbitrary number currently (based off of a similar constraint detailed in Gabriel and Scheeres (2016)) and can affect the long-term system angular momentum accuracy.

I think the best bet here may be to take another shot at the soft-sphere method if a force-based model, rather than an impulse-based one, is most likely to be successful.

Thank you for the advice!

[hannorein]

You might also want to look at Alice Quillen's code

https://arxiv.org/pdf/1901.01439v3
https://github.com/aquillen/rebound_spring