張 祺 喬 婷 季順迎 厚美瑛

* (中國科學(xué)院物理研究所,北京 100080)

? (太原理工大學(xué)機械與運載工程學(xué)院,太原 030024)

** (大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室,遼寧大連 116023)

引言

熱系統(tǒng)中無序的粒子在熱驅(qū)動下進行重新排列,自發(fā)形成有序的結(jié)構(gòu),稱為自組裝過程[1-2].作為非熱學(xué)流,盡管熱運動對顆粒物質(zhì)的影響可以忽略(mgd?kBT),但容器的振動[3-5]、旋轉(zhuǎn)[6]、敲擊[7]、剪切[8-9]以及重力驅(qū)動下自身流動[10]也可以發(fā)揮熱漲落的驅(qū)動作用,改變顆粒體系的構(gòu)型,達到整體或者部分顆粒有序堆積的效果.宏觀顆粒有序堆積過程為研究微觀粒子系統(tǒng)自組裝問題提供了很好的借鑒.

顆粒的有序堆積必然伴隨顆粒系統(tǒng)的密實化過程,這不僅涉及拓撲結(jié)構(gòu)、對稱性等顆粒系統(tǒng)結(jié)構(gòu)優(yōu)化問題,而且對提高交通運輸、包裝工業(yè)、工農(nóng)業(yè)生產(chǎn)中顆粒材料的空間利用率具有重要現(xiàn)實意義[11].研究發(fā)現(xiàn)三維單分散球體的最密隨機堆積(random close packing,RCP)體積分數(shù)約為0.64[12-13],更高體積分數(shù)尤其是最密規(guī)則堆積的三維單分散球體結(jié)構(gòu)則需要一些繁復(fù)的方法才能實現(xiàn).Carvente等[14]利用一種類似退火的方式(振動強度逐漸降低)實現(xiàn)了三維方形容器中的等大鋼球從初始的顆粒氣體狀態(tài)逐漸形成體心四方或者面心立方結(jié)構(gòu).Yu 等[15]通過一個類似晶體外延生長的方式向振動顆粒床中逐步添加球體最終實現(xiàn)了三維的最密規(guī)則堆積結(jié)構(gòu).Shinde 等[16]使用混合蒙特卡羅算法模擬等大球在振動下的自發(fā)結(jié)晶化過程,發(fā)現(xiàn)較低振動強度下結(jié)晶化只會在局部區(qū)域出現(xiàn)而不能實現(xiàn)全局結(jié)晶化,而高振動強度時結(jié)晶化過程表現(xiàn)為六方最密堆積結(jié)構(gòu)和面心立方結(jié)構(gòu)的競爭.然而現(xiàn)實世界中真實顆粒多是非球形的,非球形顆粒與球體顆粒的堆積性質(zhì)存在較大的差異.其中一個主要原因是非球形顆粒接觸形式遠比球形顆粒接觸更加復(fù)雜,包括面-面、邊-面、邊-邊和頂點-面4 種基本接觸形式[17].這促使人們開始進行關(guān)于非球形顆粒堆積性質(zhì)的研究.

Baker 等[18]通過沉降實驗研究了正四面體、正六面體、正八面體、正十二面體和正二十面體共五種柏拉圖實體的最密隨機堆積(RCP)和最松隨機堆積(random loose packing,RLP)的體積分數(shù),發(fā)現(xiàn)除了立方體可以實現(xiàn)無空隙的全空間密堆積外,其他四種實體的RCP 和RLP 體積分數(shù)均隨邊數(shù)的增加而逐漸下降,實體的摩擦系數(shù)越大其RCP 和RLP 的體積分數(shù)越低.Donev 等[19]發(fā)現(xiàn)不同長短軸的橢球顆粒最密隨機堆積體積分數(shù)不盡相同,最高可達0.74,此時橢球平均配位數(shù)遠高于球體.Torquato等[20]通過ASC (adaptive shrinking cell)數(shù)值模擬算法指出正八面體最密規(guī)則堆積的體積分數(shù)約為0.947,正十二面體約為0.904,而正二十面體約為0.836,遠高于實驗發(fā)現(xiàn)的隨機堆積體積分數(shù).B?rzs?nyi 等[21]通過環(huán)剪實驗發(fā)現(xiàn)細長棒狀顆粒會自發(fā)形成一定取向的規(guī)則排列,且同流線之間存在一個特定夾角,該夾角與剪切速率無關(guān).Hidalgo 等[22]通過筒倉中棒狀顆粒沉降實驗發(fā)現(xiàn),顆粒傾向于以不同的傾角在筒倉中形成有序排列,最可幾傾角的大小與顆粒的長徑比有關(guān),其中方形顆粒最易以對角線沿重力方向排列.Asencio 等[23]通過容器往復(fù)水平旋轉(zhuǎn)的方式而不是傳統(tǒng)的振動或者敲擊的方式對圓筒容器中的立方體骰子進行旋轉(zhuǎn)剪切,實現(xiàn)了立方體顆粒最密規(guī)則堆積,為非球形顆粒的堆積研究提供了新思路.

