圓盤體系在各向同性壓縮下的接觸結(jié)構(gòu)及晶格化

戴丹, 張興剛

(1. 貴州大學(xué)計算機(jī)科學(xué)與技術(shù)學(xué)院, 貴陽 550025; 2. 貴州大學(xué)物理學(xué)院, 貴陽 550025)

現(xiàn)采用離散元法(discrete element method, DEM)對雙分散的圓盤體系在各向同性壓縮過程中的堵塞轉(zhuǎn)變進(jìn)行模擬。通過調(diào)節(jié)兩種顆粒的粒徑比使體系成為單分散或近似單分散的體系,研究摩擦與粒徑比對壓強(qiáng)-堆積分?jǐn)?shù)變化曲線的影響。用鄰接矩陣描述圓盤體系中接觸網(wǎng)絡(luò)的拓?fù)浣Y(jié)構(gòu),采用各種接觸類型的百分比隨堆積分?jǐn)?shù)的變化表征壓縮過程中接觸結(jié)構(gòu)變化的統(tǒng)計特征。著重討論了摩擦、單分散程度等因素對接觸結(jié)構(gòu)變化曲線以及體系晶格化的影響。

1 DEM模型與臨界堵塞態(tài)

離散元法是一種類似于分子動力學(xué)的數(shù)值模擬方法,它被廣泛地用于研究顆粒物質(zhì)的幾何性質(zhì)與力學(xué)行為[24-25]。離散元法的基本思想是:首先將不連續(xù)的離散介質(zhì)看成一組離散單元的集合,每一個單元可作平動和轉(zhuǎn)動,各單元間也可以接觸產(chǎn)生相互作用,單元的運動服從經(jīng)典力學(xué)所給出的動力學(xué)方程;然后將時間離散化并且迭代求解各個離散單元服從的動力學(xué)方程,進(jìn)而得到非連續(xù)介質(zhì)的運動過程及狀態(tài)。文中用DEM模擬圓盤體系在各向同性壓縮下的堵塞轉(zhuǎn)變以及堵塞態(tài)在壓縮過程中的結(jié)構(gòu)變化,具體的模擬方法如下。

(1)在半徑為R0的圓形容器內(nèi)隨機(jī)、均勻地產(chǎn)生N個圓盤,從而生成圓盤體系的一個初始構(gòu)型。圓盤體系包含NB個半徑為rB的大圓盤,NS個半徑為rS的小圓盤。為了方便,引入三個無量綱的量來描述這個圓盤體系:粒徑比λ=rB/rS、大圓盤所占的比率kB=NB/N、小圓盤所占的比率kS=NS/N;顯然,kB+kS=1。體系的初始構(gòu)型為堆積分?jǐn)?shù)很小的疏松態(tài);此時,圓盤隨機(jī)均勻地分布于圓形容器內(nèi),圓盤間不相互重疊。模擬過程中不考慮重力的影響,也不考慮圓盤的熱運動。

(1)

(2)

(3)

在DEM模擬中,隨著半徑R減小,堆積分?jǐn)?shù)增加,體系從疏松態(tài)轉(zhuǎn)變?yōu)槎氯麘B(tài)。在此過程中,可以獲得具有不同堆積分?jǐn)?shù)的一組靜態(tài)堆積樣本,用于進(jìn)一步的分析。

(3)由于隨機(jī)因素的影響,對同一條件下的不同樣本進(jìn)行統(tǒng)計平均才有實際意義。為此,在同樣的模擬參數(shù)下生成許多不同的初始構(gòu)型,然后對各個初始構(gòu)型進(jìn)行各向同性壓縮從而得到多組靜態(tài)堆積樣本。接著,對多組靜態(tài)堆積樣本進(jìn)行統(tǒng)計平均,從而得到一種模擬參數(shù)下的統(tǒng)計結(jié)果。最后,改變模擬的參數(shù),再執(zhí)行上述DEM模擬過程,進(jìn)而分析體系的結(jié)構(gòu)與性質(zhì)如何隨著這些參數(shù)的變化而改變。DEM模擬的主要參數(shù)如表1所示。

表1 DEM模擬的主要參數(shù)Table 1 Major parameters of DEM simulations

對于散體系統(tǒng),一個重要的幾何參量是堆積分?jǐn)?shù)φ,它表征體系堆積的密集程度。在二維圓盤堆積中,它定義[10]為

(4)

式(4)中:An為第n個圓盤的面積;AC為圓形容器的面積。

壓強(qiáng)是表征各向同性體系力學(xué)狀態(tài)的一個關(guān)鍵物理量。壓強(qiáng)P可通過計算體系的平均應(yīng)力張量的跡得到;在圓盤體系中,它的計算式[8,12]為

