二維粒子運(yùn)動(dòng)模型的模擬和解析研究

樓鑫鑫,丁中俊,田 波

(1.合肥工業(yè)大學(xué) 汽車與交通工程學(xué)院,安徽 合肥 230009; 2.安徽農(nóng)業(yè)大學(xué) 工學(xué)院,安徽 合肥 230036)

驅(qū)動(dòng)擴(kuò)散系統(tǒng)因其豐富而復(fù)雜的現(xiàn)象而備受關(guān)注[1]。相分離是在驅(qū)動(dòng)擴(kuò)散系統(tǒng)[2]中觀察到的一種重要現(xiàn)象,在許多一維或準(zhǔn)一維模型中都有提及,如ABC模型[3]、具有慢啟動(dòng)規(guī)則的交通流模型[4]、雙車道模型[5]、一維硬粒子系統(tǒng)[6]等。

許多研究集中在二維驅(qū)動(dòng)的擴(kuò)散系統(tǒng),文獻(xiàn)[7]在二維周期網(wǎng)絡(luò)上研究了豐富的相圖,其中2種粒子以相反的方向運(yùn)動(dòng);文獻(xiàn)[8]在二維非對稱簡單排它過程中觀察到2種干擾相;文獻(xiàn)[9]在二維系統(tǒng)中研究了流體力學(xué)對液體混合物相分離后期動(dòng)力學(xué)的影響;文獻(xiàn)[10]在二維Lennard-Jones流體中分析了非平衡通道的形成,該流體由2個(gè)相反方向運(yùn)動(dòng)的顆粒組成;文獻(xiàn)[11]提出2種膠體顆粒的二維模型,表明在驅(qū)動(dòng)膠體系統(tǒng)中的通道形成,是宏觀范圍內(nèi)的相變,但從無序的初始條件開始,在有限時(shí)間內(nèi)不會(huì)發(fā)生宏觀相分離;文獻(xiàn)[12]的二維格子氣模型中存在2種粒子,在外力的反向驅(qū)動(dòng)下研究了“云”的粗化和動(dòng)態(tài)標(biāo)度;二維模型大多采用隨機(jī)更新規(guī)則,考慮到一些應(yīng)用背景(如行人對向流),文獻(xiàn)[13]建立了一個(gè)具有并行更新規(guī)則與存在位置交換的二維驅(qū)動(dòng)擴(kuò)散系統(tǒng),在該模型中,2種類型的粒子在二維網(wǎng)絡(luò)上以相反的方向運(yùn)動(dòng),當(dāng)粒子發(fā)生沖突時(shí)根據(jù)周圍環(huán)境判斷是否發(fā)生位置交換,然后采用簡單平均場方法對存在位置交換的二維驅(qū)動(dòng)擴(kuò)散系統(tǒng)進(jìn)行了解析,但該方法并沒有考慮粒子之間的相關(guān)性而使得誤差較大。

本文以文獻(xiàn)[13]為基礎(chǔ),建立存在位置交換的對向粒子流模型。然而粒子在實(shí)際的運(yùn)動(dòng)過程中存在相互作用,學(xué)者們將相關(guān)性[14]引入到簡單平均場理論,得到簇平均場理論[15],也稱作cluster平均場方法。本文采用簇平均場方法對模型進(jìn)行解析,在解析中考慮粒子間的相關(guān)性,解析結(jié)果與仿真結(jié)果在均勻態(tài)下較為吻合,當(dāng)出現(xiàn)堵塞帶時(shí)誤差仍然較大。

1 對向粒子流模型

本文基于元胞自動(dòng)機(jī)建立對向粒子流模型,如圖1所示。在一個(gè)網(wǎng)絡(luò)中有N×N個(gè)格點(diǎn),每個(gè)格點(diǎn)被一個(gè)粒子占據(jù)或者為空。系統(tǒng)中存在2類粒子,每個(gè)粒子可以向3個(gè)方向運(yùn)動(dòng)(圖1a)。圖1a中:PE(實(shí)心圓)為第1類粒子,向東、向北、向南跳躍的概率分別為q、(1-q)/2、(1-q)/2;PW(空心圓)為第2類粒子,向西、向北、向南跳躍的概率分別為q、(1-q)/2、(1-q)/2。