盡管實驗上已經(jīng)對旋轉(zhuǎn)圓筒中立方體顆粒的有序堆積的研究取得了重要進展,然而由于所研究顆粒體系的復(fù)雜性和控制參數(shù)的多樣性,目前對非球形顆粒有序堆積的細致過程和力學(xué)機制的理解還不夠深入,需要進一步研究討論.三維顆粒系統(tǒng)內(nèi)部結(jié)構(gòu)的實驗觀測通常是困難的,其力學(xué)行為的問題也難以通過實驗得到完全解決.離散元方法能夠跟蹤到每一個顆粒的空間位置,從而可以記錄內(nèi)部結(jié)構(gòu)在介觀尺度的細節(jié)信息,是對顆粒體系動力學(xué)行為實驗研究的有力補充[24-26].因此,本文利用離散元方法對旋轉(zhuǎn)圓筒容器內(nèi)立方體顆粒有序堆積過程進行數(shù)值模擬以獲取立方體顆粒有序化過程中的結(jié)構(gòu)演變特征,分析了容器運動方式對立方體顆粒有序堆積過程的影響規(guī)律,有助于深入理解復(fù)雜顆粒系統(tǒng)的結(jié)構(gòu)演化的微觀機理.

1 圓角立方體顆粒有序堆積過程的離散元數(shù)值方法

1.1 基于超二次曲面的圓角立方體顆粒

超二次曲面形狀是二次曲面的延伸,80%的常見非球形顆粒幾何形狀可由超二次曲面方程描述[27-28].本文所模擬的立方體顆粒具有高度的對稱性,因此選用超二次曲面連續(xù)函數(shù)表示.超二次曲面顆粒的形狀表示為

其中,a,b,c分別為顆粒沿主軸方向的半軸長,n1和n2是表征顆粒尖銳度的塊效應(yīng)參數(shù).通過改變式(1)中表征顆粒形狀的五個形狀參數(shù)(a,b,c,n1,n2),可以靈活地獲得不同形狀的顆粒.當(dāng) n1=n2=8,a=b=c 時,即可構(gòu)造帶有較小圓角的立方體顆粒單元.

1.2 圓角立方體顆粒單元間的接觸檢測

顆粒接觸檢測的主要內(nèi)容是檢測顆粒是否接觸以及確定接觸點和接觸法向.超二次曲面單元間的接觸檢測問題可以考慮為一個定義在全局坐標(biāo)系下連續(xù)且二階可導(dǎo)形狀函數(shù)的超二次曲面顆粒間最小距離的最優(yōu)化問題[25,29-30].最優(yōu)化問題通過設(shè)立拉式算子 λ,由下面的拉式函數(shù)給出

其中 X=(x,y,z)T,最優(yōu)化問題的解由?X,λL(X,λ)=0確定,引入 μ2=(1-λ)/(1+λ),最終簡化為四元非線性方程組

確定法向重疊量的牛頓迭代公式可寫為

式中,X=X+dX,μ=μ+dμ .計算結(jié)果 X0滿足Fi(X0)＜0 且 Fj(X0)＜0 則表明兩個單元間發(fā)生接觸.接觸法向 n 可表示為n=?Fi(X0)/||?Fi(X0)|| .進一步考慮Xi=X0+αn 和 Xj=X0+βn,這里未知參數(shù) α 和β 可通過一元非線性牛頓迭代得到:因此,法向重疊量可表述為 δn=Xi-Xj.

1.3 接觸模型及接觸力的計算

基于球體的赫茲接觸模型和橢球體的表面曲率接觸模型,超二次曲面單元的非線性接觸力可按照考慮等效曲率的接觸模型計算[31].在超二次曲面顆粒的相互作用中,顆粒間的相互作用力可認為是法向接觸力和切向接觸力的疊加.法向接觸力主要包括彈性力和黏滯力,此外非球形不規(guī)則顆粒的接觸力可能不通過顆粒形心,因此引起的附加彎矩也應(yīng)被考慮[32].

法向作用力可以表示為

其中,Kn為兩個接觸顆粒間法向有效剛度,δn為兩個接觸顆粒間的法向重疊量,Cn為兩個接觸顆粒間的法向黏性系數(shù),B 為兩個接觸顆粒間的黏滯系數(shù),vn為兩個接觸顆粒間的法向相對速度.

