鄧 曉 胡其圖 李發(fā)澤 張小靈

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

（2寧夏大學物理電氣信息學院，寧夏 銀川 750021）

理想氣體分子熱運動、布朗運動、氣體分子擴散等物理現(xiàn)象的模擬，都可以抽象為大量小球無規(guī)則運動的形式.使用計算機對這種運動形式進行實時模擬有著廣泛的應用.由于系統(tǒng)中小球數(shù)目眾多，小球之間存在大量的彈性碰撞，實現(xiàn)實時模擬可視化的關鍵在于如何減少計算量、提高運行效率.

一種常用的辦法是：系統(tǒng)按固定時間步長向前運動，在系統(tǒng)準備向前走一步時，先對所有的小球掃描一遍，檢測它們之間是否將發(fā)生重疊.如果將發(fā)生重疊，則認為會有碰撞發(fā)生，按照彈性碰撞所滿足的規(guī)律改變相關小球的運動速度.這種算法的每次掃描過程需要O（n2）的運算量，在小球數(shù)目較多時運算量很大.另外，由于所取的時間步長不能無限小，當小球運動速度較快或小球半徑相對較小時，會使得本來應該檢測到的碰撞沒有被檢測到，從而在實現(xiàn)實時模擬可視化時容易出現(xiàn)兩個小球相互重疊或穿透的情形.本方法由于運行效率不高，模擬的小球數(shù)目到達600個時系統(tǒng)就很難運行了，并且小球之間的重疊和穿透現(xiàn)象無法避免［1］.

Dong－Jin Kim 提出一種基于碰撞事件驅動的模擬方法［2］，我們利用此方法實現(xiàn)了大量小球無重疊、無穿透運動的計算機模擬，系統(tǒng)不再按固定時間步長運行，而是按實際碰撞發(fā)生的時間來運行.這種方法有效地解決了上述問題，并大大提高了計算效率.考慮一個擁有大量小球的二維封閉系統(tǒng)，其邊界是可移動的（例如：汽缸的活塞），因而也具有一定的運動速度，下面我們來具體討論本模擬方法.

1 單次碰撞的碰撞檢測與碰撞反應

系統(tǒng)中包含了兩種類型的碰撞：小球與容器壁的碰撞、小球與小球之間的碰撞.下面分別討論它們的碰撞檢測與碰撞反應.

1.1 小球與容器壁的碰撞

一個容器在二維空間中可用一個長方形來表示，容器壁有四個：左右兩個豎直壁和上下兩個水平壁.小球與水平壁之間的碰撞和小球與豎直壁之間的碰撞算法原理相同，本文只給出小球與豎直左壁在一個時間步長T 內發(fā)生碰撞的情況.

如圖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 為正數(shù)時，為左邊容器壁；當R 為負數(shù)時，為右邊容器壁.

以容器壁為參考系，小球相對容器壁的碰撞過程如圖2.容器壁靜止，小球在x 方向上以相對速度vxr靠近容器壁.

圖2 以容器壁為參考系的碰撞示意圖

dx表示碰撞前的相對位移，d′x為碰撞后的位移，顯然

取負號時表示沿x 的負方向.

在時間步長T 中，小球與器壁發(fā)生碰撞的條件為：

①當容器壁為左邊容器壁，即dx＞0 時，dx≤－vxrT；

②當容器壁為右邊容器壁，即dx＜0 時，－dx≤vxrT.

如果發(fā)生碰撞，則發(fā)生碰撞的時間

聯(lián)合式（1）、（2）、（3），得

1.2 小球與小球之間的碰撞

假設有兩個半徑不同、質量不同的小球相撞，如圖3展示了在T 時間步長里小球A 與小球B發(fā)生碰撞的情況.

圖3 小球與小球之間的碰撞

對于小球A，它從P1位置以vA的速度經過時間tc后運動到P2的位置并與小球B發(fā)生碰撞，速度改變?yōu)関′A，然后再以碰撞后的速度繼續(xù)運動tr的時間到達P3位置，其中tr＝T－tc.小球B 的運動也是如此.如果我們能夠計算出碰撞的發(fā)生時間tc以及碰撞后小球A、B的速度，那么整個過程就被求解了.

（1）小球碰撞后的速度

圖4展示了小球A 在位置P2處與小球B 發(fā)生碰撞時的碰前、碰后的瞬間狀態(tài).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，同理可得

②碰撞后的速度

碰撞前、碰撞后系統(tǒng)的動能相等，得

碰撞前、碰撞后系統(tǒng)的動量相等，得

將式（6）、（7）兩個方程聯(lián)立，可得兩小球在m方向上的速度

兩小球在n方向上的速度保持不變，即

