Particles sticking together after collision...sometimes

[Joseph Meyer]

Hello, I have viewed the sample gas.py code that comes with Vpython and am trying to replicate "particles in a box colliding with one another" as my foray into Vpython. However, my code is telling some of the particles to stick together after colliding. I cannot figure out why. I think I am missing some basic concept rather than there being an issue with my code. Below is my code, anyhow. Looking forward to any advice!

from visual import *

import numpy

frame = box(pos=(0,0,0), size=(10,10,1), color=color.cyan, opacity=.5)

rad=.25

numParticles=10

deltat=.005

maxSpeed=10

minSpeed=0

particles=[]

vel=[]

pos=[]

i=0

count=0

#randomly choose starting locations of particles

while i < numParticles:

    coord=True

    

    tempX = random.uniform(-frame.size.x/2. + rad, frame.size.x/2. - rad)

    tempY = random.uniform(-frame.size.y/2. + rad, frame.size.y/2. - rad)

    tempZ = 0#random.uniform(-frame.size.z/2. + rad, frame.size.z/2. - rad)

    for j in range(len(particles)):

        if pow(tempX-particles[j].pos.x, 2) + pow(tempY-particles[j].pos.y, 2) + pow(tempZ-particles[j].pos.z, 2) < pow(2*rad, 2):

            coord=False

            break

        

    if coord==True:

        i=i+1

        particles.append(sphere(pos=(tempX,tempY,tempZ), radius=rad, color=color.magenta))

        pos.append((tempX, tempY, tempZ))

    count=count+1

    if count==1000:

        print (str(i) + " balls produced. No room for any more.")

        break

    

#randomly choose starting velocities of particles

for i in range(len(particles)):

    velX=random.uniform(minSpeed, maxSpeed)

    velY=random.uniform(minSpeed, maxSpeed)

    velZ=0#random.uniform(minSpeed, maxSpeed)

    

    particles[i].velocity=vector(velX, velY, velZ)

    vel.append((velX, velY, velZ))

    

posArray=numpy.array(pos)

velArray=numpy.array(vel)

while True:

    rate (100)

    #update position of every particle by v*dt

    posArray=posArray+(velArray*deltat)

    

    #determine if a particle collides with another particle

    diff=posArray.reshape(numParticles,1,3)-posArray #difference between every particle and every other particle

    D = (diff**2).sum(2) #square the distances and sum xyz dimensions for each particle

    index=numpy.arange(numParticles) 

    D[index, index]=numpy.inf #set the diagonal = infinity to remove distance between each particle and itself (obviously zero and not a collision)

    collision=numpy.less_equal(D, (2*rad)**2) #collison is a boolean array of True if a particle is colliding with another particle and False if it is not

    collisionList = numpy.transpose(numpy.nonzero(collision)) #collisionList is a list of all indices of all particles colliding

    

    duplicate=[]

    for i in range(len(collisionList)):

        j=collisionList[i][0]

        k=collisionList[i][1]

        

        if [j, k] not in duplicate:

            duplicate.append([k, j])

            

            finalVx1=0

            finalVy1=0

            finalVz1=0

            

            finalVx2=0

            finalVy2=0

            finalVz2=0

            

            Vx1=velArray[j][0]

            Vy1=velArray[j][1]

            Vz1=velArray[j][2]

            Vx2=velArray[k][0]

            Vy2=velArray[k][1]

            Vz2=velArray[k][2]

            xcoeff=[1, -1*(Vx1+Vx2), Vx1*Vx2]

            ycoeff=[1, -1*(Vy1+Vy2), Vy1*Vy2]

            zcoeff=[1, -1*(Vz1+Vz2), Vz1*Vz2]

            xroots=numpy.roots(xcoeff)

            yroots=numpy.roots(ycoeff)

            zroots=numpy.roots(zcoeff)

            if xroots[0]!=Vx2:

                finalVx2=xroots[0]

            else:

                finalVx2=xroots[1]

            if yroots[0]!=Vy2:

                finalVy2=yroots[0]

            else:

                finalVy2=yroots[1]

            if zroots[0]!=Vz2:

                finalVz2=zroots[0]

            else:

                finalVz2=zroots[1]

                

            finalVx1 = Vx1 + Vx2 - finalVx2

            finalVy1 = Vy1 + Vy2 - finalVy2

            finalVz1 = Vz1 + Vz2 - finalVz2

            velArray[j] = [finalVx1, finalVy1, 0]#finalVz1]

            velArray[k] = [finalVx2, finalVy2, 0]#finalVz2]

                                       

    #determine if a particle hits a wall

    wallHit=numpy.less_equal(frame.size.y/2., rad+abs(posArray)) #wallHit is a boolean array of True if a particle's position is outside frame and False if inside frame

    V=velArray*wallHit #V is velArray with only velocities of particles that are outside the frame

    V=V*-1 #this reverses all the velocities to make them "bounce" off the wall

    velArray=velArray*numpy.less_equal(rad+abs(posArray), frame.size.y/2.) #this makes velArray only have values for particles within the frame. Particles outside the frame get vel of zero

    velArray=velArray+V #add the arrays to update velArray

    for i in range(numParticles):

        particles[i].pos = posArray[i]

        particles[i].velocity=vector(velArray[i])