法向有效剛度 Kn和黏滯系數(shù) B 可由顆粒材料的彈性模量 E、泊松比 ν 和接觸點處的平均曲率半徑R 通過下式得到

其中,Ri和 Rj分別表示顆粒 i 和顆粒 j 的等體積球半徑,mi和mj分別表示顆粒 i 和顆粒j 的質(zhì)量.

其中,μt為滑動摩擦系數(shù),δt為切向重疊量,

其中,Ct為切向黏滯系數(shù),vt,ij為單元間的切向相對速度.

當(dāng)顆粒間發(fā)生相對轉(zhuǎn)動時,由滾動摩擦引起的力矩可表示為

其中,μr為滾動摩擦系數(shù),=ω/|ω| 為接觸顆粒間的單位相對角速度.

1.4 圓角立方體顆粒的運動控制方程

顆粒的運動包括平動和轉(zhuǎn)動,其平動的控制方程為

其中,mi和 Xci是顆粒i 的質(zhì)量和顆粒質(zhì)心的坐標(biāo),Fi是作用在顆粒 i 上的合外力.

而顆粒的轉(zhuǎn)動遵循動量矩定理為

其中,Ii和 ωi分別是顆粒 i 全局坐標(biāo)下的慣性矩和角速度,Li是顆粒 i 的角動量,Ti是作用在顆粒 i 上對顆粒質(zhì)心的總轉(zhuǎn)動力矩.對于非球形顆粒,其慣性矩與顆粒的位向有關(guān),在全局坐標(biāo)系下,非球形顆粒的慣性矩是一個各個分量都可能不為零的二階張量,運動的描述十分復(fù)雜,而在固定于顆粒本身的局部坐標(biāo)系中,非球形顆粒的主慣性矩是一個僅對角線元素不為零的二階張量,運動描述簡便.因此,引入可以實現(xiàn)局部坐標(biāo)系和全局坐標(biāo)系相互轉(zhuǎn)化的四元數(shù) q 表示顆粒的位向

非球形顆粒的位向通常由全局坐標(biāo)系 eG向固定于顆粒本身的局部坐標(biāo)系 eL的轉(zhuǎn)換表示,而四元數(shù)q 可以通過轉(zhuǎn)換矩陣 A 實現(xiàn)這種坐標(biāo)系之間的相互轉(zhuǎn)換

這里轉(zhuǎn)換矩陣滿足 A-1=AT,且由四元數(shù)確定

根據(jù)建立在角動量定理上的歐拉運動定律可以得到局部坐標(biāo)系下的顆粒轉(zhuǎn)動方程

1.5 圓角立方體顆粒有序堆積過程的實施方案

程序中顆粒所處容器為圓柱形轉(zhuǎn)筒,半徑r=0.067 m,高度h=0.17 m.初始時,800 個邊長為0.01 m 的立方體顆粒在圓筒上方空間生成,之后以一定速度掉落至容器中并等待直至形成位置和取向均隨機的堆積體.通過邊界的“往復(fù)旋轉(zhuǎn)”向立方體顆粒堆積體傳遞邊界轉(zhuǎn)動產(chǎn)生的剪切作用,如圖1所示.離散元模擬中所需的主要計算參數(shù)如表1 所示,模擬環(huán)境重力加速度為g=9.81 m/s2.轉(zhuǎn)筒邊界以一個選定的峰值角速度 ω0進行正反向的循環(huán)旋轉(zhuǎn),邊界的切向加速度a 為

表1 轉(zhuǎn)動驅(qū)動立方體顆粒有序堆積過程的主要計算參數(shù)Table 1 DEM simulation parameters of ordering of cubes

圖1 裝置示意圖以及邊界旋轉(zhuǎn)方式Fig.1 Schematic diagram of simulation and boundary rotation mode

其中R 為圓柱形邊界的半徑,v 是邊界的切向速度,Δt是反轉(zhuǎn)過程需要的時間.使用重力加速度無量綱化后的扭轉(zhuǎn)加速度 γ 為

旋轉(zhuǎn)過程中單向轉(zhuǎn)動的角位移 β 為

在上述模型和計算條件的基礎(chǔ)上,借助自主開發(fā)的GPU 并行離散元計算程序針對立方體顆粒受轉(zhuǎn)動驅(qū)動有序堆積過程進行模擬計算.計算平臺為配備英特爾xeon 4210 CPU 和Nvidia Tesla K40 顯卡的工作站.由于與球體顆粒相比,超二次曲面單元間的接觸判斷更加復(fù)雜,其采用非線性牛頓迭代算法經(jīng)多次迭代才可計算出單元間的重疊量.因此,對于非球形顆粒系統(tǒng),超二次曲面單元比球形單元具有更低的計算效率和更長的計算時間.由于本文研究所模擬的實驗過程歷時較長,且Asencio 等[23]的實驗系統(tǒng)中顆粒數(shù)目(25 000 個)對于模擬而言計算耗費較大,因此選取立方體顆?？倲?shù)為800 以觀察其有序堆積現(xiàn)象.

