鄭麟 莫松平 李玉秀? 陳穎 徐進(jìn)良

1)(廣東工業(yè)大學(xué)材料與能源學(xué)院,廣州 510006)

2)(華北電力大學(xué)能源動力與機(jī)械工程學(xué)院,北京 102206)

1 引 言

顆粒是自然界常見的離散物質(zhì),靜止時(shí)呈現(xiàn)固體狀態(tài),而整體運(yùn)動時(shí)呈現(xiàn)類流體狀態(tài)[1].不同的顆粒在剪切或者振動等外力激勵(lì)條件下,會出現(xiàn)不同的移動分離和混合現(xiàn)象[2](花瓣形,圣誕樹形,指形,三明治[3]等).顆粒的分離與混合廣泛應(yīng)用在制藥、食品、化工和能源等領(lǐng)域[4],利用二元顆粒密度、尺寸、幾何形貌以及其他性質(zhì)的差異,施加外部作用,能夠?qū)崿F(xiàn)顆粒分離或混合[5-7].除工業(yè)領(lǐng)域,自然界的山體滑坡、泥石流、雪崩中亦有顆粒獨(dú)特的分離運(yùn)動,如大石塊聚集在頂部,小石塊沉降在底部的現(xiàn)象[6].正因?yàn)轭w粒運(yùn)動的高應(yīng)用價(jià)值和復(fù)雜性,在過去的幾十年中,顆粒流一直是研究熱點(diǎn)[8].

針對剪切二元顆粒體系分離的現(xiàn)象已經(jīng)開展了實(shí)驗(yàn)[9-12]、數(shù)值模擬[13]研究,而對分離現(xiàn)象的理論研究以三種典型機(jī)制為主導(dǎo): “滲流”機(jī)制理論[14];“對流”機(jī)制[15]; “凝聚”機(jī)制[16]; 顆粒分離的新機(jī)制[17-19],對于不同的分離運(yùn)動傾向于用多種機(jī)制競爭或主導(dǎo)機(jī)制解釋.隨著新分離現(xiàn)象的發(fā)現(xiàn)不斷提出不同機(jī)制,因此目前依舊沒有統(tǒng)一的顆粒運(yùn)動的物理解釋.主要原因在于以往分析集中于對顆粒群使用經(jīng)典流體理論類比分析,顆粒系統(tǒng)作為軟凝聚物質(zhì)系統(tǒng)不能被其套用; 其次,缺乏物理機(jī)制的揭示[11,20],包括1)缺乏底層顆粒在向上躍升過程中的運(yùn)動學(xué)以及動力學(xué)行為特征分析,2)缺乏討論底層顆粒向上起跳的控制條件,3)缺乏研究細(xì)觀顆粒群的空間分布結(jié)構(gòu)對顆粒分離過程的影響規(guī)律.本文的重要工作就是探討以上三點(diǎn).同時(shí),顆粒間的摩擦耗散作用是區(qū)別于其他物質(zhì)狀態(tài)的主要特征[21],摩擦力的變化導(dǎo)致顆粒系統(tǒng)的能量耗散變化,甚至導(dǎo)致顆粒分離的變化發(fā)生逆轉(zhuǎn)[22],因此研究摩擦系數(shù)對顆粒運(yùn)動的變化具有重要意義.

本文采用三維離散元(DEM)方法,模擬大小兩種尺度球形顆粒初始時(shí)刻處于底層的大顆粒向頂層躍升的運(yùn)動過程,研究顆粒間的摩擦作用強(qiáng)弱對顆粒上升變化的影響,獲得摩擦力對顆粒躍升行為動力學(xué)的影響規(guī)律,揭示大顆粒躍升的物理機(jī)制,為后續(xù)人們理解顆粒分離運(yùn)動理論提供參考依據(jù).

2 模擬設(shè)置

2.1 物理模型