系統(tǒng)中有數(shù)量相等的2類粒子,粒子總密度為ρ,粒子的更新規(guī)則為并行更新,即每一個(gè)時(shí)間步,所有粒子同時(shí)進(jìn)行更新;系統(tǒng)的邊界為周期性邊界,即系統(tǒng)中粒子數(shù)目保持不變。每個(gè)粒子先由其運(yùn)動(dòng)概率的大小來確定運(yùn)動(dòng)方向,再判斷粒子運(yùn)動(dòng)的目標(biāo)格點(diǎn)是否為空,若格點(diǎn)為空且同時(shí)沒有其他粒子試圖進(jìn)入該格點(diǎn),則該粒子在下一個(gè)時(shí)間步進(jìn)入該格點(diǎn);在同一時(shí)間步中,若有n個(gè)粒子可以進(jìn)入同一個(gè)空格點(diǎn),則以1/n的概率隨機(jī)選擇一個(gè)粒子進(jìn)入該格點(diǎn)。此外,模型中還存在位置交換的規(guī)則,圖1b中:在t時(shí)刻,若格點(diǎn)3的粒子PE試圖向東運(yùn)動(dòng),格點(diǎn)4的粒子PW試圖向西運(yùn)動(dòng),且1、2、5、6都被占據(jù)(含對角線的圓表示該格點(diǎn)可能被PE也可能被PW占據(jù)),則在t+1時(shí)間步這2個(gè)粒子以pe概率發(fā)生位置的交換;若1、2、5、6中有格點(diǎn)為空,則格點(diǎn)3、4的粒子不會(huì)發(fā)生位置交換。圖1b的右圖表示t+1時(shí)刻3、4格點(diǎn)的2個(gè)粒子成功發(fā)生交換,?表示1、2、5、6的狀態(tài)不確定。

圖1 對向粒子流模型

2 結(jié)果和討論

2.1 模擬仿真

本文在一個(gè)尺寸為100×100的周期邊界網(wǎng)格上進(jìn)行仿真,即系統(tǒng)中存在100×100個(gè)格子。首先舍去前107個(gè)時(shí)間步使系統(tǒng)達(dá)到穩(wěn)定,然后對后面的106個(gè)時(shí)間步進(jìn)行統(tǒng)計(jì),對每一種情況進(jìn)行20 次仿真并計(jì)算平均值,得到考慮位置交換下粒子平均速度與密度的關(guān)系,如圖2所示。

圖2中:平均速度為粒子在每個(gè)時(shí)間步往目標(biāo)方向運(yùn)動(dòng)的格點(diǎn)數(shù);密度為系統(tǒng)中粒子數(shù)量占比。從圖2可以看出,系統(tǒng)在任何密度下未出現(xiàn)完全堵塞,這是由于平均速度沒有降至0。圖2a中,粒子平均速度隨著密度增加而一直下降,且系統(tǒng)中存在一個(gè)臨界密度,當(dāng)粒子密度超過這個(gè)值時(shí),粒子的平均運(yùn)動(dòng)速度下降得更快。圖2b～圖2d中,平均速度隨密度增加而下降,達(dá)到臨界密度之后,平均速度增加,這是由于q值與pe值較大,當(dāng)粒子數(shù)目較多時(shí),相遇粒子之間的互換愈發(fā)頻繁,使粒子可以運(yùn)動(dòng)。相比于圖2b～圖2d,圖2a中粒子平均速度隨著密度增加持續(xù)下降，這是由于此時(shí)q值與pe值較小,相遇的粒子之間發(fā)生位置互換的概率較小,從而相遇的粒子可以運(yùn)動(dòng)的概率也較小;當(dāng)密度小于臨界密度時(shí),粒子處于均勻的狀態(tài),當(dāng)系統(tǒng)中粒子較多時(shí),系統(tǒng)中產(chǎn)生了一個(gè)堵塞帶,阻止了系統(tǒng)運(yùn)動(dòng),如圖3所示。

圖2 在考慮位置交換的情況下粒子運(yùn)動(dòng)平均速度與密度關(guān)系

圖3 q=0.6、pe=0.3的位形

2.2 2-cluster平均場解析法

本文使用2-cluster平均場理論方法解析對向粒子流模型,并分析粒子流的特性。在cluster平均場理論中,連續(xù)的n個(gè)格點(diǎn)即為n-cluster,相鄰的2個(gè)格點(diǎn)即為1個(gè)2-cluster,記做(σi,σi+1)。σi=0,1,2(i=1,2),σi=0表示該格點(diǎn)為空;σi=1表示該格點(diǎn)被目標(biāo)方向?yàn)闁|向的粒子占據(jù);σi=2表示該格點(diǎn)被目標(biāo)方向?yàn)槲飨虻牧Ｗ诱紦?jù)。當(dāng)vmax=1時(shí),t+1時(shí)刻2-cluster(σi,σi+1)狀態(tài)取決于t時(shí)刻4-cluster(σi-1,σi,σi+1,σi+2)的狀態(tài)。當(dāng)系統(tǒng)達(dá)到穩(wěn)態(tài)后,出現(xiàn)2-cluster的概率為p(σ1,σ2)。