[Bruce Sherwood]

You have omitted some crucial code that is in the gas example, having to do with splitting the time interval into a portion before impact and a portion after. As a result, when two slow particles collide, they can be trapped inside each other, and as long as their speeds are low, they keep being identified as closer than a diameter. Look closely at the animation and you'll see that the "molecules" fibrillate, never moving far enough away from each other to get free.

[Aaron Titus]

Bruce gave you the best answer, but you may be interested in knowing how I teach collisions to non-science majors in a general education course called “Physics for Video Games”. 

One of our games is Pong. When a ball and wall overlap (i.e. a collision), we step backwards using ball.pos = ball.pos - v*dt. Then we compute the new velocity of the ball after the collision. This ensures that in the next step, the ball and wall will no longer overlap. If we do not do this, then we have the same problem that you noticed. A ball can become stuck within a wall and never escape from the collision.

Aaron

--
You received this message because you are subscribed to the Google Groups "VPython-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to vpython-user...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Joseph Meyer]

Thank you Aaron and Bruce. It was quite a simple change and now I've got it working properly. Basically, you never really let the particles "touch" each other or the walls. Thanks again for the quick clarification and advice!

[Joseph Meyer]

I've noticed a new undesired effect with particles sticking together. If I have a small number of particles and a small number of collisions per deltat, everything seems to work fine. However, when I add more particles and it gets more crowded, particle end up sticking together and causing issues again. I think this might be because as it gets more crowded, particles "collide" with more than just one other particle within a single deltat. Could this be the reason? I haven't tried to write the code to account for this scenario yet, I was just hoping to confirm if I'm on the right track here or if there are any other suggestions out there.

[Bruce Sherwood]

I'm guessing that there's something still missing from the algorithm, because I increased the number of objects from 100 to 300 in gas.py and I don't think that I see any sticking.

[Joseph Meyer]

I've managed to eliminate the particles sticking together by paying closer attention to the gas.py example. Now, my code is very similar to this example. However, now some of the particles get stuck on the walls of the box. I noticed that in gas.py, there is no step back of the particle's position to the exact time of collision with the wall like there is with the collision of 2 particles. Why is that? 

[Bruce Sherwood]

The difference is that in the collision with the wall the magnitude of the momentum doesn't change, and in particular doesn't get small, so the atom can't get trapped, unlike the case of colliding with another atom.

[Joseph Meyer]

I must still be missing something. I can't seem to find anything missing in my code but after a minute or so, a particle will get stuck in the wall. If anyone can offer input on my code below that would be very appreciated. I'll keep thinking though...

from visual import *
import numpy as np

L=2
frame = box(pos=(0,0,0), size=(L,L,L), color=color.cyan, opacity=.5)

rad=.4
numParticles=4
deltat=.0005
maxSpeed=15
minSpeed=10

particles=[]
vel=[]
pos=[]

#randomly choose starting locations of particles

i=0
count=0
while i < numParticles:
    coord=True

    tempX = random.uniform(-L/2 + rad, L/2 - rad)
    tempY = random.uniform(-L/2 + rad, L/2 - rad)
    tempZ = random.uniform(-L/2 + rad, L/2 - rad)
       
    for j in range(len(particles)):
        if pow(tempX-particles[j].pos.x, 2) + pow(tempY-particles[j].pos.y, 2) + pow(tempZ-particles[j].pos.z, 2) <= pow(2*rad, 2):

            coord=False
            break
       
    if coord==True:
        i=i+1
        particles.append(sphere(pos=(tempX,tempY,tempZ), radius=rad, color=color.magenta))
        pos.append((tempX, tempY, tempZ))

    count=count+1
    if count==1000:
        print (str(i) + " balls produced. No room for any more.")
        break
   