本文模擬的二元顆粒系統(tǒng),周側(cè)固定的兩個(gè)內(nèi)外圓柱面和底部旋轉(zhuǎn)的圓環(huán)形底盤構(gòu)成,7501顆顆粒(7500顆直徑為0.002 m和1顆直徑為0.003 m的球形顆粒)自由堆積在圓環(huán)槽型限域空間中,如圖1所示.底面以一定的角速度旋轉(zhuǎn),對槽內(nèi)的顆粒群產(chǎn)生恒定剪切力.圓環(huán)槽內(nèi)徑為0.294 m,外徑為0.317 m,高度為顆粒群的堆積高度.模擬系統(tǒng)置于重力環(huán)境,以剪切盒的底部中心為坐標(biāo)原點(diǎn),重力反方向?yàn)閦軸正方向,坐標(biāo)軸設(shè)立滿足右手系定則.顆粒系統(tǒng)的填注方式如下: 首先將大顆粒置于剪切槽底層,然后在容器上方隨機(jī)生成小顆粒傾倒入槽,釋放完預(yù)設(shè)數(shù)目的小顆粒后再合上頂盤.系統(tǒng)在重力條件下弛豫足夠長時(shí)間以確?？倓幽芊€(wěn)定在較低的范圍內(nèi)后旋轉(zhuǎn)底盤,旋轉(zhuǎn)周期為20 s.觀察大顆粒在小顆粒群中的運(yùn)動特征.

圖1 三維剪切顆粒體系示意圖Fig.1.Schematic diagram of shear granular system.

2.2 數(shù)學(xué)模型

模擬中顆粒的運(yùn)動包括平動和轉(zhuǎn)動,分別采用(1)式和(2)式進(jìn)行描述

其中Fn是顆粒間碰撞的法向接觸力,Ft是顆粒間碰撞的切向接觸力,Fb是除碰撞力以外的力.重力及顆粒間的法向和切向接觸力使顆粒發(fā)生平動; 切向接觸力產(chǎn)生力矩,導(dǎo)致顆粒發(fā)生旋轉(zhuǎn).

顆粒間的相互碰撞采用彈簧阻尼軟球模型描述.彈簧模型考慮顆粒因彈性碰撞受到的作用力,顆粒在法向和切向分別發(fā)生彈性形變,彈性形變量與顆粒受到的彈性作用力成正比.阻尼模型模擬顆粒因非彈性形變導(dǎo)致的機(jī)械能損失,也可以分解為法向和切向分量.因此,當(dāng)球i,j發(fā)生碰撞時(shí),球i受到的切向和法向作用力為

其中kn表示法向彈性碰撞系數(shù),δnij是兩個(gè)顆粒之間的法向重疊距離;γn是法向黏彈性阻尼系數(shù),νnij是法向相對速度;kt是切向彈性碰撞系數(shù),δtij是兩個(gè)球狀顆粒的切向相對位移向量;γt是切向接觸的黏彈性阻尼系數(shù),νtij是切向相對速度.其中δtij必須滿足摩擦屈服準(zhǔn)則Ft≤XμFn,Xμ是最大靜摩擦系數(shù).以上參數(shù)的計(jì)算公式如下:

式中Ri,Rj為顆粒半徑,mi,mj是顆粒質(zhì)量,Yi,Yj是楊氏模量,νi,νj是泊松系數(shù),e為恢復(fù)系數(shù).

模擬采用美國圣地亞實(shí)驗(yàn)室開發(fā)的三維DEM程序LIGGGHTS[23].微分方程采用verlet積分算法.模擬中力學(xué)參數(shù)和其他參數(shù)的設(shè)置詳見表1.

表1 模擬參數(shù)Table 1. Simulation parameters.

本文的數(shù)值模型初始時(shí)刻大顆粒在下部,小顆粒在上部,經(jīng)過中間的摻混過程,最終達(dá)到位置顛倒的平衡狀態(tài),與May等[10]的實(shí)驗(yàn)結(jié)果一致,如圖2所示.驗(yàn)證了本文數(shù)值模型以及參數(shù)設(shè)置的可行性.

3 模擬結(jié)果與討論