2 結(jié)果與討論

2.1 立方體顆粒有序堆積的對比驗證

Asencio 等[23]的實驗證實當(dāng)圓筒扭轉(zhuǎn)加速度大于一定值時,圓筒內(nèi)的立方體顆粒會形成密實有序的堆積體.本文的模擬中也觀察到了類似的結(jié)果,如圖1 所示.這里引入兩個統(tǒng)計量表征立方體的有序堆積過程.體積分數(shù) φ (填充率)可以描述顆粒體系的密實狀態(tài).有序度參量 S4通過立方體顆粒的朝向計算得到,用于定量描述立方體顆粒體系的有序化堆積程度.

體積分數(shù) φ 的一般定義為總顆粒的實際體積除以顆粒堆積體占據(jù)的空間體積,即

其中,heq是顆?？傮w積除以容器底面積得到的等效圓柱體高度,是立方體顆粒堆積體上表面的平均高度.上表面的平均高度的計算方法是將顆粒域劃分為一系列適當(dāng)大小的網(wǎng)格,標(biāo)記每列網(wǎng)格中(沿容器高度方向的若干網(wǎng)格稱為一列)位置最高的顆粒的z 坐標(biāo),再對所有列位置最高顆粒的z 坐標(biāo)求平均值并加上顆粒的半邊長修正得到上表面的平均高度式(23)中,N 為顆?？倲?shù),單一顆粒體積按照標(biāo)準立方體體積計算.其圓角缺失的部分體積根據(jù)高斯-勒讓德積分公式計算為標(biāo)準立方體體積的6%,因此式(23)計算得到的體積分數(shù)會略大于真實的體積分數(shù).本文選用式(23)計算體積分數(shù)是為了與文獻[23]的體積分數(shù)計算方法保持一致,以便于數(shù)值模擬結(jié)果和實驗結(jié)果相對比.

有序度參量 S4[33]為

其中 μi為立方體顆粒i 的表面外法線方向,nb=(0,0,1)T為轉(zhuǎn)筒的轉(zhuǎn)軸方向.一般的超二次曲面顆粒有三個顆粒主軸,而對于具有良好對稱性的立方體顆粒,三個顆粒主軸無需進行區(qū)分.因此,在數(shù)值計算模型中,改進的 S4計算公式為

其方中向 向j=量1,.2m,3a;xe(ieji表j·n示b)顆表粒示i的3 個3顆個粒顆主粒軸主的軸單的單位位方向向量分別與 nb點乘的最大值.數(shù)值計算模型中,顆粒各主軸方向由四元數(shù)確定,因此,顆粒主軸在全局坐標(biāo)系下的方向向量為

圖2 為 γ=1.01,ω0=5 rad/s 時,φ 和 S4隨旋轉(zhuǎn)次數(shù)的變化.其中圖2(a)為文獻[23]的實驗結(jié)果.如圖所示,數(shù)值模擬得到的 φ 隨旋轉(zhuǎn)次數(shù)增加的變化趨勢同實驗結(jié)果所展示的變化趨勢基本一致,呈現(xiàn)隨著旋轉(zhuǎn)次數(shù)對數(shù)值的增加逐漸增大直至穩(wěn)定的特點.這種變化趨勢與Nowak 等[34]通過敲擊實現(xiàn)球形顆粒隨機堆積體系密實化所發(fā)現(xiàn)的 φ 的變化規(guī)律完全一致.這意味著對于立方體顆粒,容器的旋轉(zhuǎn)剪切扮演了敲擊的作用,通過這種簡單的人為設(shè)計擾動方案可以得到接近立方體顆粒最密規(guī)則堆積的結(jié)構(gòu).本文所得到的體積分數(shù)達到最終穩(wěn)定狀態(tài)時的數(shù)值同實驗結(jié)果略有差異.本文不同邊界運動條件下的數(shù)值模擬最終結(jié)果均為一個7 層的立方體顆粒堆積體,最上層內(nèi)的顆粒約為69～72 個,設(shè)置這些顆粒是為了保證下部6 層的顆粒數(shù)量達到飽和,其他各層顆粒約為119～123 個.若去掉最上層顆粒,φ 的值約為0.88.由于正方體邊長去貼合容器邊壁時難免有空隙且模擬計算的顆粒數(shù)小于實驗體系,這也使模擬所得到的 φ 略低于實驗值0.96.此外模擬得到的 φ 達到穩(wěn)定狀態(tài)后漲落相較實驗結(jié)果更為明顯.這是因為當(dāng)?shù)诹鶎宇w粒填滿后,下面六層的排布結(jié)構(gòu)就不再顯著地變化,但是最上層的顆粒還會伴隨下層顆粒的整體轉(zhuǎn)動產(chǎn)生移動和翻轉(zhuǎn),對于體積分數(shù)的計算結(jié)果帶來一定的漲落.如果去掉第七層顆粒,體積分數(shù)的漲落是非常小的.

