大量小球无规则运动计算机模拟的可视化

邓 晓 胡其图 李发泽 张小灵

（1上海交通大学物理系，上海 200240）

（2宁夏大学物理电气信息学院，宁夏 银川 750021）

理想气体分子热运动、布朗运动、气体分子扩散等物理现象的模拟，都可以抽象为大量小球无规则运动的形式.使用计算机对这种运动形式进行实时模拟有着广泛的应用.由于系统中小球数目众多，小球之间存在大量的弹性碰撞，实现实时模拟可视化的关键在于如何减少计算量、提高运行效率.

一种常用的办法是：系统按固定时间步长向前运动，在系统准备向前走一步时，先对所有的小球扫描一遍，检测它们之间是否将发生重叠.如果将发生重叠，则认为会有碰撞发生，按照弹性碰撞所满足的规律改变相关小球的运动速度.这种算法的每次扫描过程需要O（n2）的运算量，在小球数目较多时运算量很大.另外，由于所取的时间步长不能无限小，当小球运动速度较快或小球半径相对较小时，会使得本来应该检测到的碰撞没有被检测到，从而在实现实时模拟可视化时容易出现两个小球相互重叠或穿透的情形.本方法由于运行效率不高，模拟的小球数目到达600个时系统就很难运行了，并且小球之间的重叠和穿透现象无法避免［1］.

Dong－Jin Kim 提出一种基于碰撞事件驱动的模拟方法［2］，我们利用此方法实现了大量小球无重叠、无穿透运动的计算机模拟，系统不再按固定时间步长运行，而是按实际碰撞发生的时间来运行.这种方法有效地解决了上述问题，并大大提高了计算效率.考虑一个拥有大量小球的二维封闭系统，其边界是可移动的（例如：汽缸的活塞），因而也具有一定的运动速度，下面我们来具体讨论本模拟方法.

1 单次碰撞的碰撞检测与碰撞反应

系统中包含了两种类型的碰撞：小球与容器壁的碰撞、小球与小球之间的碰撞.下面分别讨论它们的碰撞检测与碰撞反应.

1.1 小球与容器壁的碰撞

一个容器在二维空间中可用一个长方形来表示，容器壁有四个：左右两个竖直壁和上下两个水平壁.小球与水平壁之间的碰撞和小球与竖直壁之间的碰撞算法原理相同，本文只给出小球与竖直左壁在一个时间步长T 内发生碰撞的情况.

如图1所示，左容器壁处于xw位置，运动速度为vw.小球初始位置为（x，y），碰撞前在x 方向上的运动速度为vx.两者经过时间tc后相撞，此时小球球心处于xc位置，而容器壁也将达到粗虚线位置，与小球相切.碰撞后小球在x 方向上的速度为v′x，小球在y 方向上的速度不变.再经过时间tr＝T－tc，小球将到达x′位置.

图1 小球与容器壁的碰撞

T 时间后小球和器壁的位置如下：

其中，xc=x+vxtc.

显然，我们需要先得到碰撞后小球在x 方向的速度v′x和碰撞时的时间tc.

（1）碰撞后小球在x 方向的速度v′x

碰撞前在x 方向上小球相对容器壁的速度为

由于碰撞后容器壁速度不变，所以碰撞后在x 方向上小球相对器壁的速度为

v′xr＝v′x－vw

显然ｖ′ｘｒ＝－ｖｘｒ

所以ｖ′ｘ＝２ｖｗ－ｖｘ

（2）碰撞时的时间tc

设定小球球心可以到达的边界位置为xb，则

R 的绝对值为小球半径.当R 为正数时，为左边容器壁；当R 为负数时，为右边容器壁.

以容器壁为参考系，小球相对容器壁的碰撞过程如图2.容器壁静止，小球在x 方向上以相对速度vxr靠近容器壁.

图2 以容器壁为参考系的碰撞示意图

dx表示碰撞前的相对位移，d′x为碰撞后的位移，显然

取负号时表示沿x 的负方向.

在时间步长T 中，小球与器壁发生碰撞的条件为：

①当容器壁为左边容器壁，即dx＞0 时，dx≤－vxrT；