模擬結(jié)果顯示,對于所有的摩擦系數(shù),大顆粒向上躍升的過程都會經(jīng)歷三個(gè)階段: 1)弛豫階段,大顆粒隨著剪切槽底面周向旋轉(zhuǎn),但高度位置不變;2)起跳階段,大顆粒經(jīng)歷極短的時(shí)間跨度迅速到達(dá)顆粒群頂層; 3)平衡階段,到達(dá)頂層的大顆粒在高度方向上下脈動,但不會再回落到底層,完成顆?？臻g位置的顛倒.以下的結(jié)果與討論,按照躍升過程經(jīng)歷的三個(gè)階段為順序展開,首先,以摩擦系數(shù)為0.3為例,介紹大顆粒運(yùn)動的整體特征; 其次,討論摩擦系數(shù)對大顆粒運(yùn)動的影響規(guī)律; 最后,探討顆粒起跳的物理機(jī)理.

3.1 大顆粒的整體運(yùn)動特性

圖3(a)為大顆粒三維空間位置隨時(shí)間的演變過程.由圖可見,大顆粒做逆時(shí)針旋轉(zhuǎn)上升運(yùn)動,三維軌跡可以分為三個(gè)階段: 1)定義t從0—29.92 s的時(shí)間跨度為起跳弛豫時(shí)間,大顆粒在剪切槽底部做旋轉(zhuǎn)運(yùn)動,在高度方向上變化很小,間歇性地發(fā)生小幅度不規(guī)則跳動; 2)t=29.92 s開始起跳上升,起跳的過程很短,經(jīng)歷3.88 s后,t=33.8 s時(shí)刻躍升過程完成; 3)大顆粒上升到頂部,在頂部平衡位置附近上下脈動.

為了進(jìn)一步探究顆粒軌跡特征,將圖3(a)所示的顆粒三維軌跡分解為水平方向的x和y位移曲線,如圖3(b)所示,以及高度方向上的z位移曲線,如圖3(c)所示.由圖3(b)可見,顆粒在水平方向的位移特征為類圓周運(yùn)動,周期為43 s,是剪切槽底面旋轉(zhuǎn)周期的一半,說明大顆粒的運(yùn)動落后于剪切槽底面的旋轉(zhuǎn)運(yùn)動,這是由顆粒與底面間的相對滑動引起的.圖3(c)所示的z方向位移時(shí)間序列清楚地描述了顆粒瞬間躍升的特征.大顆粒在整個(gè)過程中的運(yùn)動并非平穩(wěn)直線運(yùn)動,而是伴隨著回落的運(yùn)動,起跳上升運(yùn)動前大顆粒在底部發(fā)生多次假“起跳”,但均回落到底部,直到t=29.92 s后不再落回,而是逐漸上升到達(dá)頂部,最后停留在頂部做小幅度脈動,脈動幅度為0.00135 m.

大顆粒起跳的臨界條件與顆粒受力平衡以及鄰域內(nèi)小顆?？臻g排布結(jié)構(gòu)密切相關(guān),首先探查小顆?？臻g密度分布對大顆粒起跳的影響規(guī)律.選取三個(gè)關(guān)鍵時(shí)間點(diǎn),如圖3(c)中(1)、(2)和(3)點(diǎn)所示,對應(yīng)于上升過程的“起跳”起始點(diǎn),“躍升”中間點(diǎn)和“平衡”點(diǎn),考察這三個(gè)時(shí)刻大顆粒鄰域空間內(nèi)小顆粒的空間排列結(jié)構(gòu),如圖3(c)中插圖所示.發(fā)現(xiàn)在“起跳”起始點(diǎn),大顆粒上方的小顆粒排列密度較其他區(qū)域稀疏,這樣給大顆粒的“起跳”提供了上升空間; 在“躍升”中間點(diǎn),大顆粒底部小顆粒排列較頂部緊密,大顆粒似乎是被周圍顆?！皵D壓”上升的; 在“平衡”點(diǎn),大顆粒下部的小顆粒密集排列,支撐大顆粒在剪切槽頂部達(dá)到動態(tài)平衡.

