vpython数值模拟系列(1.小球间的弹性碰撞)

        笔者学习VPython的时候曾遇到了不少困难，所以在此撰写系列教程与案例，希望能对学习VPython提供相关借鉴。如有任何错误或疑问，欢迎评论指出！

        本案例将实现立体空间内两个小球的弹性碰撞，以下为详细教程：

1.安装VPython库
安装Python第三方库有多种方式，详细可参考其他教程，本教程不过多赘述。

2.导入VPython库

from vpython import*

3.初始变量的设定
本案例采用微元法以实现分子运动，时间微元设为，初始时间为0，小球半径为1。

R, t, dt = 1, 0, 1E-5

4.画布与小球的构建

通过以下代码将实现画布与小球的构建。其中画布的颜色定义为白色，宽度为500，高度为500。两个小球，一个位于（0，5，0）位置，颜色为青色，另一个小球则位于（5，0，5）位置，颜色为绿色，并使两小球均保留运动轨迹。（注：这里的位置使用3维向量vector类型传递参数，默认两小球质量相等）

animation = canvas(width=500, height=500, background=color.white)
ball_1 = sphere(pos=vector(0, 5, 0), radius=R, color=color.cyan, make_trail=True)
ball_2 = sphere(pos=vector(5, 0, 5), radius=R, color=color.green, make_trail=True)

5.速度与碰撞函数

为确保小球能发生碰撞，速度选取特定值。发生弹性碰撞时，速度变换的公式可写为：

其具体推导过程可参阅相关书籍[1]。根据以上便可编写代码，定义小球碰撞的函数：

v1, v2 = -vector(0, 9, 0), -vector(10, 0, 10)
def collision(pos1, pos2, v1, v2):
    ez = (ball_1.pos - ball_2.pos).norm()
    u1 = v1 - dot(v1 - v2, ez) * ez
    u2 = v2 + dot(v1 - v2, ez) * ez
    return u1, u2

6.不断小球更新位置

使用wihle循环不断更新小球位置。如果小球达到了可以发生碰撞的条件，则执行碰撞函数。

while t <= 1:
    rate(1/dt)
    t += dt
    ball_1.pos += v1 * dt
    ball_2.pos += v2 * dt
    if mag(ball_1.pos - ball_2.pos) <= 2*R and dot(ball_1.pos - ball_2.pos, v1 - v2) <= 0:
        v1, v2 = collision(ball_1.pos, ball_2.pos, v1, v2)

---

附录（相关物理推导[1]）：

1.速度变换公式推导：

小球的运动遵循牛顿运动规律，系统的总动量守恒，能量守恒。考虑两个等质量，等直径的小球在一维情况下的碰撞。

在时刻速度分别是和 ，考虑它们在时进行碰撞，碰撞之后速度为和。

由动量守恒：

和能量守恒，

可以得到碰撞之后的速度，

也就是说在一维运动的情况下，等质量的小球碰撞后只是交换其速度。现在考虑二维或三维的情况，当球体发生碰撞时，排斥力是沿球心连线方向，只有平行于方向的速度矢量分量改变。取该方向为一坐标系的 z 轴，即。由一维情形下的结果可知，两小球的平行于该方向的速度分量在碰撞后交换，即

而其他分量保持不变，即

这两个方程是镜面反射条件。由此可得，

2.碰撞的判据：

只有碰撞时才有相互作用，其余时间中小球是以一定速度沿直线运动的，因此模拟的直接目标不是运动轨迹，而是碰撞时间。如果已知两个相同大小的硬球在时刻的初始条件，现在我们来求它们的碰撞时间。两小球发生碰撞时有，

 （1）

由于小球是以直线运动的，有

则碰撞条件为

其中

则得到所求的时间

  （2）

上式的号对应两个小球接触时有两个几何位置，由于小球不能互相穿透，只有对应于负号的短时间解才是物理解，又大于，所以不难得到

两个小球的相对速度在其中心连线方向上的投影分量小于零，说明它们之间是相对趋近的，才有可能碰撞。另外，发生碰撞时满足(1)式，这时(2)式中的根式部分中的内容恒大于0，由此我们便得到了碰撞发生的条件。由于浮点误差，所以(1)式难取等号，在模拟中的碰撞发生条件改写为下式，

---

参考文献
[1]. Microsoft Word - 4－1.doc (ustc.edu.cn)