②当容器壁为右边容器壁，即dx＜0 时，－dx≤vxrT.

如果发生碰撞，则发生碰撞的时间

联合式（1）、（2）、（3），得

1.2 小球与小球之间的碰撞

假设有两个半径不同、质量不同的小球相撞，如图3展示了在T 时间步长里小球A 与小球B发生碰撞的情况.

图3 小球与小球之间的碰撞

对于小球A，它从P1位置以vA的速度经过时间tc后运动到P2的位置并与小球B发生碰撞，速度改变为v′A，然后再以碰撞后的速度继续运动tr的时间到达P3位置，其中tr＝T－tc.小球B 的运动也是如此.如果我们能够计算出碰撞的发生时间tc以及碰撞后小球A、B的速度，那么整个过程就被求解了.

（1）小球碰撞后的速度

图4展示了小球A 在位置P2处与小球B 发生碰撞时的碰前、碰后的瞬间状态.xy 坐标系为观测坐标系，mn坐标系的m 轴为两小球的轴心连线.将小球A、B 碰前的速度vA、vB和碰后的速度v′A、v′B沿m、n方向分解如图4所示.

图4 小球碰撞前后的速度

令两小球的坐标分别为（xA，yA）和（xB，yB），则它们之间的距离

①将速度从xy 坐标系映射到mn 坐标系

如图5所示，将在xy 坐标系中小球A 的速度vAx、vAy映射到mn 坐标系，可得

图5 从xy 坐标系映射到mn 坐标系

对小球B，同理可得

②碰撞后的速度

碰撞前、碰撞后系统的动能相等，得

碰撞前、碰撞后系统的动量相等，得

将式（6）、（7）两个方程联立，可得两小球在m方向上的速度

两小球在n方向上的速度保持不变，即

③将碰撞后的速度从mn 坐标系映射到xy坐标系

如图6所示，将在mn 坐标系中A 小球碰撞后的速度映射到xy 坐标系，可得

图6 从mn坐标系映射到xy 坐标系

对小球B，同理得

（2）两小球碰撞的发生时间tc

为了计算tc，我们以小球B 为坐标系来描述小球A 相对小球B的运动，如图7.

图7 以小球B为参考系的碰撞示意图

两小球之间的距离为D，小球A 以vr的相对速度运动，并经过时间tc后与小球B发生碰撞，此时相对位移S＝vrtc.

由于相对速度vr＝vA－vB，所以，在x 方向上，相对速度vxr＝vAx－vBx，在y 方向上，相对速度vyr＝vAy－vBy，小球A 相对小球B 的相对速率

再由图7，可得如下关系

（rA＋rB）2＝S2＋D2－2SDcosθ

即S2-2DcosθS+D2-（rA+rB）2=0

解方程得

其中，角度θ＝φ－γ；角度φ 为相对位移与水平方向的夹角；γ 则为小球A、B 球心连线与水平方向的夹角.

显然

可以求得

cosθ＝cos（φ－γ）＝cosφcosγ＋sinφsinγ

sinθ＝sin（φ－γ）＝sinφcosγ－cosφsinγ

图8解释了S－和S＋的意义.如图所示，以小球B中心为圆心，半径为rA＋rB的虚线圆轨迹是发生碰撞时小球A 的中心可能处在的位置.小球A 朝小球B相对运动的方向与此虚线圆相交有两点，到小球A 球心的距离分别为S－、S＋.

图8 相对位移的解释

显然，实际碰撞在S－时就发生了，故S＋舍去.

（3）T 时间后小球的位置

各小球T 时间后的位置可以表达如下：

x′A＝xA＋vAxtc＋v′Ax（T－tc）

y′A＝yA＋vAytc＋v′Ay（T－tc）

x′B＝xB＋vBxtc＋v′Bx（T－tc）

y′B＝yB＋vBytc＋v′By（T－tc）

代入前面求得的碰撞之后的速度和碰撞时的时间，即可得到结果.

2 碰撞信息的计算机表达

碰撞的类型有两种：小球与容器壁的碰撞，小球与小球的碰撞.

为此，我们需要定义三个类型：描述小球的Sphere类；描述容器壁的Side类；描述碰撞事件的Collision类.

