張 祺 喬 婷 季順迎 厚美瑛

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

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

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

引言

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

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

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

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

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

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

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

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

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

顆粒接觸檢測的主要內(nèi)容是檢測顆粒是否接觸以及確定接觸點(diǎn)和接觸法向.超二次曲面單元間的接觸檢測問題可以考慮為一個(gè)定義在全局坐標(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μ .計(jì)算結(jié)果 X0滿足Fi(X0)＜0 且 Fj(X0)＜0 則表明兩個(gè)單元間發(fā)生接觸.接觸法向 n 可表示為n=?Fi(X0)/||?Fi(X0)|| .進(jìn)一步考慮Xi=X0+αn 和 Xj=X0+βn,這里未知參數(shù) α 和β 可通過一元非線性牛頓迭代得到:因此,法向重疊量可表述為 δn=Xi-Xj.

1.3 接觸模型及接觸力的計(jì)算

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

法向作用力可以表示為

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

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

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

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

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

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

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

1.4 圓角立方體顆粒的運(yùn)動(dòng)控制方程

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

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

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

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

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

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

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

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

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

表1 轉(zhuǎn)動(dòng)驅(qū)動(dòng)立方體顆粒有序堆積過程的主要計(jì)算參數(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)過程需要的時(shí)間.使用重力加速度無量綱化后的扭轉(zhuǎn)加速度 γ 為

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

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

2 結(jié)果與討論

2.1 立方體顆粒有序堆積的對比驗(yàn)證

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

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

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

有序度參量 S4[33]為

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

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

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

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

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

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

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

進(jìn)一步地細(xì)致觀察發(fā)現(xiàn)初始時(shí)刻立方體顆粒完成沉降堆積后,雖然整體的有序度比較低,但是單一顆粒和其相鄰顆粒間并不會(huì)存在巨大的取向差異.Boton 等[35]研究結(jié)果表明,具有平面的非球形顆粒相互接觸時(shí),面與面的接觸方式最大程度阻礙顆粒的轉(zhuǎn)動(dòng)且最容易滿足顆粒堆積體的穩(wěn)定性要求,因此具有平面的規(guī)則非球形顆粒很容易形成特定方向的排列,對于具有高度對稱性的立方體顆粒更是如此.Bautista-Carbajal[36]給出了二維狹縫中方形顆粒的團(tuán)簇相圖,結(jié)果表明狹縫中的方形顆粒有三種穩(wěn)定的狀態(tài),分別是邊長平行狹縫邊界,對角線垂直狹縫邊界以及兩者的混合狀態(tài).對于三維立方體顆粒,Zuriguel 和Maza 等[22,37-38]系列工作表明方形顆粒在容器中的沉降極易形成面對角線沿重力方向(邊與水平面夾角45°)的團(tuán)簇.這一特性可以通過平均場近似解釋.力鏈穿過給定粒子的概率應(yīng)與其側(cè)面在水平方向上的相對投影成正比,方形顆粒在這一角度下力鏈的穿透概率取極大值,此時(shí)顆粒體系內(nèi)部的沿重力方向的力鏈分布最接近高斯分布狀態(tài),顆粒間具有最大的接觸面積.從側(cè)視圖可以看出,在圓筒邊壁最初的幾十轉(zhuǎn)內(nèi),重力誘導(dǎo)顆粒迅速的掉落,容易形成由幾個(gè)或者幾十個(gè)顆粒組成的面對角線沿重力方向的團(tuán)簇.后續(xù)邊壁的往復(fù)旋轉(zhuǎn)則是誘導(dǎo)這些團(tuán)簇整體轉(zhuǎn)動(dòng)形成立方體顆粒某一面平行容器底面的團(tuán)簇.文獻(xiàn)[23]認(rèn)為重力因素和邊壁旋轉(zhuǎn)因素對顆粒團(tuán)簇的最終形態(tài)構(gòu)成競爭關(guān)系.當(dāng)旋轉(zhuǎn)擾動(dòng)強(qiáng)度足夠劇烈時(shí),團(tuán)簇會(huì)以整體旋轉(zhuǎn)的形式完成排列取向的轉(zhuǎn)變.當(dāng)旋轉(zhuǎn)擾動(dòng)強(qiáng)度較低時(shí),體系最終為面對角線沿重力方向的顆粒團(tuán)簇和某一面平行容器底面的顆粒團(tuán)簇共存的狀態(tài),例如圖4 內(nèi)插圖所示結(jié)構(gòu),此時(shí) 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]認(rèn)為每次旋轉(zhuǎn)時(shí),對顆粒體系施加的 γ 對立方體顆粒的有序堆積起決定性的作用.如果 γ 大于0.5,那么立方體骰子必然能在容器約10 萬次的往復(fù)旋轉(zhuǎn)后達(dá)到有序堆積狀態(tài).此時(shí)每一層的立方體顆粒都排列成幾乎完美的同心環(huán).相反,如果 γ 低于0.5,那么有序堆積的速度就會(huì)變得非常緩慢,以至于很難判斷顆粒是否最終能達(dá)到最密堆積狀態(tài).即使經(jīng)過10 萬次的旋轉(zhuǎn),靠近圓柱體中心的顆粒仍然處于高度混亂的狀態(tài).本文的模擬方案中,容器 γ 由 ω0和反轉(zhuǎn)時(shí)間 Δt 決定.不同的 ω0和Δt可能對應(yīng)相同的γ .因此設(shè)計(jì)了一組 γ 相同但 ω0不同的容器運(yùn)動(dòng)條件,考察立方體顆粒的有序堆積過程,如圖4 所示.如果 γ 是立方體顆粒有序堆積的控制變量,那么可以預(yù)期只要 γ相同,盡管 ω0不同,容器各運(yùn)動(dòng)條件下顆粒有序堆積過程應(yīng)該是幾乎相同的.但是模擬結(jié)果顯示 γ 一定時(shí),隨著 ω0的增大(相應(yīng)的反轉(zhuǎn)時(shí)間 Δt 增加),顆粒體系只需更少的往復(fù)旋轉(zhuǎn)次數(shù)就可達(dá)到最密實(shí)狀態(tài)并完成顆粒的有序排列.應(yīng)當(dāng)指出,ω0=2.50 rad/s的這組模擬,當(dāng)往復(fù)旋轉(zhuǎn)次數(shù)達(dá)到4000 次時(shí)仍未完成完全的有序化轉(zhuǎn)變(有序度參數(shù) S4約為0.82),這是因?yàn)榭拷萜鬟叡诘牟糠诸w粒形成面對角線沿重力方向的團(tuán)簇,如圖4 內(nèi)插圖所示.無法確定更多的往復(fù)旋轉(zhuǎn)次數(shù)是否一定能實(shí)現(xiàn)所有顆粒均為某一面平行水平面的有序排列.這種情況與文獻(xiàn)[23]報(bào)道的低于0.5 時(shí)的旋轉(zhuǎn)剪切結(jié)果類似,容器內(nèi)兩種取向的團(tuán)簇顆粒共存.這意味著立方體顆粒通過邊界的剪切作用完成有序堆積確實(shí)需要外界給予的擾動(dòng)強(qiáng)度達(dá)到一定的值,但將作為衡量擾動(dòng)強(qiáng)弱的指標(biāo)并不完備,模擬結(jié)果顯示這個(gè)強(qiáng)度指標(biāo)與 ω0和Δt 有關(guān).

