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