圖2 體積分數(shù)及有序度隨旋轉(zhuǎn)次數(shù)的變化Fig.2 Packing fraction Φ and cubatic order parameter S4 versus number of rotation

有序度參量 S4隨旋轉(zhuǎn)次數(shù)增加而呈現(xiàn)的變化趨勢與 φ 的變化趨勢相同.S4最終穩(wěn)定在0.95 左右,標(biāo)志著體系完成了從無序向有序的轉(zhuǎn)變.圖中 S4的變化存在一個比較明顯的轉(zhuǎn)折點,可以定義一個特征旋轉(zhuǎn)次數(shù) Nr,用于標(biāo)記體系何時完成有序堆積.相較于實驗結(jié)果,模擬計算得到的 Nr小一到兩個量級,可以認為這是由于離散元模擬所用的顆粒數(shù)遠小于實驗顆粒數(shù)目所導(dǎo)致的.數(shù)值模擬選取的容器底面直徑和顆粒邊長比為13.4,而實驗中容器底面直徑和顆粒邊長比約為35.下一小節(jié)中,可以清楚地看到,邊壁的旋轉(zhuǎn)剪切作用所導(dǎo)致的顆粒體系的密實有序化過程是逐漸的從外側(cè)向內(nèi)側(cè)演化發(fā)展的.模擬僅需要完成4～5 層的顆粒同心圓排列就基本完成有序化過程,而實驗則需要完成15 層左右的同心圓排列.考慮到往復(fù)旋轉(zhuǎn)邊界的擾動影響必然是通過最靠近筒壁的顆粒逐漸向內(nèi)部傳遞,因此模擬結(jié)果得到的Nr小于實驗結(jié)果是可以理解的.盡管如此,模擬得到的 φ 和 S4的變化趨勢,可以很好地反映立方體顆粒在邊界往復(fù)旋轉(zhuǎn)作用下的密實化效果和有序化過程,這表明基于超二次曲面顆粒單元的DEM 數(shù)值模型的有效性.

2.2 立方體顆粒有序堆積過程的結(jié)構(gòu)演化

圖3 給出了 γ=4.04,ω0=10 rad/s 時,立方體顆粒堆積體在側(cè)面、半剖面以及底面視角下的有序化過程.如圖所示,初始時刻體系處于隨機堆積的狀態(tài)有序度較低,底面上顆粒分布相對松散.隨著旋轉(zhuǎn)剪切次數(shù)的增加,幾十轉(zhuǎn)內(nèi)位于頂層的顆粒(藍色)會不斷的通過縫隙掉落至下層顆粒(紅色)中,體系將迅速完成主要的密實化過程.同時,靠近底面和容器邊壁的顆粒呈現(xiàn)某一面平行容器底面的有序狀態(tài),而靠近轉(zhuǎn)軸的內(nèi)部顆粒則依然處于無序的狀態(tài).之后的幾百轉(zhuǎn)內(nèi),最低的6 層顆粒逐漸填滿并全部實現(xiàn)某一面平行容器底面的有序排列,最上一層顆粒則依然會出現(xiàn)移動和翻轉(zhuǎn).從底面層顆粒的速度分布(彩圖)可以看出,最靠近容器邊壁的一圈顆粒速度最大,越靠近轉(zhuǎn)軸的顆粒速度越小,顆粒間的速度差會造成顆粒間的切向滑動摩擦力.這證明邊壁擾動確實是以摩擦力的形式促使顆粒的有序堆積從靠近筒壁的顆粒逐漸向內(nèi)演化.

圖3 立方體顆粒的有序化過程Fig.3 The ordering process of cubic particles

