粒子碰撞模拟Python

问：

我正在尝试创建一个相对简单的粒子模拟，它应该考虑重力、阻力、与其他粒子的碰撞（非弹性碰撞）以及与墙壁的碰撞（完全弹性）。我得到了使用速度 Verlet 算法的重力和阻力部分，但截至目前，它无法将粒子设置为平衡状态。此外，如果我添加多个粒子，它们有时会相互攀爬，这是由于（我相信）它们仍然具有非常小的速度分量，渐近地驱动到零。如果粒子的能量变得足够小，我试图切断粒子的速度，但这看起来并不现实。有人可以指出一些如何解决这些问题的建议。我得到了一个粒子对象：

import pygame
import random
import numpy as np
import operator
from itertools import combinations

class Particle:
    def __init__(self):
        self.mass = 10
        self.radius = random.randint(10, 50) 
        self.width, self.height = 700, 500
        self.pos = np.array((self.width/2, self.height/2))
        self.v = np.array((0.0, 0.0)) 
        self.acc = np.array((0.0, 0.0))
        self.bounce = 0.95

我使用 Verlet-Integration 来考虑重力和阻力：

def update(self, ball, dt):
        new_pos = np.array((ball.pos[0] + ball.v[0]*dt + ball.acc[0]*(dt*dt*0.5), ball.pos[1] + ball.v[1]*dt + ball.acc[1]*(dt*dt*0.5)))
        new_acc = np.array((self.apply_forces(ball))) # only needed if acceleration is not constant
        new_v = np.array((ball.v[0] + (ball.acc[0]+new_acc[0])*(dt*0.5), ball.v[1] + (ball.acc[1]+new_acc[1])*(dt*0.5)))
        ball.pos = new_pos;
        ball.v = new_v;
        ball.acc = new_acc;
   
    def apply_forces(self, ball):
        grav_acc = [0.0, 9.81]
        drag_force = [0.5 * self.drag * (ball.v[0] * abs(ball.v[0])), 0.5 * self.drag * (ball.v[1] * abs(ball.v[1]))]  #D = 0.5 * (rho * C * Area * vel^2)
        drag_acc = [drag_force[0] / ball.mass, drag_force[1] / ball.mass] # a = F/m

        return (-drag_acc[0]),(grav_acc[1] - drag_acc[1])

在这里，我计算碰撞部分：

    def collision(self):
        pairs = combinations(range(len(self.ball_list)), 2)

        for i,j in pairs:
            part1 = self.ball_list[i]
            part2 = self.ball_list[j]
            distance = list(map(operator.sub, self.ball_list[i].pos, self.ball_list[j].pos))

            if np.hypot(*distance) < self.ball_list[i].radius + self.ball_list[j].radius:
                distance = part1.pos - part2.pos
                rad = part1.radius + part2.radius
                slength = (part1.pos[0] - part2.pos[0])**2 + (part1.pos[1] - part2.pos[1])**2
                length = np.hypot(*distance)
                factor = (length-rad)/length;
                x = part1.pos[0] - part2.pos[0]
                y = part1.pos[1] - part2.pos[1]

                part1.pos[0] -= x*factor*0.5
                part1.pos[1] -= y*factor*0.5
                part2.pos[0] += x*factor*0.5
                part2.pos[1] += y*factor*0.5

                u1 = (part1.bounce*(x*part1.v[0]+y*part1.v[1]))/slength
                u2 = (part2.bounce*(x*part2.v[0]+y*part2.v[1]))/slength
                part1.v[0] = u2*x-u1*x
                part1.v[1] = u1*x-u2*x
                part2.v[0] = u2*y-u1*y
                part2.v[1] = u1*y-u2*y

    def check_boundaries(self, ball):
        if ball.pos[0] + ball.radius > self.width:
            ball.v[0] *= -ball.bounce
            ball.pos[0] = self.width - ball.radius
                
        if ball.pos[0] < ball.radius:
            ball.v[0] *= -ball.bounce
            ball.pos[0] = ball.radius
                
        if ball.pos[1] + ball.radius > self.height:
            self.friction = True
            ball.v[1] *= -ball.bounce
            ball.pos[1] = self.height - ball.radius

        elif ball.pos[1] < ball.radius:
            ball.v[1] *= -ball.bounce
            ball.pos[1] = ball.radius