③將碰撞后的速度從mn 坐標系映射到xy坐標系

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

圖6 從mn坐標系映射到xy 坐標系

對小球B，同理得

（2）兩小球碰撞的發(fā)生時間tc

為了計算tc，我們以小球B 為坐標系來描述小球A 相對小球B的運動，如圖7.

圖7 以小球B為參考系的碰撞示意圖

兩小球之間的距離為D，小球A 以vr的相對速度運動，并經過時間tc后與小球B發(fā)生碰撞，此時相對位移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的虛線圓軌跡是發(fā)生碰撞時小球A 的中心可能處在的位置.小球A 朝小球B相對運動的方向與此虛線圓相交有兩點，到小球A 球心的距離分別為S－、S＋.

圖8 相對位移的解釋

顯然，實際碰撞在S－時就發(fā)生了，故S＋舍去.

（3）T 時間后小球的位置

各小球T 時間后的位置可以表達如下：

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

y′A＝y(tǒng)A＋vAytc＋v′Ay（T－tc）

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

y′B＝y(tǒng)B＋vBytc＋v′By（T－tc）

代入前面求得的碰撞之后的速度和碰撞時的時間，即可得到結果.

2 碰撞信息的計算機表達

碰撞的類型有兩種：小球與容器壁的碰撞，小球與小球的碰撞.

為此，我們需要定義三個類型：描述小球的Sphere類；描述容器壁的Side類；描述碰撞事件的Collision類.

對于Sphere類，記錄了小球的位置坐標（x，y）、x 方向的速度vx、y 方向的速度vy、小球的半徑radius、質量mass，調用其成員函數(shù)Process（double dt），小球以現(xiàn)有的速度向前走dt 的時間.

對于Side類，用m＿type表示容器壁運動的類型（值為零時表示豎直左壁、值為1時表示豎直右壁、值為2時表示水平上壁、值為3時表示水平下壁），記錄了容器壁的位置m＿pos和容器壁的運動速度m＿v，調用其成員函數(shù)Process（double dt），容器壁以現(xiàn)有的速度向前走dt的時間.

對于Collision類，記錄了碰撞事件的信息，繼承于它的BtoBCld類表示小球與小球的碰撞，此時m＿type為零，pSphereA 和pSphereB分別代表A、B 小球；繼承于它的BtoSCld類表示小球與容器壁的碰撞，此時m＿type為1，pSphere代表小球，pSide代表容器壁.m＿cldT 表示本碰撞的發(fā)生時間.調用成員函數(shù)ChangeV（），按照彈性碰撞規(guī)律改變所涉及小球的運動速度.

3 基于碰撞事件驅動的算法

前面我們已經討論了小球與容器壁、小球與小球之間的單次碰撞問題.但是，要處理整個系統(tǒng)中大量小球的相互碰撞運動，將有更多的問題需要考慮.

（1）同一個小球可能發(fā)生兩次甚至是多次碰撞，圖9為同一小球連續(xù)兩次碰撞的情況.

（2）碰撞傳遞現(xiàn)象，指的是若小球A 先與小球B碰撞，引起小球B的運動，之后小球B將與別的小球或容器壁碰撞，圖10列出了二次傳遞的情況.

（3）碰撞的并存性.整個系統(tǒng)中存在大量運動的小球，在任一時刻進行檢測，都可以找到很多碰撞事件將要發(fā)生，我們需要根據(jù)其發(fā)生的先后次序依次處理.

圖9 同一小球發(fā)生多次碰撞

圖10 碰撞的傳遞

（4）碰撞的相干性.指對某個先發(fā)生的碰撞事件的處理，可能會使已檢測到的發(fā)生在其后的碰撞事件失效，需要重新檢測.

基于以上的分析，我們給出一種基于碰撞事件驅動的算法，其中output＿T 為觀測周期，表示每隔多長時間輸出圖像；current＿t為當前時間，在當前時間到達觀測周期時，輸出畫面并置當前時間為零，重新累計；collision是存放碰撞事件信息的隊列.

算法描述具體如下：

（1）遍歷所有的小球，對每個小球，找出與周圍其他小球和容器壁之間的所有碰撞，對最近發(fā)生那次碰撞，如果事件隊列collision中沒有記錄，則記錄之.初始當前時間current＿t＝0.

（2）對所有已記錄在collision 中的碰撞信息，取出最近的那次碰撞，從隊列刪除之，并令其碰撞發(fā)生時間為min＿cld＿t.

（3）如果系統(tǒng)向前走min＿cld＿t的時間，沒有超過圖形輸出周期output＿T，即current＿t＋min＿cld＿t＜output＿T，到第4步.