vmax=1時(shí),2-cluster方程近似表達(dá)為:

p(σi,σi+1;t+1)=

p(τi-1,τi,τi+1,τi+2;t)

(1)

其中:p(τi-1,τi,τi+1,τi+2;t)為t時(shí)刻4-cluster(τi-1,τi,τi+1,τi+2)存在的概率;∑(σi,σi+1|τi-1,τi,τi+1,τi+2)為所有的4-cluster(τi-1,τi,τi+1,τi+2)更新成為2-cluster(σi,σi+1)的轉(zhuǎn)化概率。

4-cluster(σi-1,σi,σi+1,σi+2)一共有81種位形可更新得到2-cluster(σi,σi+1),在不考慮位置交換的情況下,只有54種情況滿足2-cluster平均場理論方程,如圖4所示。圖4中:第1列為粒子在t時(shí)刻的位形;第2列為粒子在t+1時(shí)刻的位形;第3列與第4列的乘積表示相應(yīng)的轉(zhuǎn)化概率。在系統(tǒng)中隨機(jī)選擇1個(gè)格點(diǎn),格點(diǎn)為空的概率為1-ρ,被占據(jù)的概率為ρ。三角形和圓圈分別對應(yīng)(σi,σi+1)中的σi與σi+1格點(diǎn),且狀態(tài)為空。箭頭表示粒子運(yùn)動(dòng)的目標(biāo)方向,→表示粒子目標(biāo)方向?yàn)闁|向,←表示粒子目標(biāo)方向?yàn)槲飨??表示格點(diǎn)可能為空也可能被占據(jù);ρ1表示東向粒子的密度,ρ2表示西向粒子的密度,2類粒子數(shù)相同時(shí)ρ1=ρ2;qs為粒子向其他方向運(yùn)動(dòng)的概率。

圖4 能轉(zhuǎn)化為p(1,0)的 4-cluster 的示意圖及對應(yīng)的轉(zhuǎn)化概率

將格點(diǎn)進(jìn)行1—12的編號,然后通過p(0,0,0,0)的例子介紹轉(zhuǎn)化概率,如圖5所示,其他概率以此類推。

圖5中:t時(shí)刻1—4格點(diǎn)為空,5—12格點(diǎn)情況未知,可能被占據(jù)也可能為空;在t+1時(shí)刻有1個(gè)東向的粒子進(jìn)入格點(diǎn)2,格點(diǎn)3依舊為空,即本文所求的p(1,0)。在模型中偏移概率的影響較小,為了解析的方便與清晰,忽略部分轉(zhuǎn)移概率中偏移方向粒子的影響。

圖5 圖4中子圖位形演化過程示意圖

2ρ1qs表示格點(diǎn)2在t時(shí)刻為空,在t+1時(shí)刻被東向粒子占據(jù)的概率。在t時(shí)刻只有格點(diǎn)6或者格點(diǎn)10的東向粒子能在t+1時(shí)刻進(jìn)入格點(diǎn)2號。格點(diǎn)6或者格點(diǎn)10被東向粒子占據(jù)的概率為ρ1,該粒子在t+1時(shí)刻進(jìn)入格點(diǎn)2的概率為qs。因此格點(diǎn)2在t+1時(shí)刻被東向粒子占據(jù)概率為2ρ1qs。

1-2ρ1qs-2ρ2qs表示格點(diǎn)3在t時(shí)刻與t+1時(shí)刻都為空的概率。因?yàn)樵趖時(shí)刻只有格點(diǎn)7或者格點(diǎn)11的粒子能在t+1時(shí)刻進(jìn)入格點(diǎn)3，所以存在2種情況:在t時(shí)刻格點(diǎn)7或者格點(diǎn)11被東向的粒子占據(jù),該粒子在t+1時(shí)刻進(jìn)入格點(diǎn)3,概率為2ρ1qs;在t時(shí)刻格點(diǎn)7或者格點(diǎn)11被西向粒子占據(jù),該粒子在t+1時(shí)刻進(jìn)入格點(diǎn)3,概率為2ρ2qs。因此格點(diǎn)3在t和t+1時(shí)刻都為空的概率為1-2ρ1qs-2ρ2qs。

其他情況以此類推。在2-cluster方法中,一共需要求解9個(gè)未知數(shù),這9個(gè)未知數(shù)分別為p(0,0)、p(0,1)、p(0,2)、p(1,0)、p(1,1)、p(1,2)、p(2,0)、p(2,1)、p(2,2)。