(5)

式(5)中:S為圓盤體系的總力矩;Zn為第n個圓盤周圍的接觸點數(shù);Fni為第n個圓盤在其第i個接觸點(定義為圓盤間重疊部分的中心)處所受的法向力的大小;rni為該圓盤中心到該接觸點的距離。壓強(qiáng)P也可通過將容器所受的法向力大小的總和除以容器邊長而得到。

使用式(4)和式(5),可以對DEM模擬中得到的靜態(tài)堆積樣本進(jìn)行數(shù)據(jù)處理,并且給出壓縮過程中表征壓強(qiáng)與堆積分?jǐn)?shù)關(guān)系的P-φ曲線。圖1是粒徑比λ=1(此時為單分散體系)、摩擦系數(shù)μ=0、法向剛度系數(shù)Kn=1×104N/m、圓盤數(shù)N取不同數(shù)值時體系的P-φ變化曲線。從圖1中可以看到,隨著堆積分?jǐn)?shù)的增大,體系由疏松態(tài)轉(zhuǎn)變?yōu)槎氯麘B(tài),J點就是臨界轉(zhuǎn)變點,對應(yīng)的狀態(tài)為臨界堵塞態(tài),相應(yīng)的堆積分?jǐn)?shù)用φJ(rèn)表示。可以看到:當(dāng)N=500時φJ(rèn)≈0.825;N=2 500時φJ(rèn)≈0.845;N=5 000時φJ(rèn)≈0.855;N=6 000時也有φJ(rèn)≈0.855;隨著圓盤數(shù)的增多,φJ(rèn)逐漸增大,圓盤數(shù)足夠多時φJ(rèn)幾乎不再變。這是因為圓盤數(shù)較少時,由于受到邊界形狀的影響,圓盤不易在其內(nèi)部形成較致密的堆積;圓盤數(shù)較多時,容器內(nèi)部遠(yuǎn)離邊界的圓盤容易形成較致密的堆積;圓盤數(shù)足夠多時,邊界形狀的影響可以忽略不計,這使φJ(rèn)幾乎不再變化。為了避免邊界的影響,文中主要對N足夠大的情形(N≥5 000)進(jìn)行模擬。進(jìn)入堵塞態(tài)后,隨著堆積分?jǐn)?shù)的增加,壓強(qiáng)幾乎呈線性增加。改變剛度系數(shù)重新進(jìn)行模擬,也會得到類似的P-φ曲線。只不過剛度系數(shù)的增大會增加壓強(qiáng)線性變化時的斜率。

圖1 單分散體系中壓強(qiáng)與堆積分?jǐn)?shù)的關(guān)系Fig.1 Pressure vs packing fraction for monodisperse system

本文中主要考慮圓盤體系從臨界堵塞態(tài)開始,被施以各向同性壓縮的過程。這時,可以采用接觸網(wǎng)描述體系的堆積結(jié)構(gòu);接觸網(wǎng)中的結(jié)點是圓盤的中心,只有相互接觸的兩個圓盤才被認(rèn)為是相近鄰的,并用直線段把它們的中心連起來。接觸網(wǎng)被廣泛地用于研究堵塞態(tài)顆粒體系的堆積結(jié)構(gòu)和力學(xué)性質(zhì)[3,16,26]。采用接觸網(wǎng)的原因有:①接觸網(wǎng)的產(chǎn)生比其他網(wǎng)絡(luò)(如Voronoi網(wǎng))要簡單很多;②可以通過接觸網(wǎng)把體系的堆積結(jié)構(gòu)與力學(xué)性質(zhì)聯(lián)系在一起?？梢杂脠D論中的鄰接矩陣或關(guān)聯(lián)矩陣對接觸網(wǎng)的拓?fù)浣Y(jié)構(gòu)進(jìn)行數(shù)學(xué)描述。設(shè)接觸網(wǎng)的鄰接矩陣為A,將它定義[16]為

(6)

在DEM模擬的數(shù)據(jù)中,每個圓盤的位矢、半徑是已知的,利用這些數(shù)據(jù)可以畫出圓盤的堆積情況,同時也可以畫出圓盤體系的接觸網(wǎng)。圖2中給出了單分散體系和雙分散體系在臨界堵塞態(tài)的堆積結(jié)構(gòu)?？梢钥吹?圖2(a)所示的單分散體系由很多個三角或四角晶格團(tuán)簇組成。三角晶格中每個圓盤周圍有六個近鄰圓盤,四角晶格中每個圓盤周圍有四個近鄰圓盤。不過由于初始隨機(jī)堆積的條件和圓形邊界條件的影響,單分散體系也很難完全晶格化,體系中總存在一些晶格裂紋。而圖2(b)所示的雙分散體系則是無序的結(jié)構(gòu)?？梢钥吹?接觸網(wǎng)能夠直觀地體現(xiàn)圓盤體系的堆積結(jié)構(gòu)。