進一步地細致觀察發(fā)現(xiàn)初始時刻立方體顆粒完成沉降堆積后,雖然整體的有序度比較低,但是單一顆粒和其相鄰顆粒間并不會存在巨大的取向差異.Boton 等[35]研究結(jié)果表明,具有平面的非球形顆粒相互接觸時,面與面的接觸方式最大程度阻礙顆粒的轉(zhuǎn)動且最容易滿足顆粒堆積體的穩(wěn)定性要求,因此具有平面的規(guī)則非球形顆粒很容易形成特定方向的排列,對于具有高度對稱性的立方體顆粒更是如此.Bautista-Carbajal[36]給出了二維狹縫中方形顆粒的團簇相圖,結(jié)果表明狹縫中的方形顆粒有三種穩(wěn)定的狀態(tài),分別是邊長平行狹縫邊界,對角線垂直狹縫邊界以及兩者的混合狀態(tài).對于三維立方體顆粒,Zuriguel 和Maza 等[22,37-38]系列工作表明方形顆粒在容器中的沉降極易形成面對角線沿重力方向(邊與水平面夾角45°)的團簇.這一特性可以通過平均場近似解釋.力鏈穿過給定粒子的概率應(yīng)與其側(cè)面在水平方向上的相對投影成正比,方形顆粒在這一角度下力鏈的穿透概率取極大值,此時顆粒體系內(nèi)部的沿重力方向的力鏈分布最接近高斯分布狀態(tài),顆粒間具有最大的接觸面積.從側(cè)視圖可以看出,在圓筒邊壁最初的幾十轉(zhuǎn)內(nèi),重力誘導(dǎo)顆粒迅速的掉落,容易形成由幾個或者幾十個顆粒組成的面對角線沿重力方向的團簇.后續(xù)邊壁的往復(fù)旋轉(zhuǎn)則是誘導(dǎo)這些團簇整體轉(zhuǎn)動形成立方體顆粒某一面平行容器底面的團簇.文獻[23]認為重力因素和邊壁旋轉(zhuǎn)因素對顆粒團簇的最終形態(tài)構(gòu)成競爭關(guān)系.當(dāng)旋轉(zhuǎn)擾動強度足夠劇烈時,團簇會以整體旋轉(zhuǎn)的形式完成排列取向的轉(zhuǎn)變.當(dāng)旋轉(zhuǎn)擾動強度較低時,體系最終為面對角線沿重力方向的顆粒團簇和某一面平行容器底面的顆粒團簇共存的狀態(tài),例如圖4 內(nèi)插圖所示結(jié)構(gòu),此時 S4穩(wěn)定在0.82 左右.

圖4 γ 一定,不同峰值角速度 ω0 下,填充率 Φ 和有序度參量 S 4 的變化Fig.4 Volume fraction Φ and cubatic order parameter S4 versus number of rotation N for γ=1.01 and differentω0

2.3 立方體顆粒有序堆積的影響因素

Asencio 等[23]認為每次旋轉(zhuǎn)時,對顆粒體系施加的 γ 對立方體顆粒的有序堆積起決定性的作用.如果 γ 大于0.5,那么立方體骰子必然能在容器約10 萬次的往復(fù)旋轉(zhuǎn)后達到有序堆積狀態(tài).此時每一層的立方體顆粒都排列成幾乎完美的同心環(huán).相反,如果 γ 低于0.5,那么有序堆積的速度就會變得非常緩慢,以至于很難判斷顆粒是否最終能達到最密堆積狀態(tài).即使經(jīng)過10 萬次的旋轉(zhuǎn),靠近圓柱體中心的顆粒仍然處于高度混亂的狀態(tài).本文的模擬方案中,容器 γ 由 ω0和反轉(zhuǎn)時間 Δt 決定.不同的 ω0和Δt可能對應(yīng)相同的γ .因此設(shè)計了一組 γ 相同但 ω0不同的容器運動條件,考察立方體顆粒的有序堆積過程,如圖4 所示.如果 γ 是立方體顆粒有序堆積的控制變量,那么可以預(yù)期只要 γ相同,盡管 ω0不同,容器各運動條件下顆粒有序堆積過程應(yīng)該是幾乎相同的.但是模擬結(jié)果顯示 γ 一定時,隨著 ω0的增大(相應(yīng)的反轉(zhuǎn)時間 Δt 增加),顆粒體系只需更少的往復(fù)旋轉(zhuǎn)次數(shù)就可達到最密實狀態(tài)并完成顆粒的有序排列.應(yīng)當(dāng)指出,ω0=2.50 rad/s的這組模擬,當(dāng)往復(fù)旋轉(zhuǎn)次數(shù)達到4000 次時仍未完成完全的有序化轉(zhuǎn)變(有序度參數(shù) S4約為0.82),這是因為靠近容器邊壁的部分顆粒形成面對角線沿重力方向的團簇,如圖4 內(nèi)插圖所示.無法確定更多的往復(fù)旋轉(zhuǎn)次數(shù)是否一定能實現(xiàn)所有顆粒均為某一面平行水平面的有序排列.這種情況與文獻[23]報道的低于0.5 時的旋轉(zhuǎn)剪切結(jié)果類似,容器內(nèi)兩種取向的團簇顆粒共存.這意味著立方體顆粒通過邊界的剪切作用完成有序堆積確實需要外界給予的擾動強度達到一定的值,但將作為衡量擾動強弱的指標(biāo)并不完備,模擬結(jié)果顯示這個強度指標(biāo)與 ω0和Δt 有關(guān).