通過系統(tǒng)中粒子的密度關(guān)系可得:

p(1,0)+p(1,1)+p(1,2)+p(2,0)+

p(2,1)+p(2,2)=ρ

(2)

p(0,0)+p(0,1)+p(0,2)=1-ρ

(3)

p(1,0)+p(1,1)+p(1,2)=

p(2,0)+p(2,1)+p(2,2)

(4)

模型中東向和西向粒子數(shù)量相等,根據(jù)對稱性可得:

p(1,1)=p(2,2)

(5)

p(1,0)=p(0,2)

(6)

p(2,0)=p(0,1)

(7)

根據(jù)p(2,1)、p(2,0)粒子之間較弱的相關(guān)性假設(shè),忽略相關(guān)性,即

p(2,1)=ρ2/4

(8)

p(2,0)=ρ(1-ρ)/2

(9)

將圖4中所有子圖的54種情況整合即可得到求解p(1,0)的總方程,即

p(1,0)=[p(0,0,0,0)+p(0,0,0,1)+p(2,0,0,0)+p(2,0,0,1)]2ρ1qs(1-2ρ1qs-2ρ2qs)+

[p(0,0,0,2)+p(2,0,0,2)]2ρ1qs(1-2ρ1qs-2ρ2qs-q)+

[p(1,0,0,0)+p(1,0,0,1)]q(1-2ρ1qs-2ρ2qs)+p(1,0,0,2)q(1-2ρ1qs-2ρ2qs-q)+

[p(0,1,0,0)+p(0,1,0,1)+p(1,1,0,0)+p(1,1,0,1)+p(2,1,0,0)+

p(2,1,0,1)][2ρ(1-ρ)qs+ρ2(1-q)](1-2ρ1qs-2ρ2qs)+

[p(0,1,0,2)+p(1,1,0,2)+p(2,1,0,2][2ρ(1-ρ)qs+

ρ2(1-q)](1-2ρ1qs-2ρ2qs-q)+

[p(0,0,2,0)+p(0,0,2,1)+p(0,0,2,2)+p(2,0,2,0)+p(2,0,2,1)+

p(2,0,2,2)](2ρ1qs)[2(1-ρ)qs]+

[p(1,0,2,0)+p(1,0,2,1)+p(1,0,2,2)]q[2(1-ρ)qs]+

[p(0,0,1,0)+p(2,0,1,0)](2ρ1qs)[ρ2q(1-q)+(1/2)ρ2q2+(1-ρ2)q]+

[p(0,0,1,1)+p(0,0,1,2)+p(2,0,1,1)+p(2,0,1,2)](2ρ1qs)[2(1-ρ)qs]+

p(1,0,1,0)q[ρ2q(1-q)+(1/2)ρ2q2+(1-ρ2)q]+

[p(1,0,1,1)+p(1,0,1,2)]q[2(1-ρ)qs]+

[p(0,1,1,0)+p(1,1,1,0)+p(2,1,1,0)][1-2(1-ρ)qs][ρ2q(1-q)+(1/2)ρ2q2+(1-ρ2)q]+

[p(0,1,1,1)+p(0,1,1,2)+p(1,1,1,1)+p(1,1,1,2)+p(2,1,1,1)+

p(2,1,1,2)][1-2(1-ρ)qs][2(1-ρ)qs]+

[p(0,1,2,0)+p(1,1,2,0)+p(2,1,2,0)+p(0,1,2,1)+p(1,1,2,1)+p(2,1,2,1)+

p(0,1,2,2)+p(1,1,2,2)+p(2,1,2,2)][1-2(1-ρ)qs][2(1-ρ)qs]

(10)

聯(lián)立方程(2)～方程(10)可以解得9個(gè)未知數(shù),得到p(1,0)計(jì)算模型中東向粒子的流率為:

J(ρ)=p(1,0)q1

(11)

其中,q1為p(1,0)中的東向粒子(1)在下一個(gè)時(shí)間步往目標(biāo)方向東側(cè)空格(0)運(yùn)動(dòng)的概率。q1的表達(dá)式為:

q1=[ρ23+2ρ1ρ22+ρ12ρ2][q(1-q)(1-qs)2+2(1/2)qqs(1-q)(1-qs)+

q(1-q)(1-qs)][(1/3)q2qs+(1/2)q2(1-qs)+(1/2)qqs(1-q)+q(1-q)(1-qs)]+

[2(1-ρ)ρ1(1-ρ2)+2(1-ρ)ρ2(1-ρ2)][q(1-qs)+(1/2)qqs]+

(1-ρ)2ρ2[q(1-q)+(1/2)q2]+(1-ρ)(1-ρ2)2q

(12)

為得到較為準(zhǔn)確的速度,q1的計(jì)算要考慮p(1,0)中東向粒子(1)往目標(biāo)格點(diǎn)(0)運(yùn)動(dòng)可能發(fā)生的競爭情況。當(dāng)目標(biāo)格點(diǎn)(0)南側(cè)、北側(cè)格點(diǎn)被任一種粒子占據(jù),或目標(biāo)格點(diǎn)(0)東側(cè)格點(diǎn)被目標(biāo)方向?yàn)槲飨虻牧Ｗ诱紦?jù)時(shí),就可能與目標(biāo)粒子(1)發(fā)生競爭,因此方程(12)考慮了所有競爭的情況。