圖2 臨界堵塞態(tài)的堆積結(jié)構(gòu)Fig.2 The packing structure of disk systems at the critical jammed state

2 摩擦系數(shù)對結(jié)構(gòu)變化的影響

在單分散圓盤體系中,通過改變摩擦系數(shù)μ來討論其對臨界堵塞態(tài)以及接觸結(jié)構(gòu)的影響。圖3所示為不同摩擦系數(shù)μ下的P-φ曲線。從圖3中可以看出,當(dāng)堆積分?jǐn)?shù)φ<φJ(rèn)(μ)時,體系處于非堵塞態(tài);當(dāng)φ>φJ(rèn)(μ)時體系處于堵塞態(tài)。圓盤體系在臨界堵塞態(tài)處的堆積分?jǐn)?shù)φJ(rèn)(μ)隨著摩擦系數(shù)μ的增大而逐漸減小;這是由于摩擦系數(shù)的增大會增加接觸網(wǎng)絡(luò)的穩(wěn)定性,從而使得較疏松的堆積也能開始承受外部載荷。從臨界堵塞態(tài)開始,壓強(qiáng)與堆積分?jǐn)?shù)近似呈線性增長關(guān)系。不同的摩擦系數(shù)對于進(jìn)入堵塞態(tài)后的壓縮過程也有一定影響,因為將不同摩擦系數(shù)的P-φ曲線平移后一般不能完全重合。顯然,μ=5和μ=10時的P-φ曲線很好地重合在一起,這說明當(dāng)摩擦系數(shù)足夠大后,摩擦系數(shù)的增加很難進(jìn)一步地提升接觸網(wǎng)絡(luò)的穩(wěn)定性。

圖3 不同摩擦系數(shù)下壓強(qiáng)與堆積分?jǐn)?shù)的關(guān)系Fig.3 Pressure vs packing fraction for different friction

接下來,從臨界堵塞態(tài)出發(fā),對圓盤體系進(jìn)行各向同性地壓縮,并且研究摩擦對壓縮過程中接觸結(jié)構(gòu)變化的影響。在接觸網(wǎng)中,將圓盤周圍與之接觸的圓盤個數(shù)稱為該圓盤的接觸數(shù),它可以粗略地刻畫局部接觸結(jié)構(gòu)的穩(wěn)定性。一般而言,接觸數(shù)越大的局部接觸結(jié)構(gòu)越穩(wěn)定。在網(wǎng)絡(luò)科學(xué)的術(shù)語中,接觸數(shù)就是接觸網(wǎng)絡(luò)中節(jié)點的度?？梢愿鶕?jù)每個圓盤的接觸數(shù)對圓盤進(jìn)行分類;將接觸數(shù)為Z的圓盤稱為Z型圓盤。設(shè)NZ是Z型圓盤的總數(shù),那么可以使用接觸類型的百分比QZ反映各種接觸類型在堆積中所占的比率,即

QZ=NZ/N

(7)

圖4所示為摩擦系數(shù)μ取不同數(shù)值時的QZ-φ變化曲線。首先討論在各種摩擦系數(shù)下都有的規(guī)律特征?？梢钥吹?0型、1型、2型圓盤的百分比始終都很小,并且隨著體系的壓縮,它們的百分比都在逐漸地減小為0;原因是進(jìn)行壓縮要求體系能承受更大的壓強(qiáng),但這三種類型的結(jié)構(gòu)是不平衡或不穩(wěn)定的,難以承受較大的壓強(qiáng)。從圖4中還可以看出,3型圓盤的百分比隨著φ的增加單調(diào)地減小,5型、6型圓盤的百分比隨著φ的增加幾乎在單調(diào)地增大,4型圓盤的QZ-φ曲線會隨著摩擦系數(shù)的增大由單調(diào)遞減的模式變?yōu)橄仍黾釉贉p小的模式。這體現(xiàn)了3、4、5、6這幾種類型的圓盤在壓縮過程中進(jìn)行著比較復(fù)雜的轉(zhuǎn)換。其中3型圓盤以轉(zhuǎn)出為主,因此是單調(diào)遞減;6型圓盤以轉(zhuǎn)入為主,因此是單調(diào)遞增;4型、5型圓盤的變化曲線應(yīng)該轉(zhuǎn)入與轉(zhuǎn)出的競爭關(guān)系共同確定。另外,對于各種摩擦系數(shù)都有Q3+Q5≈0.5的關(guān)系。利用這個關(guān)系可知,Q3的單調(diào)遞減會使得Q5單調(diào)遞增。