為了進(jìn)一步研究 ω0和 Δt 對立方體顆粒有序堆積過程的影響,首先設(shè)計(jì)了一組 ω0相同但 γ 不同的容器運(yùn)動(dòng)條件,模擬結(jié)果如圖5 所示.ω0一定時(shí),隨著的 γ 增大,立方體顆粒需要更多的往復(fù)旋轉(zhuǎn)次數(shù)才能達(dá)到最密實(shí)狀態(tài)并完成顆粒的有序化轉(zhuǎn)變.γ 對顆粒的密實(shí)化和有序化呈現(xiàn)消極影響,這一結(jié)果與文獻(xiàn)[23]矛盾.同時(shí)設(shè)計(jì)了一組 Δt 相同但 ω0不同的運(yùn)動(dòng)方案,此時(shí)隨ω0的增加 γ等比例增大,如圖6 所示.結(jié)果顯示更高的 γ 下立方體顆粒確實(shí)僅需要更少的往復(fù)旋轉(zhuǎn)次數(shù)就可以完成顆粒的密實(shí)化和有序排列,這一結(jié)果與文獻(xiàn)[23]一致.文獻(xiàn)[23]并未給出轉(zhuǎn)筒詳細(xì)的工作條件,但提及到實(shí)驗(yàn)采用的轉(zhuǎn)筒僅在反轉(zhuǎn)時(shí)對顆粒有剪切作用.可以認(rèn)為實(shí)驗(yàn)中容器的運(yùn)動(dòng)模式接近為方波形式,角速度反向所需要的時(shí)間可能為定值,如此 γ 的增大完全由峰值角速度的增大所導(dǎo)致,類似圖6 結(jié)果所對應(yīng)的運(yùn)動(dò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)時(shí)間一定,不同 γ 下,填充率 Φ 和有序度參量 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)角位移β作為衡量擾動(dòng)強(qiáng)弱的指標(biāo)更為合理.圖7 給出了旋轉(zhuǎn)角位移 β 和特征旋轉(zhuǎn)次數(shù) Nr的關(guān)系,圖中所有模擬組的 S4最終的穩(wěn)定值均大于0.94,可認(rèn)為體系完成了有序堆積.對 β 與 Nr繪制雙對數(shù)散點(diǎn)圖并進(jìn)行線性回歸,得到 β 與 Nr之間的關(guān)系為

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

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

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

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

圖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)隨機(jī)堆積立方體顆粒在容器邊界往復(fù)旋轉(zhuǎn)作用下的有序化堆積行為開展數(shù)值模擬研究,復(fù)現(xiàn)了實(shí)驗(yàn)中觀察到的立方體顆粒實(shí)現(xiàn)最密規(guī)則堆積的過程.通過本文研究,得到如下結(jié)論.

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

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

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

致謝

本項(xiàng)目承蒙中國載人航天工程國際合作項(xiàng)目支持.