对于Sphere类，记录了小球的位置坐标（x，y）、x 方向的速度vx、y 方向的速度vy、小球的半径radius、质量mass，调用其成员函数Process（double dt），小球以现有的速度向前走dt 的时间.

对于Side类，用m＿type表示容器壁运动的类型（值为零时表示竖直左壁、值为1时表示竖直右壁、值为2时表示水平上壁、值为3时表示水平下壁），记录了容器壁的位置m＿pos和容器壁的运动速度m＿v，调用其成员函数Process（double dt），容器壁以现有的速度向前走dt的时间.

对于Collision类，记录了碰撞事件的信息，继承于它的BtoBCld类表示小球与小球的碰撞，此时m＿type为零，pSphereA 和pSphereB分别代表A、B 小球；继承于它的BtoSCld类表示小球与容器壁的碰撞，此时m＿type为1，pSphere代表小球，pSide代表容器壁.m＿cldT 表示本碰撞的发生时间.调用成员函数ChangeV（），按照弹性碰撞规律改变所涉及小球的运动速度.

3 基于碰撞事件驱动的算法

前面我们已经讨论了小球与容器壁、小球与小球之间的单次碰撞问题.但是，要处理整个系统中大量小球的相互碰撞运动，将有更多的问题需要考虑.

（1）同一个小球可能发生两次甚至是多次碰撞，图9为同一小球连续两次碰撞的情况.

（2）碰撞传递现象，指的是若小球A 先与小球B碰撞，引起小球B的运动，之后小球B将与别的小球或容器壁碰撞，图10列出了二次传递的情况.

（3）碰撞的并存性.整个系统中存在大量运动的小球，在任一时刻进行检测，都可以找到很多碰撞事件将要发生，我们需要根据其发生的先后次序依次处理.

图9 同一小球发生多次碰撞

图10 碰撞的传递

（4）碰撞的相干性.指对某个先发生的碰撞事件的处理，可能会使已检测到的发生在其后的碰撞事件失效，需要重新检测.

基于以上的分析，我们给出一种基于碰撞事件驱动的算法，其中output＿T 为观测周期，表示每隔多长时间输出图像；current＿t为当前时间，在当前时间到达观测周期时，输出画面并置当前时间为零，重新累计；collision是存放碰撞事件信息的队列.

算法描述具体如下：

（1）遍历所有的小球，对每个小球，找出与周围其他小球和容器壁之间的所有碰撞，对最近发生那次碰撞，如果事件队列collision中没有记录，则记录之.初始当前时间current＿t＝0.

（2）对所有已记录在collision 中的碰撞信息，取出最近的那次碰撞，从队列删除之，并令其碰撞发生时间为min＿cld＿t.

（3）如果系统向前走min＿cld＿t的时间，没有超过图形输出周期output＿T，即current＿t＋min＿cld＿t＜output＿T，到第4步.

如果超过图形输出周期，即current＿t＋min＿cld＿t＞＝output＿T 时，系统向前走output＿T－current＿t的时间达到输出周期，输出系统当前状态下的图形，当前时间current＿t置为0，当前要处理的碰撞的发生时间还剩min＿cld＿t＝min＿cld＿t－（output＿T－current＿t），循环直到当前碰撞发生时间处于下一个图形输出周期内，即current＿t＋min＿cld＿t＜output＿T.

（4）系统向前走min＿cld＿t，当前时间改为current＿t＝current＿t＋min＿cld＿t，当前碰撞已经发生，改变当前碰撞中所涉及小球的速度.

（5）由于当前碰撞所涉及小球的运动速度发生改变，需要重新检测它们的碰撞情况.如果该碰撞为小球与容器壁的碰撞，所涉及小球只有一个；如果是小球与小球之间的碰撞，所涉及小球就将有两个.

对每个这样的小球由于其速度发生改变，以前记录在collision中的所有包含该小球的碰撞事件都失效了，将其删除.注意，如果要删除的碰撞是小球与小球之间的碰撞，删除该碰撞还将影响到其中的另一个小球（即受干扰小球），还需要重新检测该小球与周围的碰撞.