圖4 不同摩擦系數(shù)下Z型顆粒的百分比Fig.4 The percentage of particles with Z contacts vs packing fraction for different friction

接下來討論摩擦對各種類型圓盤的QZ-φ曲線的影響。從圖4中可看出,摩擦系數(shù)較小(μ<0.1)時,Q6(φ)的增長率比較大;壓縮的持續(xù)會使得6型圓盤的百分比超過其他所有類型圓盤的百分比;這表明摩擦系數(shù)較小時,體系在壓縮的過程中容易形成三角晶格團(tuán)簇。隨著摩擦系數(shù)的增大,Q6(φ)的增長率在減小。當(dāng)摩擦系數(shù)比較大時(μ≥0.8),6型圓盤的百分比始終都很小;這表明大的摩擦系數(shù)下難以形成三角晶格團(tuán)簇。對于4型、5型圓盤,隨著摩擦系數(shù)的增大,Q5(φ)的平均增長率有所增加,Q4(φ)的平均變化率由負(fù)向正轉(zhuǎn)變。另外,可以關(guān)注4、5、6型圓盤在臨界堵塞態(tài)以及深度堵塞態(tài)(定義為堆積分?jǐn)?shù)比較大的堵塞態(tài))處百分比的排序關(guān)系。從圖4還可以看出,臨界堵塞態(tài)在各種摩擦系數(shù)下都有Q4>Q5>Q6,這使得臨界堵塞態(tài)的平均接觸數(shù)比較小。對于深度堵塞態(tài),當(dāng)摩擦系數(shù)很小時,有Q6>Q5>Q4;隨著摩擦系數(shù)的增加,會變?yōu)镼5>Q6>Q4,然后變?yōu)镼5>Q4>Q6,最后變?yōu)镼4>Q5>Q6。這使得深度堵塞態(tài)的平均接觸數(shù)比臨界堵塞態(tài)時大,也使得深度堵塞態(tài)的平均接觸數(shù)隨摩擦系數(shù)的增加而減小。

3 粒徑比對結(jié)構(gòu)變化的影響

在雙分散圓盤體系中,保持摩擦系數(shù)μ=0不變,討論粒徑比λ對臨界堵塞態(tài)以及接觸結(jié)構(gòu)的影響。顯然,大小圓盤的粒徑比λ越接近1,體系的單分散程度就越高,反之就越低。這里主要考慮λ比較接近1,也就是近似單分散的情形。

首先,當(dāng)圓盤體系的粒徑比發(fā)生變化時,它的臨界堵塞態(tài)也會發(fā)生改變。圖5所示為不同粒徑比λ下的P-φ曲線。從圖5中可以看出,圓盤體系在臨界堵塞態(tài)處的堆積分?jǐn)?shù)φJ(rèn)(λ)隨著粒徑比λ的增大而逐漸減小。在進(jìn)入堵塞態(tài)后壓強(qiáng)與堆積分?jǐn)?shù)近似呈線性增長關(guān)系。不同的粒徑比對于進(jìn)入堵塞態(tài)后的壓縮過程也有一定影響,因為將不同粒徑比的P-φ曲線平移后一般不能完全重合。顯然,λ≥1.15后的幾條P-φ曲線基本重合在一起;這是由于λ=1.15時體系處于無序堆積的狀態(tài),在一定范圍內(nèi)增大粒徑比不會更多地改變體系的無序堆積結(jié)構(gòu)。不過如果在λ=1.4的基礎(chǔ)上增大粒徑比,這時小顆?？梢赃M(jìn)入大顆粒之間的縫隙,這一效應(yīng)會增加堆積分?jǐn)?shù)。

圖5 不同粒徑比下壓強(qiáng)與堆積分?jǐn)?shù)的關(guān)系Fig.5 Pressure vs packing fraction for different particle size ratio

