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

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

樓鑫鑫,丁中俊,田波(1.合肥工業(yè)大學(xué)汽車與交通工程學(xué)院,安徽合肥230009;2.安徽農(nóng)業(yè)大學(xué)工學(xué)院,安徽合肥230036)驅(qū)動擴(kuò)散系統(tǒng)因其豐富而復(fù)雜的現(xiàn)象而備受關(guān)注[1]。相分離是在驅(qū)動擴(kuò)散系統(tǒng)[2]中觀察到的一種重要現(xiàn)象,在許多一維或準(zhǔn)一維模型中都有提及,如ABC模型[3]、具有慢啟動規(guī)則的交通流模型[4]、雙車道模型[5]、一維硬粒子系統(tǒng)[6]等。許多研究集中在二維驅(qū)動的擴(kuò)散系統(tǒng),文獻(xiàn)[7]在二維周期網(wǎng)絡(luò)上研究了豐富的相圖,其中2種粒子以相反的方向運(yùn)動;文獻(xiàn)[8]在二維非對稱簡單排它過程中觀察到2種干擾相;文獻(xiàn)[9]在二維系統(tǒng)中研究了流體力學(xué)對液體混合物相分離后期動力學(xué)的影響;文獻(xiàn)[10]在二維Lennard-Jones流體中分析了非平衡通道的形成,該流體由2個相反方向運(yùn)動的顆粒組成;文獻(xiàn)[11]提出2種膠體顆粒的二維模型,表明在驅(qū)動膠體系統(tǒng)中的通道形成,是宏觀范圍內(nèi)的相變,但從無序的初始條件開始,在有限時間內(nèi)不會發(fā)生宏觀相分離;文獻(xiàn)[12]的二維格子氣模型中存在2種粒子,在外力的反向驅(qū)動下研究了“云”的粗化和動態(tài)標(biāo)度;二維模型大多采用隨機(jī)更新規(guī)則,考慮到一些應(yīng)用背景(如行人對向流),文獻(xiàn)[13]建立了一個具有并行更新規(guī)則與存在位置交換的二維驅(qū)動擴(kuò)散系統(tǒng),在該模型中,2種類型的粒子在二維網(wǎng)絡(luò)上以相反的方向運(yùn)動,當(dāng)粒子發(fā)生沖突時根據(jù)周圍環(huán)境判斷是否發(fā)生位置交換,然后采用簡單平均場方法對存在位置交換的二維驅(qū)動擴(kuò)散系統(tǒng)進(jìn)行了解析,但該方法并沒有考慮粒子之間的相關(guān)性而使得誤差較大。本文以文獻(xiàn)[13]為基礎(chǔ),建立存在位置交換的對向粒子流模型。然而粒子在實(shí)際的運(yùn)動過程中存在相互作用,學(xué)者們將相關(guān)性[14]引入到簡單平均場理論,得到簇平均場理論[15],也稱作cluster平均場方法。本文采用簇平均場方法對模型進(jìn)行解析,在解析中考慮粒子間的相關(guān)性,解析結(jié)果與仿真結(jié)果在均勻態(tài)下較為吻合,當(dāng)出現(xiàn)堵塞帶時誤差仍然較大。1對向粒子流模型本文基于元胞自動機(jī)建立對向粒子流模型,如圖1所示。在一個網(wǎng)絡(luò)中有N×N個格點(diǎn),每個格點(diǎn)被一個粒子占據(jù)或者為空。系統(tǒng)中存在2類粒子,每個粒子可以向3個方向運(yùn)動(圖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ī)則為并行更新,即每一個時間步,所有粒子同時進(jìn)行更新;系統(tǒng)的邊界為周期性邊界,即系統(tǒng)中粒子數(shù)目保持不變。每個粒子先由其運(yùn)動概率的大小來確定運(yùn)動方向,再判斷粒子運(yùn)動的目標(biāo)格點(diǎn)是否為空,若格點(diǎn)為空且同時沒有其他粒子試圖進(jìn)入該格點(diǎn),則該粒子在下一個時間步進(jìn)入該格點(diǎn);在同一時間步中,若有n個粒子可以進(jìn)入同一個空格點(diǎn),則以1/n的概率隨機(jī)選擇一個粒子進(jìn)入該格點(diǎn)。此外,模型中還存在位置交換的規(guī)則,圖1b中:在t時刻,若格點(diǎn)3的粒子PE試圖向東運(yùn)動,格點(diǎn)4的粒子PW試圖向西運(yùn)動,且1、2、5、6都被占據(jù)(含對角線的圓表示該格點(diǎn)可能被PE也可能被PW占據(jù)),則在t+1時間步這2個粒子以pe概率發(fā)生位置的交換;若1、2、5、6中有格點(diǎn)為空,則格點(diǎn)3、4的粒子不會發(fā)生位置交換。圖1b的右圖表示t+1時刻3、4格點(diǎn)的2個粒子成功發(fā)生交換,?表示1、2、5、6的狀態(tài)不確定。圖1對向粒子流模型2結(jié)果和討論2.1模擬仿真本文在一個尺寸為100×100的周期邊界網(wǎng)格上進(jìn)行仿真,即系統(tǒng)中存在100×100個格子。首先舍去前107個時間步使系統(tǒng)達(dá)到穩(wěn)定,然后對后面的106個時間步進(jìn)行統(tǒng)計(jì),對每一種情況進(jìn)行20次仿真并計(jì)算平均值,得到考慮位置交換下粒子平均速度與密度的關(guān)系,如圖2所示。圖2中:平均速度為粒子在每個時間步往目標(biāo)方向運(yùn)動的格點(diǎn)數(shù);密度為系統(tǒng)中粒子數(shù)量占比。從圖2可以看出,系統(tǒng)在任何密度下未出現(xiàn)完全堵塞,這是由于平均速度沒有降至0。圖2a中,粒子平均速度隨著密度增加而一直下降,且系統(tǒng)中存在一個臨界密度,當(dāng)粒子密度超過這個值時,粒子的平均運(yùn)動速度下降得更快。圖2b～圖2d中,平均速度隨密度增加而下降,達(dá)到臨界密度之后,平均速度增加,這是由于q值與pe值較大,當(dāng)粒子數(shù)目較多時,相遇粒子之間的互換愈發(fā)頻繁,使粒子可以運(yùn)動。相比于圖2b～圖2d,圖2a中粒子平均速度隨著密度增加持續(xù)下降，這是由于此時q值與pe值較小,相遇的粒子之間發(fā)生位置互換的概率較小,從而相遇的粒子可以運(yùn)動的概率也較小;當(dāng)密度小于臨界密度時,粒子處于均勻的狀態(tài),當(dāng)系統(tǒng)中粒子較多時,系統(tǒng)中產(chǎn)生了一個堵塞帶,阻止了系統(tǒng)運(yùn)動,如圖3所示。圖2在考慮位置交換的情況下粒子運(yùn)動平均速度與密度關(guān)系圖3q=0.6、pe=0.3的位形2.22-cluster平均場解析法本文使用2-cluster平均場理論方法解析對向粒子流模型,并分析粒子流的特性。在cluster平均場理論中,連續(xù)的n個格點(diǎn)即為n-cluster,相鄰的2個格點(diǎn)即為1個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時,t+1時刻2-cluster(σi,σi+1)狀態(tài)取決于t時刻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時,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時刻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時刻的位形;第2列為粒子在t+1時刻的位形;第3列與第4列的乘積表示相應(yīng)的轉(zhuǎn)化概率。在系統(tǒng)中隨機(jī)選擇1個格點(diǎn),格點(diǎn)為空的概率為1-ρ,被占據(jù)的概率為ρ。三角形和圓圈分別對應(yīng)(σi,σi+1)中的σi與σi+1格點(diǎn),且狀態(tài)為空。箭頭表示粒子運(yùn)動的目標(biāo)方向,→表示粒子目標(biāo)方向?yàn)闁|向,←表示粒子目標(biāo)方向?yàn)槲飨??表示格點(diǎn)可能為空也可能被占據(jù);ρ1表示東向粒子的密度,ρ2表示西向粒子的密度,2類粒子數(shù)相同時ρ1=ρ2;qs為粒子向其他方向運(yùn)動的概率。圖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時刻1—4格點(diǎn)為空,5—12格點(diǎn)情況未知,可能被占據(jù)也可能為空;在t+1時刻有1個東向的粒子進(jìn)入格點(diǎn)2,格點(diǎn)3依舊為空,即本文所求的p(1,0)。在模型中偏移概率的影響較小,為了解析的方便與清晰,忽略部分轉(zhuǎn)移概率中偏移方向粒子的影響。圖5圖4中子圖位形演化過程示意圖2ρ1qs表示格點(diǎn)2在t時刻為空,在t+1時刻被東向粒子占據(jù)的概率。在t時刻只有格點(diǎn)6或者格點(diǎn)10的東向粒子能在t+1時刻進(jìn)入格點(diǎn)2號。格點(diǎn)6或者格點(diǎn)10被東向粒子占據(jù)的概率為ρ1,該粒子在t+1時刻進(jìn)入格點(diǎn)2的概率為qs。因此格點(diǎn)2在t+1時刻被東向粒子占據(jù)概率為2ρ1qs。1-2ρ1qs-2ρ2qs表示格點(diǎn)3在t時刻與t+1時刻都為空的概率。因?yàn)樵趖時刻只有格點(diǎn)7或者格點(diǎn)11的粒子能在t+1時刻進(jìn)入格點(diǎn)3，所以存在2種情況:在t時刻格點(diǎn)7或者格點(diǎn)11被東向的粒子占據(jù),該粒子在t+1時刻進(jìn)入格點(diǎn)3,概率為2ρ1qs;在t時刻格點(diǎn)7或者格點(diǎn)11被西向粒子占據(jù),該粒子在t+1時刻進(jìn)入格點(diǎn)3,概率為2ρ2qs。因此格點(diǎn)3在t和t+1時刻都為空的概率為1-2ρ1qs-2ρ2qs。其他情況以此類推。在2-cluster方法中,一共需要求解9個未知數(shù),這9個未知數(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個未知數(shù),得到p(1,0)計(jì)算模型中東向粒子的流率為:J(ρ)=p(1,0)q1(11)其中,q1為p(1,0)中的東向粒子(1)在下一個時間步往目標(biāo)方向東側(cè)空格(0)運(yùn)動的概率。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)動可能發(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ù)時,就可能與目標(biāo)粒子(1)發(fā)生競爭,因此方程(12)考慮了所有競爭的情況。通過東向、西向粒子的對稱性,解得系統(tǒng)中無位置交換部分的粒子平均運(yùn)動速度為:v(ρ)=2J/ρ(13)因?yàn)榘l(fā)生位置交換是2個粒子在目標(biāo)方向上的運(yùn)動,所以通過p(1,2)的位形可以計(jì)算位置交換的粒子流率,之前2-cluster方法已經(jīng)得到p(1,2)的值,此時系統(tǒng)中所有粒子的流率計(jì)算方程為:J1(ρ)=2J(ρ)+2p(1,2)ρ4q2pe(14)其中:ρ4為2個對向粒子的南北兩側(cè)相鄰格點(diǎn)都被其他粒子占據(jù)的概率;q2為2個對向的粒子朝著目標(biāo)方向運(yùn)動;pe為位置成功交換的概率。通過對稱性,得到粒子在位置交換時的流率。系統(tǒng)中所有粒子的平均速度為:v1(ρ)=J1/ρ(15)聯(lián)立方程(2)～方程(15)得到存在位置交換規(guī)則下的解析結(jié)果(圖2)。從圖2可以看出,解析結(jié)果與模擬結(jié)果較為吻合，隨著密度的增加,解析與仿真之間的誤差先增加再減小。這是由于當(dāng)密度較小時,系統(tǒng)中實(shí)際發(fā)生位置交換的粒子較少