[FIXED] What on earth is going in in my particle simulation?

Issue

I've got my colliding balls simulation kind of working but the particles go crazy.

Some of them stick to each other and spin around. others kind of get stuck at the edges.

i can't think of why... could anyone help me out here? thanks

This code uses the Taichi Gui. I intend to use the taichi decorators in the future but for now I'm just testing it with regular python

import taichi as ti
import numpy as np

ti.init(arch=ti.vulkan)

gui = ti.GUI('Circles Test', res=(600, 600))

numpoints = 20
centres = np.random.random((numpoints, 2))
velocities = np.random.uniform(-1,1, (numpoints,2))
radius = 20
guiradius = 0.04
damping = 0.9
timestep = 1/60

'''
@ti.kernel
def randompoints():
    for i, j in centres:
        centres[i,j] =ti.random()
'''
#randompoints()
print(centres[1])

def wallcollision(centres, velocities):
    for i in range(len(centres)):
        if centres[i][0] - guiradius <=0:
            velocities[i][0] = -damping*velocities[i][0]
        if centres[i][0] + guiradius >=1:
            velocities[i][0] = -damping*velocities[i][0]
        
        if centres[i][1] - guiradius <=0:
            velocities[i][1] = -damping*velocities[i][1]
        if centres[i][1] + guiradius >=1:
            velocities[i][1] = -damping*velocities[i][1]
    return velocities

def selfcollision(centres, velocities):
    for i in range(len(centres)):
        for j in range(i+1,len(centres)):
            if np.linalg.norm(centres[i]-centres[j]) < 2*guiradius:
                print("Collision!")
                vitemp = np.copy(velocities[i])
                vitemp = velocities[i] - (np.dot(velocities[i]-velocities[j], centres[i]-centres[j])*(centres[i]-centres[j])/(np.linalg.norm(centres[i]-centres[j])**2))
                velocities[j] = velocities[j] - (np.dot(velocities[j]-velocities[i], centres[j]-centres[i])*(centres[j]-centres[i])/(np.linalg.norm(centres[j]-centres[i])**2))
                velocities[i] = vitemp
                
                ##velocities[i] = -velocities[i]
    return velocities

def update(centres, velocities):
    centres += velocities*timestep
    velocities = wallcollision(centres, velocities)
    velocities = selfcollision(centres, velocities)

while gui.running:
    gui.circles(centres, radius, color=0xEEEEF0)
    gui.show()
    update(centres, velocities)

Solution

The issue is that when a collision is detected, the velocity is modified but not the centre, and the velocity is not always big enough for particles to escape the collision area in the next frame. This means particles can stay in collision during many frames. This causes the velocity components to oscillate between negative and positive values making the collision escape sometimes impossible. This is why particles are somehow glued together. The same issue happens on the wall : when the centre of a particle is beyond the walls, the velocity components can oscillate so the particle is glued to the walls. This kind of issue arises quite frequently in basic physics simulations. A simple fix is simply to correct the centres so for the particles not to collide in the next frames.

Here is the corrected code:

# My parameters:
# guiradius = 0.04
# radius = 600*guiradius # More accurate visually

def wallcollision(centres, velocities):
    for i in range(len(centres)):
        if centres[i][0] - guiradius <=0:
            velocities[i][0] = -damping*velocities[i][0]
            centres[i][0] = guiradius
        if centres[i][0] + guiradius >=1:
            velocities[i][0] = -damping*velocities[i][0]
            centres[i][0] = 1 - guiradius
        if centres[i][1] - guiradius <=0:
            velocities[i][1] = -damping*velocities[i][1]
            centres[i][1] = guiradius
        if centres[i][1] + guiradius >=1:
            velocities[i][1] = -damping*velocities[i][1]
            centres[i][1] = 1 - guiradius
    return velocities

def selfcollision(centres, velocities):
    for i in range(len(centres)):
        for j in range(i+1,len(centres)):
            distance = np.linalg.norm(centres[i]-centres[j])
            if distance < guiradius * 2:
                print("Collision!")
                vitemp = np.copy(velocities[i])
                vitemp = velocities[i] - (np.dot(velocities[i]-velocities[j], centres[i]-centres[j])*(centres[i]-centres[j])/distance**2)
                velocities[j] = velocities[j] - (np.dot(velocities[j]-velocities[i], centres[j]-centres[i])*(centres[j]-centres[i])/distance**2)
                velocities[i] = vitemp

                # Fix the centers
                overlapping_distance = guiradius * 2 - distance
                direction = (centres[j] - centres[i]) / distance
                epsilon = 1e-6 # Additional safe gap to avoid particles to collide again
                shift = overlapping_distance * direction * (0.5 + epsilon)
                centres[i] -= shift
                centres[j] += shift
    return velocities

Please note that in more realistic simulations, the force is dependent of the collision surface (using a non-trivial formula if one want to simulate electrons responsible for the ball repulsion) so the collisions should not occur for a long time.

Answered By - Jérôme Richard