如果超過圖形輸出周期，即current＿t＋min＿cld＿t＞＝output＿T 時，系統(tǒng)向前走output＿T－current＿t的時間達到輸出周期，輸出系統(tǒng)當前狀態(tài)下的圖形，當前時間current＿t置為0，當前要處理的碰撞的發(fā)生時間還剩min＿cld＿t＝min＿cld＿t－（output＿T－current＿t），循環(huán)直到當前碰撞發(fā)生時間處于下一個圖形輸出周期內，即current＿t＋min＿cld＿t＜output＿T.

（4）系統(tǒng)向前走min＿cld＿t，當前時間改為current＿t＝current＿t＋min＿cld＿t，當前碰撞已經發(fā)生，改變當前碰撞中所涉及小球的速度.

（5）由于當前碰撞所涉及小球的運動速度發(fā)生改變，需要重新檢測它們的碰撞情況.如果該碰撞為小球與容器壁的碰撞，所涉及小球只有一個；如果是小球與小球之間的碰撞，所涉及小球就將有兩個.

對每個這樣的小球由于其速度發(fā)生改變，以前記錄在collision中的所有包含該小球的碰撞事件都失效了，將其刪除.注意，如果要刪除的碰撞是小球與小球之間的碰撞，刪除該碰撞還將影響到其中的另一個小球（即受干擾小球），還需要重新檢測該小球與周圍的碰撞.

（6）回到第（2）步.

下面舉例來說明這種基于碰撞事件驅動的算法的執(zhí)行過程.本例展示了小球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），則此時事件隊列狀態(tài)為：C（AB）、C（CD）.

從事件隊列中取出最先發(fā)生的事件，即C（CD），從事件隊列中刪除之，并讓系統(tǒng)中各小球按各自軌跡向前運動，直到C（CD）發(fā)生.此時，C小球獲得速度，而D 小球停下來，如圖11（b）.由于小球C速度發(fā)生改變，刪除隊列collision中所有與C小球有關的碰撞事件，并重新檢測小球C與系統(tǒng)中其他小球的碰撞.我們會找到小球C 與小球B之間的碰撞事件，將其存入事件隊列中，此時事件隊列狀態(tài)為：C（AB）、C（BC）.

從事件隊列中取出最先發(fā)生的事件，即C（BC），從事件隊列中刪除之，并讓系統(tǒng)繼續(xù)向前運動，直到C（BC）發(fā)生，并改變小球B、C 的運動速度，如圖11（c）.小球B的速度改變之后，需要從隊列中刪除所有與小球B 有關的事件，從而把先前圖11（a）中保存在隊列中的C（AB）事件給刪除了.由于該事件的刪除，隊列中不再存有與小球A相關的碰撞事件了，所以稱小球A 為受干擾小球，此時需要重新檢測小球A 的最近碰撞事件并加入隊列.同時檢測小球B 的最近碰撞事件和小球C的最近碰撞事件，加入事件隊列.如此運行下去.

由算法第（1）步，可以看到，對每個小球，只記錄了其最近發(fā)生的那次碰撞的信息，所以我們在隊列中只需要保存O（n）的碰撞事件.本算法在第（1）步初始化系統(tǒng)時，需要O（n2）的時間，但在運行起來之后，處理每個碰撞，只需O（n）的計算量，處理nc個碰撞事件所需要的時間僅為nc（k0＋k1n），k0，k1為常數(shù).

上述算法的特點：（1）本算法是以碰撞事件為驅動的算法，即對任一速度不為零的小球，在事件隊列中總能找到與之相關的碰撞事件.為了保證這一點，需要在進入動態(tài)模擬之前，檢測出所有小球的最近碰撞信息.在動態(tài)模擬過程中，小球如果由于碰撞而使速度發(fā)生變化，需要重新檢測它以及受影響的小球，找出其與周圍小球或容器壁的最近碰撞事件并記錄在事件隊列里.（2）本算法對給定初始狀態(tài)的小球系統(tǒng)，能準確計算出經過一段時間T 后各小球的狀態(tài).（3）本算法的運行效率為nc（k0＋k1n），處理單次碰撞的時間與系統(tǒng)中小球數(shù)目呈線性關系，能夠實時模擬大量小球的相互碰撞運動.

將上述方法應用于模擬大量氣體分子的無規(guī)則運動，如圖12，其中折線表示跟蹤單個分子的運動軌跡，圖中有800個半徑為4像素的小球，實時模擬效果流暢，小球之間作彈性碰撞，沒有觀察到穿透和重疊現(xiàn)象.圖13為相應的速率分布，可以看到它與麥克斯韋速率分布曲線吻合.在主頻2.66GHz、內存512M 的計算機上進行的運行測試中，小球數(shù)目設置為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