（6）回到第（2）步.

下面举例来说明这种基于碰撞事件驱动的算法的执行过程.本例展示了小球A、B、C、D 以及容器壁之间的碰撞，并显示了碰撞的相干性，指出了受干扰小球，如图11所示.

图11 小球A、B、C、D 以及容器壁之间的碰撞

图11（a）中，各小球的运动速度如箭头所示，其中小球C 的速度为零，按算法第（1）步，此时只能检测到两个有效的碰撞：小球A 与小球B在O点的碰撞，小球C 与小球D 的碰撞，而小球A 与容器壁E的碰撞将舍去.将这两个有效的碰撞事件记录在事件队列collision中，记物体A 与B 之间的碰撞事件为C（AB），则此时事件队列状态为：C（AB）、C（CD）.

从事件队列中取出最先发生的事件，即C（CD），从事件队列中删除之，并让系统中各小球按各自轨迹向前运动，直到C（CD）发生.此时，C小球获得速度，而D 小球停下来，如图11（b）.由于小球C速度发生改变，删除队列collision中所有与C小球有关的碰撞事件，并重新检测小球C与系统中其他小球的碰撞.我们会找到小球C 与小球B之间的碰撞事件，将其存入事件队列中，此时事件队列状态为：C（AB）、C（BC）.

从事件队列中取出最先发生的事件，即C（BC），从事件队列中删除之，并让系统继续向前运动，直到C（BC）发生，并改变小球B、C 的运动速度，如图11（c）.小球B的速度改变之后，需要从队列中删除所有与小球B 有关的事件，从而把先前图11（a）中保存在队列中的C（AB）事件给删除了.由于该事件的删除，队列中不再存有与小球A相关的碰撞事件了，所以称小球A 为受干扰小球，此时需要重新检测小球A 的最近碰撞事件并加入队列.同时检测小球B 的最近碰撞事件和小球C的最近碰撞事件，加入事件队列.如此运行下去.

由算法第（1）步，可以看到，对每个小球，只记录了其最近发生的那次碰撞的信息，所以我们在队列中只需要保存O（n）的碰撞事件.本算法在第（1）步初始化系统时，需要O（n2）的时间，但在运行起来之后，处理每个碰撞，只需O（n）的计算量，处理nc个碰撞事件所需要的时间仅为nc（k0＋k1n），k0，k1为常数.

上述算法的特点：（1）本算法是以碰撞事件为驱动的算法，即对任一速度不为零的小球，在事件队列中总能找到与之相关的碰撞事件.为了保证这一点，需要在进入动态模拟之前，检测出所有小球的最近碰撞信息.在动态模拟过程中，小球如果由于碰撞而使速度发生变化，需要重新检测它以及受影响的小球，找出其与周围小球或容器壁的最近碰撞事件并记录在事件队列里.（2）本算法对给定初始状态的小球系统，能准确计算出经过一段时间T 后各小球的状态.（3）本算法的运行效率为nc（k0＋k1n），处理单次碰撞的时间与系统中小球数目呈线性关系，能够实时模拟大量小球的相互碰撞运动.

将上述方法应用于模拟大量气体分子的无规则运动，如图12，其中折线表示跟踪单个分子的运动轨迹，图中有800个半径为4像素的小球，实时模拟效果流畅，小球之间作弹性碰撞，没有观察到穿透和重叠现象.图13为相应的速率分布，可以看到它与麦克斯韦速率分布曲线吻合.在主频2.66GHz、内存512M 的计算机上进行的运行测试中，小球数目设置为2000，小球半径为1像素，运动效果仍然流畅.本方法已经在《大学物理学V4.0》计算机模拟教学软件［3］中得到应用，取得了良好的效果.

［1］胡其图.大学物理学V3.0［M］.北京：高等教育出版社，2003

［2］Dong－Jin Kim，Leonidas J.Guibas and Sung－Yong Shin.Fast Collision Detection Among Multiple Moving Spheres.IEEE Transactions on Visualization and Computer Graphics，1998，4（3）：230～242

［3］胡其图，张偶利，邓晓，张超，张小灵.大学物理学V4.0［M］.北京：高等教育出版社，2006