圖2 實(shí)驗(yàn)與模擬結(jié)果對比圖(a)May等[10]實(shí)驗(yàn)獲得大小顆粒位置互相顛倒的現(xiàn)象;(b)本文的模擬結(jié)果Fig.2.Comparison between the present simulation and the previous experimental results:(a)The Brazil-nut effect phenomenon from May[10];(b)the present numerical simulation results.

引起大顆粒起跳的另外一個(gè)可能因素是大顆粒在高度方向的受力特征,特別需要關(guān)注大顆粒起跳前后在高度方向上的受力Fz是否具有明顯的變化,為此,我們對大顆粒進(jìn)行z方向的力學(xué)分析,提取速度Vz,受力Fz,以及轉(zhuǎn)動角速度ωz的時(shí)間序列,如圖3(d)所示,圖中的灰色條帶對應(yīng)起跳過程.從圖中可以看出,速度、受力和轉(zhuǎn)動角速度在“起跳”時(shí)刻以及“躍升”過程中,呈現(xiàn)均值為0的脈動分布狀態(tài),沒有呈現(xiàn)出直覺預(yù)期的上升力大于下沉力的規(guī)律.因此,在起跳前后,顆粒在高度方向上的受力沒有特殊性,受力特征本身不是顆粒起跳的決定因素.

圖3 摩擦系數(shù)為0.3的大顆粒運(yùn)動及力學(xué)特征圖(a)大顆粒軌跡隨時(shí)間演變圖;(b)大顆粒x,y方向位移隨時(shí)間變化圖;(c)大顆粒z方向位移隨時(shí)間演變圖;(d)大顆粒在z方向線速度、受力以及旋轉(zhuǎn)角速度隨時(shí)間變化規(guī)律,灰色條帶對應(yīng)起跳過程Fig.3.Characteristics of the large particle kinetics and kinematics with friction coefficient of 0.3:(a)Large particle trajectory evolution over time;(b)changes of large particle trajectory components in x and y directions with time;(c)changes of large particle trajectory components in z direction with time;(d)changes of translational velocity,force and rotational velocity with time.

3.2 摩擦系數(shù)對顆粒群動量傳輸以及整體堆積密度的影響

本文模擬中,顆粒群的堆積高度代表顆粒堆積的松散程度,通過大顆粒到達(dá)頂部后的平衡高度表征.如圖4(a)所示,大顆粒的平衡高度隨著摩擦系數(shù)的增大而增大,這也意味著小顆粒群空間排列更為松散,有利于大顆粒由下至上的躍升過程,對應(yīng)于大顆粒的起跳弛豫時(shí)間將會更短,這一推論由圖4(b)所示的大顆粒起跳弛豫時(shí)間隨摩擦系數(shù)的變化規(guī)律.

驗(yàn)證表明,總體上大顆粒的起跳弛豫時(shí)間隨摩擦系數(shù)的增大而減小,即隨著摩擦系數(shù)的增大,大顆粒能夠更快地達(dá)到起跳臨界點(diǎn).當(dāng)摩擦系數(shù)為0.47時(shí),起跳時(shí)間點(diǎn)較為特殊,這個(gè)特殊現(xiàn)象的原因會在后續(xù)研究中進(jìn)一步分析.

在模擬過程中,因小顆粒注入滿剪切槽的隨機(jī)性,可能出現(xiàn)大顆粒被小顆粒擊中彈起無法回落到底部而是其他小顆粒上方的情況,這種情況下大顆粒的躍升不是本文探討的內(nèi)容.本文模擬中摩擦系數(shù)為0.45和0.6兩種情況下,會出現(xiàn)以上情況,因此,在分析顆粒上升過程隨摩擦系數(shù)的變化規(guī)律時(shí),將這幾個(gè)數(shù)據(jù)點(diǎn)剔除.

圖4 (a)平衡高度隨摩擦系數(shù)變化圖;(b)起跳弛豫時(shí)間隨摩擦系數(shù)變化圖Fig.4.(a)Equilibrium height versus friction coefficient;(b)relaxation time changes with friction.

3.3 摩擦系數(shù)對大顆粒動力學(xué)特征的影響規(guī)律