為了進一步研究 ω0和 Δt 對立方體顆粒有序堆積過程的影響,首先設(shè)計了一組 ω0相同但 γ 不同的容器運動條件,模擬結(jié)果如圖5 所示.ω0一定時,隨著的 γ 增大,立方體顆粒需要更多的往復(fù)旋轉(zhuǎn)次數(shù)才能達到最密實狀態(tài)并完成顆粒的有序化轉(zhuǎn)變.γ 對顆粒的密實化和有序化呈現(xiàn)消極影響,這一結(jié)果與文獻[23]矛盾.同時設(shè)計了一組 Δt 相同但 ω0不同的運動方案,此時隨ω0的增加 γ等比例增大,如圖6 所示.結(jié)果顯示更高的 γ 下立方體顆粒確實僅需要更少的往復(fù)旋轉(zhuǎn)次數(shù)就可以完成顆粒的密實化和有序排列,這一結(jié)果與文獻[23]一致.文獻[23]并未給出轉(zhuǎn)筒詳細的工作條件,但提及到實驗采用的轉(zhuǎn)筒僅在反轉(zhuǎn)時對顆粒有剪切作用.可以認為實驗中容器的運動模式接近為方波形式,角速度反向所需要的時間可能為定值,如此 γ 的增大完全由峰值角速度的增大所導(dǎo)致,類似圖6 結(jié)果所對應(yīng)的運動方案.

圖5 角速度 ω0 一定,不同 γ 下,填充率 Φ 和有序度參量 S 4 的變化Fig.5 Volume fraction Φ and cubatic order parameter S4 versus number of rotation N for ω0=5.0 rad/s and differentγ

圖6 反轉(zhuǎn)時間一定,不同 γ 下,填充率 Φ 和有序度參量 S 4 的變化Fig.6 Volume fraction Φ and cubatic order parameter S4 versus number of rotation N for the constant inversion time and differentγ

綜合以上結(jié)果,本文提出用容器的旋轉(zhuǎn)角位移β作為衡量擾動強弱的指標(biāo)更為合理.圖7 給出了旋轉(zhuǎn)角位移 β 和特征旋轉(zhuǎn)次數(shù) Nr的關(guān)系,圖中所有模擬組的 S4最終的穩(wěn)定值均大于0.94,可認為體系完成了有序堆積.對 β 與 Nr繪制雙對數(shù)散點圖并進行線性回歸,得到 β 與 Nr之間的關(guān)系為

圖7 特征旋轉(zhuǎn)次數(shù)與旋轉(zhuǎn)角位移的關(guān)系Fig.7 Characteristics rotation number versus angular displacement

考慮到實驗誤差,可以認為 Nr和 β 之間為簡單的反比例關(guān)系即

其中 Nr0為系數(shù)取決于顆粒體系的數(shù)量,容器能容納的單層顆粒越多則 Nr0越大.利用這種反比例關(guān)系,很容易理解模擬中發(fā)現(xiàn)的S4曲線重合的現(xiàn)象.比如邊界 γ=1.01,ω0=5 rad/s 和 γ=4.04,ω0=10 rad/s 的兩個邊界運動方案,對應(yīng)的旋轉(zhuǎn)角位移均為0.338 rad,因此兩者的 Nr幾乎完全相同,如圖8 所示.

圖8 相同 β 值下,填充率 Φ 和有序度參量 S 4 的變化曲線.Fig.8 Volume fraction Φ and cubatic order parameter S4 versus number of rotation N for the constantβ

這種反比例關(guān)系是因為,容器旋轉(zhuǎn)擾動作為類似熱漲落的驅(qū)動力,提供了立方體顆粒有序堆積的能量.考慮圓形底面的對稱性,可以近似將所有顆粒受到的來自邊界的摩擦力簡化為對全體顆粒的力偶,單次旋轉(zhuǎn)時邊界對顆粒體系做功等于力偶矩乘以旋轉(zhuǎn)角位移.旋轉(zhuǎn)角位移越大,容器邊界輸入能量越多,顆粒越容易產(chǎn)生運動,直至完成有序排列,此時體系的勢能最低,顆粒單元與相鄰顆粒間的面-面接觸約束最強烈.當(dāng)然式(28)的適用范圍一定存在上下限.對于適用上限可以合理想象即使旋轉(zhuǎn)角位移非常大,但是完成有序堆積依然需要一定次數(shù)的往復(fù)旋轉(zhuǎn)而不僅僅是一兩次的旋轉(zhuǎn).對于適用下限上文已經(jīng)闡述當(dāng)旋轉(zhuǎn)角位移小于某一臨界值時,體系最終會形成兩種形式的團簇顆粒共存的結(jié)構(gòu),只是下限的臨界值所對應(yīng)的力學(xué)機理有待進一步研究.

