Search code examples
pythonnumpyphysics

What on earth is going in in my particle simulation?


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.