圓盤體系在各向同性壓縮下的接觸結構及晶格化

戴丹, 張興剛

(1. 貴州大學計算機科學與技術學院, 貴陽 550025; 2. 貴州大學物理學院, 貴陽 550025)

現采用離散元法(discrete element method, DEM)對雙分散的圓盤體系在各向同性壓縮過程中的堵塞轉變進行模擬。通過調節兩種顆粒的粒徑比使體系成為單分散或近似單分散的體系,研究摩擦與粒徑比對壓強-堆積分數變化曲線的影響。用鄰接矩陣描述圓盤體系中接觸網絡的拓撲結構,采用各種接觸類型的百分比隨堆積分數的變化表征壓縮過程中接觸結構變化的統計特征。著重討論了摩擦、單分散程度等因素對接觸結構變化曲線以及體系晶格化的影響。

1 DEM模型與臨界堵塞態

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

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

(1)

(2)

(3)

在DEM模擬中,隨著半徑R減小,堆積分數增加,體系從疏松態轉變為堵塞態。在此過程中,可以獲得具有不同堆積分數的一組靜態堆積樣本,用于進一步的分析。

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

表1 DEM模擬的主要參數Table 1 Major parameters of DEM simulations

對于散體系統,一個重要的幾何參量是堆積分數φ,它表征體系堆積的密集程度。在二維圓盤堆積中,它定義[10]為

(4)

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

壓強是表征各向同性體系力學狀態的一個關鍵物理量。壓強P可通過計算體系的平均應力張量的跡得到;在圓盤體系中,它的計算式[8,12]為

(5)

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

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

圖1 單分散體系中壓強與堆積分數的關系Fig.1 Pressure vs packing fraction for monodisperse system

本文中主要考慮圓盤體系從臨界堵塞態開始,被施以各向同性壓縮的過程。這時,可以采用接觸網描述體系的堆積結構;接觸網中的結點是圓盤的中心,只有相互接觸的兩個圓盤才被認為是相近鄰的,并用直線段把它們的中心連起來。接觸網被廣泛地用于研究堵塞態顆粒體系的堆積結構和力學性質[3,16,26]。采用接觸網的原因有:①接觸網的產生比其他網絡(如Voronoi網)要簡單很多;②可以通過接觸網把體系的堆積結構與力學性質聯系在一起。可以用圖論中的鄰接矩陣或關聯矩陣對接觸網的拓撲結構進行數學描述。設接觸網的鄰接矩陣為A,將它定義[16]為

(6)

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

圖2 臨界堵塞態的堆積結構Fig.2 The packing structure of disk systems at the critical jammed state

2 摩擦系數對結構變化的影響

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

圖3 不同摩擦系數下壓強與堆積分數的關系Fig.3 Pressure vs packing fraction for different friction

接下來,從臨界堵塞態出發,對圓盤體系進行各向同性地壓縮,并且研究摩擦對壓縮過程中接觸結構變化的影響。在接觸網中,將圓盤周圍與之接觸的圓盤個數稱為該圓盤的接觸數,它可以粗略地刻畫局部接觸結構的穩定性。一般而言,接觸數越大的局部接觸結構越穩定。在網絡科學的術語中,接觸數就是接觸網絡中節點的度。可以根據每個圓盤的接觸數對圓盤進行分類;將接觸數為Z的圓盤稱為Z型圓盤。設NZ是Z型圓盤的總數,那么可以使用接觸類型的百分比QZ反映各種接觸類型在堆積中所占的比率,即

QZ=NZ/N

(7)

圖4所示為摩擦系數μ取不同數值時的QZ-φ變化曲線。首先討論在各種摩擦系數下都有的規律特征。可以看到,0型、1型、2型圓盤的百分比始終都很小,并且隨著體系的壓縮,它們的百分比都在逐漸地減小為0;原因是進行壓縮要求體系能承受更大的壓強,但這三種類型的結構是不平衡或不穩定的,難以承受較大的壓強。從圖4中還可以看出,3型圓盤的百分比隨著φ的增加單調地減小,5型、6型圓盤的百分比隨著φ的增加幾乎在單調地增大,4型圓盤的QZ-φ曲線會隨著摩擦系數的增大由單調遞減的模式變為先增加再減小的模式。這體現了3、4、5、6這幾種類型的圓盤在壓縮過程中進行著比較復雜的轉換。其中3型圓盤以轉出為主,因此是單調遞減;6型圓盤以轉入為主,因此是單調遞增;4型、5型圓盤的變化曲線應該轉入與轉出的競爭關系共同確定。另外,對于各種摩擦系數都有Q3+Q5≈0.5的關系。利用這個關系可知,Q3的單調遞減會使得Q5單調遞增。

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

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

3 粒徑比對結構變化的影響

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

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

圖5 不同粒徑比下壓強與堆積分數的關系Fig.5 Pressure vs packing fraction for different particle size ratio

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

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

4 摩擦及單分散程度對晶格結構的影響

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

圖7 6型顆粒的百分比隨堆積分數增量的變化Fig.7 The percentage of particles for 6 contacts with excess packing fraction

5 結論

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