通過東向、西向粒子的對稱性,解得系統(tǒng)中無位置交換部分的粒子平均運(yùn)動(dòng)速度為:

v(ρ)=2J/ρ

(13)

因?yàn)榘l(fā)生位置交換是2個(gè)粒子在目標(biāo)方向上的運(yùn)動(dòng),所以通過p(1,2)的位形可以計(jì)算位置交換的粒子流率,之前2-cluster方法已經(jīng)得到p(1,2)的值,此時(shí)系統(tǒng)中所有粒子的流率計(jì)算方程為:

J1(ρ)=2J(ρ)+2p(1,2)ρ4q2pe

(14)

其中:ρ4為2個(gè)對向粒子的南北兩側(cè)相鄰格點(diǎn)都被其他粒子占據(jù)的概率;q2為2個(gè)對向的粒子朝著目標(biāo)方向運(yùn)動(dòng);pe為位置成功交換的概率。通過對稱性,得到粒子在位置交換時(shí)的流率。

系統(tǒng)中所有粒子的平均速度為:

v1(ρ)=J1/ρ

(15)

聯(lián)立方程(2)～方程(15)得到存在位置交換規(guī)則下的解析結(jié)果(圖2)。

從圖2可以看出,解析結(jié)果與模擬結(jié)果較為吻合，隨著密度的增加,解析與仿真之間的誤差先增加再減小。這是由于當(dāng)密度較小時(shí),系統(tǒng)中實(shí)際發(fā)生位置交換的粒子較少,對粒子平均速度的影響小;隨著密度的繼續(xù)增加,粒子之間的相關(guān)性越來越強(qiáng),而位置交換是影響速度大小的主要因素,2-cluster方法恰好考慮了該相關(guān)性,使求得的p(1,2)與實(shí)際值較為接近,從而使得解析與仿真結(jié)果較為吻合。

3 結(jié) 論

本文以元胞自動(dòng)機(jī)為基礎(chǔ)建立周期邊界的對向粒子流模型,在模型中存在2種目標(biāo)方向不同的粒子,每個(gè)粒子有3個(gè)運(yùn)動(dòng)方向、1個(gè)目標(biāo)方向與2個(gè)偏移方向,采用并行更新的規(guī)則。為了反映生活中粒子運(yùn)動(dòng)發(fā)生沖突的現(xiàn)象,在該模型中加入了一個(gè)位置交換規(guī)則。首先本文進(jìn)行計(jì)算機(jī)模擬,然后通過2-cluster的平均場方法將二維問題劃分成為2個(gè)一維問題,對是否存在位置交換的對向粒子流模型進(jìn)行解析,并在該方法中考慮粒子之間的相關(guān)性。

在該模型中,模擬和解析均發(fā)現(xiàn)粒子的平均速度隨著密度增加不會(huì)下降至0而發(fā)生完全堵塞;隨著q或者pe的增大,則出現(xiàn)了平均速度隨著密度增加先下降后上升的現(xiàn)象,因?yàn)樵诘兔芏葧r(shí),系統(tǒng)中粒子數(shù)對平均速度產(chǎn)生較大的影響,而發(fā)生位置交換的粒子較少,所以位置交換對粒子的平均速度影響較小;而高密度下發(fā)生位置交換的粒子增多,成為影響速度的主要因素,從而使得粒子的平均速度增加。

2-cluster平均場方法在解析均勻態(tài)的系統(tǒng)時(shí)結(jié)果與仿真結(jié)果較為吻合,而出現(xiàn)堵塞帶時(shí)誤差較大,這個(gè)現(xiàn)象可能是粒子間相關(guān)性造成的,處于堵塞帶的粒子相關(guān)性較強(qiáng),同時(shí)阻礙了其他粒子的運(yùn)動(dòng),導(dǎo)致解析與仿真誤差增大,在未來的研究中有待進(jìn)一步探索解析堵塞時(shí)更加精確的方法。