圖6為雙分散圓盤體系在粒徑比逐漸增大的情況下Z型圓盤百分比的變化,可知在雙分散圓盤體系中,由于存在一個大圓盤周圍有許多小圓盤的可能性,因此會有7型圓盤或其他類型的圓盤。從圖6中可以看出,1型、2型圓盤的百分比幾乎為0。對于0型圓盤,其百分比隨著堆積分?jǐn)?shù)的增加單調(diào)減小。在粒徑比λ接近1時,Q0(φ)接近于0,但是當(dāng)λ遠(yuǎn)離1時,Q0(φ)的數(shù)值并不是很小,這說明半徑差異越大的雙分散體系越容易形成可以容納小顆粒的孔洞結(jié)構(gòu)。當(dāng)粒徑比λ接近1時,由于無法形成7型圓盤,因此Q7(φ)=0;當(dāng)λ偏離1較大時才可以形成7型圓盤。對于各種粒徑比,Q3(φ)都呈現(xiàn)單調(diào)遞減,Q6(φ)都呈現(xiàn)單調(diào)遞增的變化趨勢。不過對于Q4(φ)與Q5(φ),它們在不同的粒徑比有不同的變化特征。

圖6 不同粒徑比下Z型顆粒的百分比Fig.6 The percentage of particles with Z contacts for different particle size ratio

4 摩擦及單分散程度對晶格結(jié)構(gòu)的影響

接下來,討論摩擦及單分散程度對圓盤體系晶格結(jié)構(gòu)及其演化的影響。這里主要討論圓盤體系中的三角晶格。可以通過統(tǒng)計6型顆粒的百分比Q6來討論圓盤體系的晶格化程度。圖7所示為不同條件下6型圓盤百分比隨堆積分?jǐn)?shù)增量Δφ=φ-φJ(rèn)的變化曲線。可以看到,在同一摩擦系數(shù)或粒徑比的情況下,壓縮過程的進(jìn)行都會提高晶格化的程度。在圖7(a)中,隨著摩擦系數(shù)的減小,同一堆積分?jǐn)?shù)下的Q6值越來越大,Q6(φ)的增長率也在變大。這說明摩擦系數(shù)越小,體系越容易晶格化。當(dāng)摩擦系數(shù)μ很接近0時,6型圓盤的百分比可超過50%,這時體系由許多三角晶格團(tuán)簇構(gòu)成。當(dāng)摩擦系數(shù)μ>0.8之后,很難形成三角晶格的結(jié)構(gòu)。由圖7(b)可以看到,隨著λ從1出發(fā)開始增大,Q6(φ)的曲線在開始時會發(fā)生很大的變化;但是當(dāng)λ>1.05之后,Q6(φ)曲線的變化就很小。這說明圓盤體系的晶格化很容易受到粒徑比的影響。只有單分散程度很高(粒徑比很接近1)時,壓縮過程的持續(xù)才會導(dǎo)致明顯的晶格化。當(dāng)單分散程度較低時,很難有三角晶格結(jié)構(gòu)的形成。

圖7 6型顆粒的百分比隨堆積分?jǐn)?shù)增量的變化Fig.7 The percentage of particles for 6 contacts with excess packing fraction

5 結(jié)論

結(jié)合離散元模擬及接觸網(wǎng)絡(luò)的統(tǒng)計分析對雙分散圓盤體系在各向同性壓縮下的接觸結(jié)構(gòu)及其晶格化進(jìn)行研究。模擬的結(jié)果說明摩擦及粒徑比對體系的臨界阻塞態(tài)會產(chǎn)生重要的影響。隨著圓盤數(shù)的增多,體系的臨界堆積分?jǐn)?shù)φJ(rèn)逐漸增大;隨著摩擦系數(shù)的增大,φJ(rèn)逐漸減小;而隨著粒徑比的增大,φJ(rèn)也逐漸減小。通過接觸類型的百分比QZ更細(xì)致地討論了接觸結(jié)構(gòu)隨摩擦系數(shù)及粒徑比的變化,結(jié)果表明各種接觸類型的百分比隨摩擦系數(shù)和粒徑比有較為復(fù)雜的變化。同時,討論了摩擦和單分散程度對晶格結(jié)構(gòu)的影響。在單分散圓盤體系中,隨著壓縮的不斷進(jìn)行,晶格化程度也不斷增加;當(dāng)摩擦系數(shù)較小(μ<0.1)時,體系在壓縮的過程中容易形成三角晶格團(tuán)簇。隨著摩擦系數(shù)的增大,Q6(φ)的增長率在減小,當(dāng)摩擦系數(shù)比較大(μ>0.8)時則難以形成三角晶格團(tuán)簇。在雙分散圓盤體系中,粒徑比的變化會對晶格結(jié)構(gòu)的形成產(chǎn)生重要的影響。只有單分散程度很高(粒徑比很接近1)時,壓縮過程的持續(xù)才會導(dǎo)致明顯的晶格化。當(dāng)單分散程度較低(λ>1.05)時,很難形成三角晶格結(jié)構(gòu)。