大量小球無(wú)規(guī)則運(yùn)動(dòng)計(jì)算機(jī)模擬的可視化

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

（1上海交通大學(xué)物理系，上海 200240）

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

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

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

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

1 單次碰撞的碰撞檢測(cè)與碰撞反應(yīng)

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

1.1 小球與容器壁的碰撞

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

如圖1所示，左容器壁處于xw位置，運(yùn)動(dòng)速度為vw.小球初始位置為（x，y），碰撞前在x 方向上的運(yùn)動(dòng)速度為vx.兩者經(jīng)過(guò)時(shí)間tc后相撞，此時(shí)小球球心處于xc位置，而容器壁也將達(dá)到粗虛線位置，與小球相切.碰撞后小球在x 方向上的速度為v′x，小球在y 方向上的速度不變.再經(jīng)過(guò)時(shí)間tr＝T－tc，小球?qū)⒌竭_(dá)x′位置.

圖1 小球與容器壁的碰撞

T 時(shí)間后小球和器壁的位置如下：

其中，xc=x+vxtc.

顯然，我們需要先得到碰撞后小球在x 方向的速度v′x和碰撞時(shí)的時(shí)間tc.

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

碰撞前在x 方向上小球相對(duì)容器壁的速度為

由于碰撞后容器壁速度不變，所以碰撞后在x 方向上小球相對(duì)器壁的速度為

v′xr＝v′x－vw

顯然ｖ′ｘｒ＝－ｖｘｒ

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

（2）碰撞時(shí)的時(shí)間tc

設(shè)定小球球心可以到達(dá)的邊界位置為xb，則

R 的絕對(duì)值為小球半徑.當(dāng)R 為正數(shù)時(shí)，為左邊容器壁；當(dāng)R 為負(fù)數(shù)時(shí)，為右邊容器壁.

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

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

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

取負(fù)號(hào)時(shí)表示沿x 的負(fù)方向.

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

①當(dāng)容器壁為左邊容器壁，即dx＞0 時(shí)，dx≤－vxrT；

②當(dāng)容器壁為右邊容器壁，即dx＜0 時(shí)，－dx≤vxrT.

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

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

1.2 小球與小球之間的碰撞

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

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

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

（1）小球碰撞后的速度

圖4展示了小球A 在位置P2處與小球B 發(fā)生碰撞時(shí)的碰前、碰后的瞬間狀態(tài).xy 坐標(biāo)系為觀測(cè)坐標(biāo)系，mn坐標(biāo)系的m 軸為兩小球的軸心連線.將小球A、B 碰前的速度vA、vB和碰后的速度v′A、v′B沿m、n方向分解如圖4所示.

圖4 小球碰撞前后的速度

令兩小球的坐標(biāo)分別為（xA，yA）和（xB，yB），則它們之間的距離

①將速度從xy 坐標(biāo)系映射到mn 坐標(biāo)系

如圖5所示，將在xy 坐標(biāo)系中小球A 的速度vAx、vAy映射到mn 坐標(biāo)系，可得

圖5 從xy 坐標(biāo)系映射到mn 坐標(biāo)系

對(duì)小球B，同理可得

②碰撞后的速度

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

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

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

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

③將碰撞后的速度從mn 坐標(biāo)系映射到xy坐標(biāo)系

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

圖6 從mn坐標(biāo)系映射到xy 坐標(biāo)系

對(duì)小球B，同理得

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

為了計(jì)算tc，我們以小球B 為坐標(biāo)系來(lái)描述小球A 相對(duì)小球B的運(yùn)動(dòng)，如圖7.

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

兩小球之間的距離為D，小球A 以vr的相對(duì)速度運(yùn)動(dòng)，并經(jīng)過(guò)時(shí)間tc后與小球B發(fā)生碰撞，此時(shí)相對(duì)位移S＝vrtc.

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

再由圖7，可得如下關(guān)系

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

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

解方程得

其中，角度θ＝φ－γ；角度φ 為相對(duì)位移與水平方向的夾角；γ 則為小球A、B 球心連線與水平方向的夾角.

顯然

可以求得

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

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

圖8解釋了S－和S＋的意義.如圖所示，以小球B中心為圓心，半徑為rA＋rB的虛線圓軌跡是發(fā)生碰撞時(shí)小球A 的中心可能處在的位置.小球A 朝小球B相對(duì)運(yùn)動(dòng)的方向與此虛線圓相交有兩點(diǎn)，到小球A 球心的距離分別為S－、S＋.

圖8 相對(duì)位移的解釋

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

（3）T 時(shí)間后小球的位置

各小球T 時(shí)間后的位置可以表達(dá)如下：

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）

代入前面求得的碰撞之后的速度和碰撞時(shí)的時(shí)間，即可得到結(jié)果.