圖5(a)為不同摩擦系數(shù)下z方向分力的全程變化曲線,可以看出,針對不同的摩擦系數(shù),Fz均隨時(shí)間表現(xiàn)為0值附近上下脈動的特征,且Fz在對應(yīng)的起跳時(shí)間點(diǎn)附近無明顯變化,因此顆粒的突跳并非是受力的突然變化造成的.對于不同的摩擦系數(shù),我們觀察到,隨著摩擦系數(shù)的增加,Fz上下脈動的幅值范圍變大.為了提取這種幅值變化特征,采用如下方法: 首先針對整體脈動數(shù)據(jù)去掉時(shí)均值,然后分別求出正負(fù)脈動波的均值,最后對兩均值絕對值求和即為脈動域?qū)捴?圖5(b)所示即為Fz脈動波域?qū)掚S摩擦系數(shù)的變化曲線圖,由圖可看出,隨著摩擦系數(shù)的增加,Fz域?qū)挸尸F(xiàn)增大的趨勢,即每次脈動,大顆粒獲得了更大的彈跳的能力.因此,摩擦系數(shù)的增大,有利于大顆粒在更短的時(shí)間內(nèi)實(shí)現(xiàn)起跳.

圖5 (a)不同摩擦系數(shù)下z方向受力Fz時(shí)間變化序列;(b)Fz脈動域?qū)掚S摩擦系數(shù)的變化規(guī)律Fig.5.(a)Fz changes with time for different friction coefficients;(b)Fz oscillation magnitude changes with friction coefficients.

3.4 大顆粒鄰域空間內(nèi)小顆粒群空間排列結(jié)構(gòu)

小顆粒群在空間的排列狀態(tài)對大顆粒的行為特征具有重要的影響,甚至影響大顆粒的起跳過程.Zhang和Campbell[24]提出沿用至今的思路:顆粒群體運(yùn)動行為存在兩種狀態(tài),即類液態(tài)和類固態(tài),類固態(tài)行為下,顆粒群整體呈現(xiàn)出固體特性,顆粒的相對位置以及間距變化很小; 類液態(tài)行為下,顆粒群呈現(xiàn)流體流動特征,顆粒間相對位置呈現(xiàn)較大變化.更為重要的是,這兩種狀態(tài)能夠同時(shí)存在于顆粒群體系中,中間存在一個(gè)明顯的界面劃分類液態(tài)和類固態(tài)區(qū)域.聯(lián)系本文的模擬結(jié)果,提出以下假設(shè)理論: 大顆粒在躍升過程中,處于大顆粒上部的小顆粒呈現(xiàn)類流體運(yùn)動態(tài),小顆粒間的間距變化較大,給大顆粒向上躍升提供了可能的空間;另一方面,處于大顆粒下部的小顆粒呈現(xiàn)類固體運(yùn)動態(tài),小顆粒間的間距變化較小,使得大顆粒能夠獲得底部的支撐作用,不再回落到剪切槽的底部,從而完成躍升過程,在頂部達(dá)到動態(tài)平衡.

為了驗(yàn)證以上理論,我們引入浮升因子的概念,表征大顆粒周圍球形鄰域空間內(nèi),上半部分和下半部分小顆粒個(gè)數(shù)的差異.浮升因子計(jì)算方法示意圖如圖6所示,考慮到本文模擬薄層顆粒體系,我們以大顆粒球心為原點(diǎn)、0.007 m(大顆粒半徑與小顆粒直徑之和)為半徑劃定球形鄰域.設(shè)立上下區(qū)域分界面AB,分界面穿過大顆粒球心并且平行于xy平面,分別統(tǒng)計(jì)分界面上下的小顆粒數(shù)目,定義浮升因子為上區(qū)域顆粒數(shù)目與下區(qū)域顆粒數(shù)目的比值,浮升因子越小,說明下部顆粒比上部顆粒排列越密集,越有利于大顆粒躍升.

針對低、中、高三種摩擦系數(shù),鑒別浮升因子隨時(shí)間的變化過程,其中低摩擦系數(shù)為0.35、中摩擦系數(shù)為0.49、高摩擦系數(shù)為0.55.