2.4 不同重力加速度下的立方體顆粒剪切排序

深空環(huán)境中,固體星球表面聚集著大量的顆粒物質(zhì).人類對于深空環(huán)境下顆粒物質(zhì)的物理規(guī)律和操控方式的認知仍比較淺顯.比如太陽系中大行星表面的風(fēng)化層典型顆粒粒徑小于小行星表面風(fēng)化層的典型顆粒粒徑,即使是這樣簡單的現(xiàn)象,人們對于其物理機制的認知尚未完全統(tǒng)一.深空環(huán)境的顆粒堆積問題亦是如此.為研究重力加速度對立方體顆粒受轉(zhuǎn)動驅(qū)動實現(xiàn)有序堆積的影響,設(shè)計了三種不同重力加速度下(地球、火星和月球的重力加速度)的邊界運動條件.鑒于顆粒物質(zhì)具有非線性和對微小作用力敏感的特性,初始堆積狀態(tài)將影響后續(xù)堆積演變特征,所以本文三種重力加速度下的顆粒堆積結(jié)構(gòu)保持不變.且同時保持 γ=4.04,ω0=10 rad/s 不變,模擬結(jié)果如圖9 所示.地球和火星的重力加速度下,立方體顆粒都實現(xiàn)了無序向有序的轉(zhuǎn)變并達到了最密堆積和完全有序狀態(tài),但火星重力加速度下達到有序狀態(tài)需要的往復(fù)旋轉(zhuǎn)次數(shù)更多.而月球重力加速度下顆粒體系雖然 φ 和 S4也隨著邊界不斷旋轉(zhuǎn)而變大,但無法達到完全有序和最密實狀態(tài).上述結(jié)果表明,亞重力加速度條件下,容器的往復(fù)旋轉(zhuǎn)依然是實現(xiàn)立方體顆粒有序堆積的有效方式,但重力加速度的減少會對立方體顆粒的有序堆積過程產(chǎn)生抑制效果.這是因為隨著重力加速度的減小,顆粒間的摩擦力相應(yīng)減小,不利于顆粒受到擾動后尋找下沉路徑.本文也曾嘗試過模擬十分之一地球重力加速度以及更低重力加速度條件下的立方體顆粒有序堆積過程,但上層顆粒極容易形成類似顆粒氣體的漂浮狀態(tài),無法再將其視為一個類似固體的堆積體.

圖9 不同重力加速度下,填充率 Φ 和有序度參量 S4 的變化曲線(模擬參數(shù):邊界 γ=4.04,峰值角速度 ω0=10.0 rad/s)Fig.9 Volume fraction Φ and order parameter S4 versus number of rotation N for different acceleration of gravity (γ=4.04,ω0=10.0 rad/s)

3 結(jié)論

基于超二次曲面方程構(gòu)建立方體顆粒形態(tài),采用離散元方法針對圓筒內(nèi)隨機堆積立方體顆粒在容器邊界往復(fù)旋轉(zhuǎn)作用下的有序化堆積行為開展數(shù)值模擬研究,復(fù)現(xiàn)了實驗中觀察到的立方體顆粒實現(xiàn)最密規(guī)則堆積的過程.通過本文研究,得到如下結(jié)論.

(1)立方體顆粒系統(tǒng)的體積分數(shù)和有序度隨著容器往復(fù)旋轉(zhuǎn)次數(shù)增加而逐漸增大直至穩(wěn)定.顆粒系統(tǒng)的有序堆積過程表現(xiàn)為由外層顆粒向內(nèi)部顆粒逐漸擴展的演化趨勢.

(2)容器旋轉(zhuǎn)角位移是立方體顆粒實現(xiàn)有序堆積過程的控制變量.顆粒完成有序堆積過程所需的特征旋轉(zhuǎn)次數(shù)與容器旋轉(zhuǎn)角位移呈現(xiàn)反比例關(guān)系.旋轉(zhuǎn)角位移過低,體系只會形成某一面平行水平面的顆粒團簇與面對角線平行重力的顆粒團簇共存的結(jié)構(gòu).

(3)在亞重力環(huán)境下立方體顆粒同樣可以通過容器的往復(fù)旋轉(zhuǎn)實現(xiàn)最密規(guī)則堆積.重力加速度的減少會抑制立方體顆粒從無序堆積向有序堆積的轉(zhuǎn)變.

致謝

本項目承蒙中國載人航天工程國際合作項目支持.