#randomly choose starting velocities of particles
for i in range(len(particles)):
   
    velX=random.uniform(minSpeed, maxSpeed)
    velY=random.uniform(minSpeed, maxSpeed)

    velZ=random.uniform(minSpeed, maxSpeed)

   
    particles[i].velocity=vector(velX, velY, velZ)
    vel.append((velX, velY, velZ))
   

posArray=np.array(pos)
velArray=np.array(vel)


while True:
    rate (250)

   
    #update position of every particle by v*dt
    posArray=posArray+(velArray*deltat)

    #determine if a particle collides with another particle

    diff=posArray.reshape(numParticles,1,3)-posArray #difference in distance between every particle and every other particle
    index=np.arange(numParticles)
    D = (diff**2).sum(2) #square the difference in xyz dimensions between each particle and sum xyz dimensions for each particle
    D[index, index]=np.inf #set the diagonal = infinity to remove distance between each particle and itself (obviously zero and not a collision)
    collision=np.less_equal(D, (2*rad)**2) #collison is a boolean array of True if a particle is colliding with another particle and False if it is not
    collisionList = np.transpose(np.nonzero(collision)).tolist()
   
    for i,j in collisionList:
        collisionList.remove([j,i])
        Vi=velArray[i]
        Vj=velArray[j]

        #calculate exact time particles collided
        a = mag(Vj-Vi)**2
        if a == 0: print("same velocities!") # exactly same velocities
        b = 2*dot(diff[i][j],Vj-Vi)
        c = mag(diff[i][j])**2-(2*rad)**2
        d = b**2-4.*a*c
        if d < 0: print ("something went wrong") # something wrong; ignore this rare case
        t = (-b+sqrt(d))/(2.*a) # deltat-t is when they made contact
   
        posArray[i]=posArray[i]-(Vi*t) #back up the particles to the exact point collision
        posArray[j]=posArray[j]-(Vj*t) # back up the particles to the exact point of collision

        #calculate new velocities
        V1f= Vi - (np.dot(Vi-Vj, posArray[i]-posArray[j])/pow(mag(posArray[i]-posArray[j]),2)) * (posArray[i]-posArray[j])
        V2f= Vj - (np.dot(Vj-Vi, posArray[j]-posArray[i])/pow(mag(posArray[j]-posArray[i]),2)) * (posArray[j]-posArray[i])

        velArray[i]=V1f
        velArray[j]=V2f
        posArray[i]=posArray[i]+(velArray[i]*t) #step forward the particles
        posArray[j]=posArray[j]+(velArray[j]*t) #step forward the particles
           
    #determine if a particle hits a wall and bounce it off
    outsideBool=np.less_equal(L/2., rad+abs(posArray))
    insideBool=np.less(rad+abs(posArray), L/2.)
    outside=velArray*outsideBool #array of only vel values that are outside the frame
    inside=velArray*insideBool
    outside=-1*outside
    velArray=inside + outside

    for i in range(numParticles):
        particles[i].pos = posArray[i]

[Bruce Sherwood]

I don't know what the problem is, but maybe the behavior that differs from the example program gas.py is due to some difference in the parameters, not a difference in the algorithm. Two suggestions:

1) Do for the walls what you do for the particles, including moving the particle inward from the wall at the same time you force the momentum inward.

2) Keep a copy of key data from the preceding time step, find a way to identify trapping at a wall, and in that case print the preceding and current data.

[Joseph Meyer]

I did #1 per your suggestion below and I am now going to consider this solved. I wish I understood which parameters that are different in my code vs. the example code that required me to do this. Nevertheless, it works and I am satisfied. Thank you for all the support. Here's the addition to the wall collision code I made:

 #determine if a particle hits a wall and bounce it off

    outsideBool=np.less_equal(L/2, rad+abs(posArray))
    insideBool=np.less(rad+abs(posArray), L/2)
    outsideVel=velArray*outsideBool #array of only vel values that are outside the frame
    insideVel=velArray*insideBool
    wallList=np.transpose(np.nonzero(outsideBool)).tolist()

    #update velocities in tempVel
    outsideVel=-1*outsideVel
    tempVel=insideVel+outsideVel
   
    for i,j in wallList:
        wallList.remove([i,j])
        #exact time particle hit the wall
        dd=(abs(posArray[i][j])+rad)-(L/2) #distance beyond the wall the particle has gone
        dt=dd/abs(velArray[i][j]) #time in which it took the particle to go that distance (deltat-dt is the exact time of the collision)

        #back up the particle
        posArray[i]=posArray[i]-(velArray[i]*dt)

        #new velocity is already updated in tempVel
        #step forward the particle with new velocities
        posArray[i]=posArray[i]+(tempVel[i]*dt)

    #update the velArray
    velArray=tempVel