2 碰撞信息的計(jì)算機(jī)表達(dá)

碰撞的類(lèi)型有兩種：小球與容器壁的碰撞，小球與小球的碰撞.

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

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

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

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

3 基于碰撞事件驅(qū)動(dòng)的算法

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

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

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

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

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

圖10 碰撞的傳遞

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

基于以上的分析，我們給出一種基于碰撞事件驅(qū)動(dòng)的算法，其中output＿T 為觀測(cè)周期，表示每隔多長(zhǎng)時(shí)間輸出圖像；current＿t為當(dāng)前時(shí)間，在當(dāng)前時(shí)間到達(dá)觀測(cè)周期時(shí)，輸出畫(huà)面并置當(dāng)前時(shí)間為零，重新累計(jì)；collision是存放碰撞事件信息的隊(duì)列.

算法描述具體如下：

（1）遍歷所有的小球，對(duì)每個(gè)小球，找出與周?chē)渌∏蚝腿萜鞅谥g的所有碰撞，對(duì)最近發(fā)生那次碰撞，如果事件隊(duì)列collision中沒(méi)有記錄，則記錄之.初始當(dāng)前時(shí)間current＿t＝0.

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

（3）如果系統(tǒng)向前走min＿cld＿t的時(shí)間，沒(méi)有超過(guò)圖形輸出周期output＿T，即current＿t＋min＿cld＿t＜o(jì)utput＿T，到第4步.

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

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

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

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

（6）回到第（2）步.

下面舉例來(lái)說(shuō)明這種基于碰撞事件驅(qū)動(dòng)的算法的執(zhí)行過(guò)程.本例展示了小球A、B、C、D 以及容器壁之間的碰撞，并顯示了碰撞的相干性，指出了受干擾小球，如圖11所示.

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

圖11（a）中，各小球的運(yùn)動(dòng)速度如箭頭所示，其中小球C 的速度為零，按算法第（1）步，此時(shí)只能檢測(cè)到兩個(gè)有效的碰撞：小球A 與小球B在O點(diǎn)的碰撞，小球C 與小球D 的碰撞，而小球A 與容器壁E的碰撞將舍去.將這兩個(gè)有效的碰撞事件記錄在事件隊(duì)列collision中，記物體A 與B 之間的碰撞事件為C（AB），則此時(shí)事件隊(duì)列狀態(tài)為：C（AB）、C（CD）.

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

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

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

上述算法的特點(diǎn)：（1）本算法是以碰撞事件為驅(qū)動(dòng)的算法，即對(duì)任一速度不為零的小球，在事件隊(duì)列中總能找到與之相關(guān)的碰撞事件.為了保證這一點(diǎn)，需要在進(jìn)入動(dòng)態(tài)模擬之前，檢測(cè)出所有小球的最近碰撞信息.在動(dòng)態(tài)模擬過(guò)程中，小球如果由于碰撞而使速度發(fā)生變化，需要重新檢測(cè)它以及受影響的小球，找出其與周?chē)∏蚧蛉萜鞅诘淖罱鲎彩录⒂涗浽谑录?duì)列里.（2）本算法對(duì)給定初始狀態(tài)的小球系統(tǒng)，能準(zhǔn)確計(jì)算出經(jīng)過(guò)一段時(shí)間T 后各小球的狀態(tài).（3）本算法的運(yùn)行效率為nc（k0＋k1n），處理單次碰撞的時(shí)間與系統(tǒng)中小球數(shù)目呈線性關(guān)系，能夠?qū)崟r(shí)模擬大量小球的相互碰撞運(yùn)動(dòng).

將上述方法應(yīng)用于模擬大量氣體分子的無(wú)規(guī)則運(yùn)動(dòng)，如圖12，其中折線表示跟蹤單個(gè)分子的運(yùn)動(dòng)軌跡，圖中有800個(gè)半徑為4像素的小球，實(shí)時(shí)模擬效果流暢，小球之間作彈性碰撞，沒(méi)有觀察到穿透和重疊現(xiàn)象.圖13為相應(yīng)的速率分布，可以看到它與麥克斯韋速率分布曲線吻合.在主頻2.66GHz、內(nèi)存512M 的計(jì)算機(jī)上進(jìn)行的運(yùn)行測(cè)試中，小球數(shù)目設(shè)置為2000，小球半徑為1像素，運(yùn)動(dòng)效果仍然流暢.本方法已經(jīng)在《大學(xué)物理學(xué)V4.0》計(jì)算機(jī)模擬教學(xué)軟件［3］中得到應(yīng)用，取得了良好的效果.

［1］胡其圖.大學(xué)物理學(xué)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］胡其圖，張偶利，鄧曉，張超，張小靈.大學(xué)物理學(xué)V4.0［M］.北京：高等教育出版社，2006