我的粒子模拟到底发生了什么？

Question

我的碰撞球模拟已经可以运行，但粒子变得疯狂。

其中一些相互粘在一起并旋转。其他人则陷入边缘。

我想不出为什么......有人可以帮助我吗？谢谢

此代码使用 Taichi Gui。我打算将来使用 taichi 装饰器，但现在我只是用常规 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)

Answer 1

问题是，当检测到碰撞时，速度会被修改，但中心不会被修改，并且 速度并不总是足够大，以便粒子在下一帧中逃脱碰撞。这意味着粒子可以在许多帧中保持碰撞。这会导致速度分量在负值和正值之间振荡，从而有时无法避免碰撞。这就是为什么颗粒会粘在一起。同样的问题也发生在墙壁上：当粒子中心超出墙壁时，速度分量可能会振荡，因此粒子会粘在墙壁上。此类问题在简单的物理粒子模拟中经常出现。一个简单的修复方法就是校正中心，这样粒子就不会碰撞下一帧。

这是更正后的代码：

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))
guiradius = 0.04
radius = 600*guiradius # More accurate visually
damping = 0.9
timestep = 1/60

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]
            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
                if distance < guiradius * 2:
                    overlapping_distance = guiradius * 2 - distance
                    direction = (centres[j] - centres[i]) / distance
                    epsilon = 1e-6 # Additional safe gap to avoid particles to colide again
                    shift = overlapping_distance * direction * (0.5 + epsilon)
                    centres[i] -= shift
                    centres[j] += shift
    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)

请注意，在更现实的模拟中，力取决于碰撞表面（如果想模拟负责球排斥的电子，请使用重要的公式），因此碰撞不应长时间发生。