圖7(a)—圖7(c)分別描述了三種摩擦系數(shù)下顆粒起跳軌跡與浮升因子隨時(shí)間變化的對應(yīng)關(guān)系:起跳前,浮升因子的數(shù)值大于1,大顆粒在底部; 起跳過程中,浮升因子減小,下部小顆粒的數(shù)目大于上部小顆粒的數(shù)目,大顆粒獲得了起跳必要的上升空間,下部顆粒群提供了必要的底部支撐; 起跳結(jié)束,浮升因子小于1,大顆粒到達(dá)頂部.由圖7(a)—圖7(c)可見,大顆粒的起跳軌跡與浮升因子的突降轉(zhuǎn)變點(diǎn)具有很好的對應(yīng)關(guān)系,在起跳時(shí)刻,浮升因子的陡降,大顆粒上層的小顆粒排列疏松,為大顆粒的躍升提供至為關(guān)鍵的“窗口時(shí)間”.圖7(d)比較了浮升因子突降轉(zhuǎn)變點(diǎn)隨摩擦系數(shù)增大更快出現(xiàn),發(fā)現(xiàn)隨摩擦系數(shù)的增大,大顆粒起跳弛豫時(shí)間變短,對應(yīng)的浮升因子更快地下降到平衡狀態(tài).

圖6 浮升因子計(jì)算模型Fig.6.Floating factor calculation model.

大顆粒的平均起跳過程為3 s左右,在這個(gè)時(shí)間跨度內(nèi)顆粒受力脈動為15—30次,即大顆粒受力的高速脈動特性為大顆粒提供多次躍升的力學(xué)可能性,當(dāng)且僅當(dāng)大顆粒上部和下部的小顆粒數(shù)目比減小到臨界值(即頂部小顆粒逐漸稀疏,為大顆粒上升提供可能空間,底部小顆粒逐漸緊密,為大顆粒上升提供底部支撐),大顆粒才能夠開始躍升.

綜合來看,大顆粒起跳的控制因素是高度方向受力的高頻率脈動和起跳過程中浮升因子逐漸減小共同作用的結(jié)果.摩擦系數(shù)的增加,使顆粒群的運(yùn)動加劇,整體顆粒密度呈下降趨勢,即顆粒群高度上升,由0.0075 m升至0.00817 m,因而大顆粒能夠到達(dá)更高的平衡高度; 摩擦系數(shù)增大,大顆粒受到顆粒群更大的碰撞力,由原來的低摩擦系數(shù)(0.3)受力域?qū)?.000714N到高摩擦系數(shù)(0.55)受力域?qū)?0.00199N,致使大顆粒能夠更快上升；同時(shí)，浮升因子與大顆粒運(yùn)動軌跡對應(yīng)說明大顆粒周圍空間排布對其運(yùn)動具有重要作用.

4 結(jié) 論

本文針對薄層剪切槽中的二元顆粒分離過程開展了離散元數(shù)值模擬,發(fā)現(xiàn)初始時(shí)刻位于剪切槽底層的大顆粒最終在剪切槽頂部脈動平衡.通過對大顆粒的運(yùn)動學(xué)、動力學(xué)分析,獲得了大顆粒的運(yùn)動總體特征,發(fā)現(xiàn)了大顆粒起跳的臨界條件.結(jié)論如下:

1)大顆粒躍升過程分為三個(gè)典型的階段: 起跳前的弛豫階段,快速的起跳階段和動態(tài)平衡階段;

2)顆粒在平衡位置的高度隨著摩擦系數(shù)的增加而增大,即顆粒群更為松散,顆粒群的排布松散有利于顆粒的上升.同時(shí)顆粒的受力幅度隨著摩擦系數(shù)的增加而增大,導(dǎo)致大顆粒的起跳弛豫時(shí)間隨著摩擦系數(shù)的增大而減小;

3)提出浮升因子概念,為顆粒運(yùn)動提出局域分析思想.大顆粒的起跳控制因素是顆粒受力高頻脈動和浮升因子逐漸減小綜合作用的結(jié)果.