重力粒子模拟：粒子静止处理及碰撞穿透/粘连问题求助

问题核心分析 #

你遇到的两个现象本质是同一问题：

• 重力场中小球末期穿透地面而非静止
• 无重力下多粒子碰撞粘连

根源在于离散时间步长的局限性：当粒子速度过小时，碰撞后速度不足以让粒子在下一帧脱离边界，导致反复触发碰撞、速度翻转，最终出现穿透/嵌入的异常状态。当前的handleBoxCollision()仅修改速度，未修正粒子位置，当粒子已经嵌入边界内部时，单纯反转速度无法让它回到边界外，陷入死循环。

解决方案 #

需要从两个关键维度修正：

1. 碰撞时强制修正位置：将嵌入边界的粒子推回至边界外部，避免重复触发碰撞
2. 添加静止阈值逻辑：当粒子速度足够小时，直接将其速度设为0并关闭重力，模拟现实中静止状态

修正后的碰撞处理函数 #

def handleBoxCollision():
    cor = 0.8
    # 处理左右边界碰撞
    if p.left <= box.left:
        # 计算重叠量并将粒子推回边界外
        overlap = box.left - p.left
        p.position[0] += overlap
        p.velocity[0] = -cor * p.velocity[0]
        # 速度过小则直接静止
        if abs(p.velocity[0]) < 0.1:
            p.velocity[0] = 0
    elif p.right >= box.right:
        overlap = p.right - box.right
        p.position[0] -= overlap
        p.velocity[0] = -cor * p.velocity[0]
        if abs(p.velocity[0]) < 0.1:
            p.velocity[0] = 0

    # 处理上下边界碰撞（重点优化地面静止逻辑）
    if p.bottom <= box.bottom:
        overlap = box.bottom - p.bottom
        p.position[1] += overlap
        # 反弹速度低于阈值时，直接静止并关闭重力
        if abs(p.velocity[1]) < 0.2:
            p.velocity[1] = 0
            p.acceleration[1] = 0
        else:
            p.velocity[1] = -cor * p.velocity[1]
    elif p.top >= box.top:
        overlap = p.top - box.top
        p.position[1] -= overlap
        p.velocity[1] = -cor * p.velocity[1]
        if abs(p.velocity[1]) < 0.1:
            p.velocity[1] = 0

    # 更新粒子边界坐标
    p.left = p.position[0] - p.radius
    p.right = p.position[0] + p.radius
    p.top = p.position[1] + p.radius
    p.bottom = p.position[1] - p.radius

额外优化点 #

• 格式化输出的时间和速度数值，提升可读性
• 多粒子碰撞场景可复用相同逻辑：计算粒子间重叠量、推回位置、添加速度阈值避免粘连

修正后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt
import math
import os

class Particle:
    def __init__(self, mass, position, velocity, acceleration):
        self.mass = mass
        self.radius = math.sqrt(self.mass/(math.pi*1.5))
        self.position = position
        self.velocity = velocity
        self.acceleration = acceleration
        self.KE = (1/2)*self.mass*np.dot(self.velocity,self.velocity)
        
        self.left = position[0]-self.radius
        self.right = position[0]+self.radius
        self.top = position[1]+self.radius
        self.bottom = position[1]-self.radius
        
class Box:
    def __init__(self):
        self.left = -10
        self.right = 10
        self.bottom = -10
        self.top = 10
        
def handleBoxCollision():
    cor = 0.8
    # 左右边界碰撞处理
    if p.left <= box.left:
        overlap = box.left - p.left
        p.position[0] += overlap
        p.velocity[0] = -cor * p.velocity[0]
        if abs(p.velocity[0]) < 0.1:
            p.velocity[0] = 0
    elif p.right >= box.right:
        overlap = p.right - box.right
        p.position[0] -= overlap
        p.velocity[0] = -cor * p.velocity[0]
        if abs(p.velocity[0]) < 0.1:
            p.velocity[0] = 0

    # 上下边界碰撞处理
    if p.bottom <= box.bottom:
        overlap = box.bottom - p.bottom
        p.position[1] += overlap
        if abs(p.velocity[1]) < 0.2:
            p.velocity[1] = 0
            p.acceleration[1] = 0
        else:
            p.velocity[1] = -cor * p.velocity[1]
    elif p.top >= box.top:
        overlap = p.top - box.top
        p.position[1] -= overlap
        p.velocity[1] = -cor * p.velocity[1]
        if abs(p.velocity[1]) < 0.1:
            p.velocity[1] = 0

    # 更新粒子边界坐标
    p.left = p.position[0] - p.radius
    p.right = p.position[0] + p.radius
    p.top = p.position[1] + p.radius
    p.bottom = p.position[1] - p.radius
        
def update(dt):
    p.velocity = p.velocity + (p.acceleration * dt)
    p.position = p.position + (p.velocity * dt)
    p.left = p.position[0]-p.radius
    p.right = p.position[0]+p.radius
    p.top = p.position[1]+p.radius
    p.bottom = p.position[1]-p.radius
    
    handleBoxCollision()

# 清空旧帧文件
dir = 'C:/PathToImageFile'
for f in os.listdir(dir):
    os.remove(os.path.join(dir, f))

# 初始化模拟参数
mass = 10
position = np.array([0,0])
velocity = np.array([5,0])
acceleration = np.array([0,-9.8])
dt = 1/60

p = Particle(mass, position, velocity, acceleration)
box = Box()

# 运行模拟并生成帧
for i in range(600):
    figure, axes = plt.subplots()
    
    update(dt)
    
    # 绘制盒子边界
    plt.plot([-10,-10],[-10,10],color='black')
    plt.plot([10,10],[-10,10],color='black')
    plt.plot([-10,10],[10,10],color='black')
    plt.plot([-10,10],[-10,-10],color='black')
    
    # 绘制粒子
    cc = plt.Circle((p.position[0] ,p.position[1]), p.radius)
    plt.scatter(p.position[0],p.position[1])
    plt.xlim(-11,11)
    plt.ylim(-11,11)
    # 格式化显示时间和速度
    plt.text(-10,-10.7,f"time={i*dt:.2f}")
    plt.text(-10,10.3,f"velocity={p.velocity[1]:.2f}")
    axes=plt.gca()
    axes.add_artist(cc)
    axes.set_aspect(1)
    
    figure.savefig(f'/PathToImageFile{i}.png')
    plt.close('all')

关键逻辑说明 #

• 位置修正：通过计算粒子与边界的重叠量，强制将粒子推回至边界外部，从根源避免重复触发碰撞
• 静止阈值：设置合理的速度阈值（比如0.2），当反弹速度低于该值时，直接将速度设为0并关闭重力，完美模拟粒子静止在地面的物理状态
• 多粒子扩展：粒子间碰撞可复用相同逻辑，计算粒子间重叠量、推回位置、添加速度阈值即可解决粘连问题

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