粒子碰撞模拟异常求助：球体粘连、边缘卡顿问题排查

问题分析 #

当前实现存在两个核心异常：

• 边界卡滞：仅反转球体速度但未修正穿透边界的位置，导致球体持续触发边界碰撞检测，反复反转速度，表现为卡在边界无法移动。
• 球体粘连旋转：球体碰撞时已发生位置穿透，未修正重叠状态；同时碰撞速度计算逻辑存在误差，导致球体持续处于碰撞判定状态，反复更新速度出现粘连旋转现象。

修复方案 #

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  # 修正位置到边界内侧
        elif 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
        elif centres[i][1] + guiradius >= 1:
            velocities[i][1] = -damping * velocities[i][1]
            centres[i][1] = 1 - guiradius
    return velocities

2. 球体间碰撞修复 #

• 先修正球体重叠位置：将两个球体沿碰撞法线方向推开，消除穿透状态
• 简化并修正弹性碰撞速度计算逻辑（假设所有球体质量相同）：

def selfcollision(centres, velocities):
    for i in range(len(centres)):
        for j in range(i+1, len(centres)):
            delta_pos = centres[i] - centres[j]
            dist = np.linalg.norm(delta_pos)
            # 避免除以0，仅处理真实碰撞
            if dist < 2 * guiradius and dist > 1e-8:
                # 修正重叠位置，将两个球沿法线方向分开
                overlap = 2 * guiradius - dist
                norm = delta_pos / dist
                centres[i] += norm * overlap / 2
                centres[j] -= norm * overlap / 2
                
                # 计算弹性碰撞后的速度
                vel_diff = velocities[i] - velocities[j]
                dot_product = np.dot(vel_diff, norm)
                # 仅处理相向运动的碰撞，避免反向触发
                if dot_product > 0:
                    velocities[i] -= norm * dot_product
                    velocities[j] += norm * dot_product
    return velocities

3. 更新逻辑调整 #

调整处理顺序，先修正边界位置再处理球体碰撞，避免位置修正后再次触发边界问题：

def update(centres, velocities):
    centres += velocities * timestep
    wallcollision(centres, velocities)  # 优先处理边界，修正位置
    selfcollision(centres, velocities)  # 再处理球体间碰撞

完整修正代码 #

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

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
        elif 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
        elif 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)):
            delta_pos = centres[i] - centres[j]
            dist = np.linalg.norm(delta_pos)
            if dist < 2 * guiradius and dist > 1e-8:
                overlap = 2 * guiradius - dist
                norm = delta_pos / dist
                centres[i] += norm * overlap / 2
                centres[j] -= norm * overlap / 2
                
                vel_diff = velocities[i] - velocities[j]
                dot_product = np.dot(vel_diff, norm)
                if dot_product > 0:
                    velocities[i] -= norm * dot_product
                    velocities[j] += norm * dot_product
    return velocities

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

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

后续优化建议 #

后续切换到Taichi装饰器时，建议将centres和velocities改为Taichi场（ti.Vector.field），用Taichi内核处理碰撞和更新逻辑，彻底规避Python循环的性能瓶颈，同时进一步优化碰撞检测精度。

内容的提问来源于stack exchange，提问作者ollie d