黃仕昭，鄭 麟，李玉秀，陳 穎

(1. 廣東工業(yè)大學(xué) 材料與能源學(xué)院，廣東 廣州 510006；2. 浙江大學(xué) 能源工程學(xué)院，浙江 杭州 310013)

顆粒系統(tǒng)在靜止時(shí)表現(xiàn)出固體的性質(zhì)，但是因內(nèi)部耗散效應(yīng)的存在而區(qū)別于固體。在稀疏顆粒流中，顆粒碰撞和作用時(shí)間較短，表現(xiàn)出流體的性質(zhì)；而在致密顆粒流中，顆粒持續(xù)接觸碰撞，既表現(xiàn)出固體的特征，又存在流體的性質(zhì)。在工業(yè)生產(chǎn)中，因?qū)︻w粒相互作用機(jī)理缺乏認(rèn)識(shí)而產(chǎn)生了許多問(wèn)題。例如，在材料制備過(guò)程中，壓制的材料因均勻性缺陷而不具備理想的性能；在藥物運(yùn)輸過(guò)程中，震動(dòng)導(dǎo)致藥物顆粒分離；在粉體材料的管道輸運(yùn)中，顆粒在靜電作用下聚集，引起管道阻塞。除此之外，谷倉(cāng)堆積、采礦、泥石流、雪崩、河床泥沙堆積、管道擁堵等一系列問(wèn)題的解決都需要掌握顆粒的作用機(jī)制，從而解決工業(yè)生產(chǎn)中實(shí)際問(wèn)題和應(yīng)對(duì)自然災(zāi)害的發(fā)生[1-2]。

顆粒系統(tǒng)可以分為準(zhǔn)靜態(tài)流、致密顆粒流和稀疏顆粒流。自巴西果效應(yīng)發(fā)現(xiàn)以來(lái)，研究者在準(zhǔn)靜態(tài)流和稀疏顆粒流這2個(gè)方面均取得了較好的研究成果，而致密顆粒流因相互作用復(fù)雜而尚未建立起統(tǒng)一的解釋機(jī)制。Gray[3]在關(guān)于致密顆粒流的綜述中表明，目前至少有10種可解釋的顆粒分離機(jī)制，包括自然滲透、對(duì)流、流態(tài)化、慣性、碰撞凝結(jié)、密度、空氣阻力差、軌跡分離、聚類(lèi)和有序沉降。顆粒的尺寸分離產(chǎn)生條件分為剪切誘導(dǎo)[4]和振動(dòng)激勵(lì)[5]。用于解釋振動(dòng)分離的機(jī)制無(wú)法解釋剪切分離的現(xiàn)象，而大部分用于解釋巴西果效應(yīng)的機(jī)制無(wú)法解釋反巴西果效應(yīng)。Hong等[6]、Breu等[7]分別利用定義的顆粒溫度和密度比預(yù)測(cè)反巴西果效應(yīng)的發(fā)生。有學(xué)者通過(guò)實(shí)驗(yàn)指出，反巴西果效應(yīng)依賴(lài)于間隙空氣，在真空狀態(tài)下不會(huì)發(fā)生；但是有其他學(xué)者模擬出了反巴西果效應(yīng)，因此反巴西果效應(yīng)依賴(lài)于間隙空氣的說(shuō)法存在爭(zhēng)議。離散元法(DEM)的提出和計(jì)算機(jī)技術(shù)的日益成熟有助于模擬各種復(fù)雜的顆粒流[8]和工業(yè)顆粒模型，比如，類(lèi)庫(kù)埃特顆粒流[9]、環(huán)形剪切模型、滾筒剪切模型[10-12]、斜槽流動(dòng)模型等，但是這些根據(jù)工業(yè)實(shí)際應(yīng)用模擬的模型具有局限性，基于某種模型得出的結(jié)論在其他模型中并不適用。

早期通過(guò)實(shí)驗(yàn)觀察得到的2種主流解釋機(jī)制是對(duì)流和孔隙填充。Knight等[13]觀察到顆粒系統(tǒng)在圓筒中產(chǎn)生對(duì)流，當(dāng)顆粒與壁面的摩擦力較大時(shí)，顆粒流動(dòng)方向從容器底的中部向上，再沿壁面向下，摩擦力較小時(shí)反之。Savage等[14]通過(guò)斜槽剪切流實(shí)驗(yàn)，提出了重力誘導(dǎo)、尺寸相關(guān)的空隙填充機(jī)制，以及因單個(gè)粒子上的接觸力不平衡而將粒子從自身層擠壓到相鄰層的擠壓-排出機(jī)制。對(duì)流和孔隙填充機(jī)制對(duì)顆粒流的研究具有重要的指導(dǎo)意義；但是這2種機(jī)制都是由宏觀觀察得出的結(jié)論，尚不足以揭示顆粒分離深刻的內(nèi)在機(jī)理，因此模擬技術(shù)的進(jìn)步有助于學(xué)者研究顆粒分離的微觀機(jī)制。Jing等[15]采用DEM模擬了大、小顆粒個(gè)數(shù)均等且大顆粒在下方的二元顆粒系統(tǒng)，揭示了充分摩擦與旋轉(zhuǎn)對(duì)促進(jìn)大顆粒遷移的重要性；小顆粒在沒(méi)有長(zhǎng)期接觸的情況下通過(guò)空隙進(jìn)行滲透，大顆粒在剪切作用下通過(guò)具有各向異性接觸網(wǎng)絡(luò)的擁擠鄰域爬升，并且受到持續(xù)剪切作用。

綜上所述，在致密顆粒流方面亟需能夠解釋顆粒分離的新機(jī)制，為工業(yè)實(shí)際應(yīng)用提供指導(dǎo)。采用DEM模擬分析顆粒分離的微觀機(jī)制是一種很好的研究思路。Tirapelle等[16]指出，滲透實(shí)際上并不是一種基本的偏析機(jī)制，而是空間、摩擦和慣性效應(yīng)綜合作用的結(jié)果，空間、摩擦、慣性和阻力的結(jié)合產(chǎn)生了豐富的偏析現(xiàn)象。

本文中采用Golick等[17]、May等[18]的實(shí)驗(yàn)?zāi)Ｐ?，在鄭麟等[19-20]研究的基礎(chǔ)上，將大、小顆粒個(gè)數(shù)均等的環(huán)形剪切薄層系統(tǒng)進(jìn)一步改進(jìn)為單一大顆粒在二元系統(tǒng)中的爬升模型。從空間排列結(jié)構(gòu)出發(fā)，通過(guò)分析鄰域小顆粒系統(tǒng)對(duì)大顆粒的作用，解釋大顆粒上升的機(jī)理。

1 物理模型與模擬方法

環(huán)形剪切系統(tǒng)、顆粒層結(jié)構(gòu)及大顆粒的鄰域空間排列結(jié)構(gòu)如圖1所示。環(huán)形剪切系統(tǒng)由內(nèi)、外圓壁面，旋轉(zhuǎn)的底盤(pán)和施加正壓力的頂盤(pán)構(gòu)成。模擬開(kāi)始前，在初始化坐標(biāo)位置固定1個(gè)大顆粒，然后隨機(jī)釋放15 000個(gè)小顆粒自然下落，堆積在圓槽道內(nèi)，弛豫一段時(shí)間，直到所有小顆粒大致穩(wěn)定下來(lái)，然后在模擬時(shí)間為0.60 s時(shí)解除大顆粒的固定，蓋上頂板并開(kāi)始以20 s為周期旋轉(zhuǎn)底盤(pán)。以圓槽底部圓心為坐標(biāo)原點(diǎn)，建立直角坐標(biāo)系，重力反方向?yàn)閦軸正方向，頂盤(pán)在豎直方向上的高度可變化，底盤(pán)逆時(shí)針旋轉(zhuǎn)，圓槽內(nèi)徑為0.588 m，外徑為0.634 m，模擬的時(shí)間步長(zhǎng)為2×10-6s，總步數(shù)為2×107，總模擬時(shí)間為40 s。在垂直方向，根據(jù)小顆粒直徑將顆粒層分為5層。小顆粒在底層的排列較有序，越往上越無(wú)序。大顆粒在模擬時(shí)間為1.50 s時(shí)開(kāi)始上升。

xyz—直角坐標(biāo)系。(a)環(huán)形剪切系統(tǒng)

基于DEM，采用開(kāi)源軟件LIGGGHTS中的Hertz軟球模型模擬環(huán)形剪切系統(tǒng)。Hertz軟球模型中2個(gè)顆粒間的力為

F=Fn+Ft=(knαijnij-γnvn，ij)+(ktδtij-γtvt，ij)，

(1)

式中：Fn為法向力，包括彈簧力和阻尼力；Ft為切向力，包括剪切力和阻尼力；kn為法向彈性碰撞系數(shù)；αij為第i、j個(gè)顆粒的法向重疊量；nij為顆粒i球心到顆粒j球心的單位矢量；γn為法向黏彈性阻尼系數(shù)；vn，ij為第i、j個(gè)顆粒的法向相對(duì)速度；kt為切向彈性碰撞系數(shù)；δt,ij為第i、j個(gè)顆粒的切向相對(duì)位移向量；γt為切向接觸的黏彈性阻尼系數(shù)；vt,ij為第i、j個(gè)顆粒的切向相對(duì)速度。

Fn與Ft必須滿足摩擦屈服準(zhǔn)則，即|Ft|≤μ|Fn|，其中|·|為向量的模運(yùn)算，μ為最大靜摩擦系數(shù)。kn、γn、kt、γt由顆粒材料性質(zhì)決定，即

(2)

(3)

(4)

(5)

其中

式中：Y為彈性模量；R為顆粒半徑；e為恢復(fù)系數(shù)；m為顆粒質(zhì)量；G為剪切模量；vi、vj分別為第i、j個(gè)顆粒的泊松系數(shù)；Yi、Yj分別為第i、j個(gè)顆粒的彈性模量；Ri、Rj分別為第i、j個(gè)顆粒的半徑；mi、mj分別為第i、j個(gè)顆粒的質(zhì)量。

模擬中顆粒的材料性質(zhì)與模擬條件參數(shù)如表1所示，其中初始化坐標(biāo)y的3個(gè)數(shù)值代表3種不同的工況。

2 結(jié)果與討論

2.1 軌跡

摩擦力越大，大顆粒上升越快，但在施加滾動(dòng)摩擦阻力抑制顆粒滾動(dòng)的同時(shí)，也會(huì)抑制顆粒分離的發(fā)生，可見(jiàn)摩擦力和滾動(dòng)對(duì)大顆粒的上升十分重要。在前期工作中，已模擬了11種不同的摩擦系數(shù)與3種不同初始化位置的組合共計(jì)33種工況，其中11種摩擦系數(shù)分別為0.3、0.33、0.36、0.39、0.42、0.45、0.48、0.51、0.54、0.57、0.6，顆粒-顆粒與顆粒-壁面的摩擦系數(shù)相同；3種初始化位置分別為x=0，y=0.297 0，0.305 5，0.314 0，z=0.299 0。當(dāng)摩擦系數(shù)為0.6時(shí)，初始化位置處于內(nèi)壁和中間的大顆粒才能在模擬設(shè)置的時(shí)間40 s內(nèi)爬升到頂部，因此選用0.6作為摩擦系數(shù)參數(shù)。

圖2所示為不同初始化位置的大顆粒上升軌跡及大顆粒與壁面沿徑向的相對(duì)位置。從圖2(a)中可以看出，3種不同初始化位置的大顆粒上升軌跡都呈現(xiàn)出階梯式的爬升現(xiàn)象，隨著高度的增加，弛豫時(shí)間也逐漸延長(zhǎng)，而初始化位置在中間的大顆粒首先完成了第1層的爬升。外壁的大顆粒在相同的時(shí)間內(nèi)沒(méi)有上升到頂部，因此大顆粒的上升與壁面有關(guān)。通過(guò)分析大顆粒相對(duì)于原點(diǎn)的半徑軌跡，可以觀察大顆粒與壁面沿徑向的相對(duì)位置，如圖2(b)所示。從圖2(b)中可以看出，位于中間的大顆粒圍繞初始化位置即圖中虛線上下波動(dòng)；初始位置在壁面的大顆粒始終接近壁面，位于外壁的大顆粒因離心力的作用而始終受到壁面的擠壓。由此可見(jiàn)，大顆粒初始位置緊貼底板，最終運(yùn)動(dòng)到顆粒群頂部，在轉(zhuǎn)盤(pán)的帶動(dòng)下也呈現(xiàn)出離心力的作用，說(shuō)明環(huán)形剪切系統(tǒng)模型與模擬參數(shù)設(shè)置具有可行性。

(a)大顆粒上升軌跡

2.2 速度和空間排列結(jié)構(gòu)

顆粒系統(tǒng)由圓槽底部的底板旋轉(zhuǎn)施加恒定的剪切力，再向上傳遞，有逐漸減小的趨勢(shì)。壁面也會(huì)影響槽道內(nèi)顆粒群的速度分布，處在壁面附近的顆粒很難逃離壁面運(yùn)動(dòng)到中間。顆粒群沿圓槽徑向可分為6層顆粒，按照大顆粒初始化位置，將圓槽道顆粒系統(tǒng)沿徑向分為5個(gè)圓環(huán)區(qū)域，以此為條件計(jì)算包含在區(qū)域內(nèi)顆粒的平均速度。同理，將顆粒系統(tǒng)沿豎直方向分為5個(gè)小顆粒層，計(jì)算每層的平均速度。圖3所示為顆粒系統(tǒng)中顆粒群沿水平、豎直方向的速度分布。從圖3(a)中可以看出，中間的顆粒運(yùn)動(dòng)較快，越靠近壁面的顆粒運(yùn)動(dòng)越慢，壁面阻滯了附近顆粒的運(yùn)動(dòng)，對(duì)中間顆粒的運(yùn)動(dòng)影響較小。顆粒系統(tǒng)的動(dòng)量輸入來(lái)源于底板，在相鄰層的傳遞中造成損失。從圖3(b)中可以看出，第1層顆粒運(yùn)動(dòng)比第2層顆粒運(yùn)動(dòng)快3倍多，然后逐層減小；由于第1層小顆粒群的速度遠(yuǎn)大于第2層的小顆粒群的速度，因此層與層之間的動(dòng)量傳遞依靠滑動(dòng)的形式進(jìn)行；大顆粒的初始化位置緊貼底板，直徑是小顆粒直徑的1.5倍，尺寸差異使得處于第1層的大顆粒入侵到第2層，改變了大顆粒所處位置周?chē)木植拷Y(jié)構(gòu)，大顆粒作為媒介，加大了周?chē)w粒層之間的動(dòng)量傳遞幅度，因此在大顆粒鄰域，第1層顆粒層對(duì)第2層的傳遞作用也是最大的。

(a)水平方向

在環(huán)形剪切槽模型中，大顆粒的上升除了受摩擦力的影響，還受到壁面的直接或間接作用，另外，動(dòng)量傳遞向上逐層減弱，大顆粒爬升越來(lái)越慢，因此選用初始化位置在圓槽中間底部的大顆粒來(lái)分析大顆粒上升到第1層的過(guò)程，從而使得到的結(jié)果更具有普適性。

在編號(hào)為1440的小顆粒持續(xù)接觸條件下大顆粒上升到第1層的軌跡如圖4所示。從圖中可以看出，在模擬時(shí)間為1.50～1.56 s時(shí)，大顆粒經(jīng)過(guò)一小段上升后進(jìn)入模擬時(shí)間為1.56～1.66 s的弛豫階段而緩慢上升，最終持續(xù)上升到第1層。在大顆粒上升到第1層小顆粒層之上的過(guò)程中，大顆粒與后方編號(hào)為1440的小顆粒持續(xù)接觸，因此編號(hào)為1440的小顆粒對(duì)大顆粒的上升起到推動(dòng)作用。以大顆粒質(zhì)心為球心，大、小顆粒半徑之和為條件半徑，搜索包含在此球域的小顆粒，得到大顆粒上升過(guò)程的鄰域空間排列結(jié)構(gòu)演變及接觸的顆粒層層數(shù)，如圖5所示，其中顆粒受到來(lái)自底板的剪切力方向向右。大顆粒的上升與周?chē)☆w粒的相對(duì)位置變化有關(guān)。從圖5中可以看出，在初始化完成的0.60 s時(shí)，大顆粒主要與第2層小顆粒接觸；經(jīng)過(guò)一段時(shí)間的弛豫，轉(zhuǎn)變?yōu)?.50 s時(shí)同時(shí)接觸第1、2層的小顆粒，并開(kāi)始上升過(guò)程；完成上升過(guò)程后，在1.82 s時(shí)與3層小顆粒都接觸。

圖4 在編號(hào)為1440的小顆粒持續(xù)接觸條件下大顆粒上升到第1層的軌跡

圖5 大顆粒上升過(guò)程的鄰域空間排列結(jié)構(gòu)演變及接觸的顆粒層層數(shù)

在大顆粒上升到第1層的過(guò)程中，空間排列結(jié)構(gòu)的變化有3個(gè)特點(diǎn)：1)當(dāng)大顆粒與底板接觸時(shí)，底板剪切作用帶動(dòng)整個(gè)空間排列結(jié)構(gòu)移動(dòng)，發(fā)生微小變化；2)底板剪切傳導(dǎo)和受壁面摩擦的間接作用引起的速度不同，導(dǎo)致空間排列結(jié)構(gòu)中的個(gè)別小顆粒被替換；3)大顆粒的慣性作用使前、后接觸的小顆粒不對(duì)稱(chēng)。

2.3 上升分析

大顆粒與底板接觸時(shí)，直接受到底板的剪切作用。當(dāng)大顆粒離開(kāi)底板時(shí)，剪切力經(jīng)第1層小顆粒傳導(dǎo)至大顆粒，第1層5個(gè)小顆粒與大顆粒接觸的個(gè)數(shù)及結(jié)構(gòu)演變過(guò)程如圖6所示，其中顆粒受到底板的剪切力方向向右。從圖中可以看出，大顆粒在模擬時(shí)間為1.50 s時(shí)與5個(gè)小顆粒接觸，隨后編號(hào)為1440的小顆粒始終與大顆粒接觸，另外的4個(gè)小顆粒出現(xiàn)間斷接觸，在1.80 s以后，大顆粒只與編號(hào)為1440的小顆粒接觸。圖7所示為大顆粒與5個(gè)小顆粒的法向重疊量變化。在軟球模型中，顆粒法向重疊量大小直接反映了顆粒接觸力的大小。結(jié)合圖6、7可以看出，在模擬時(shí)間為1.42～1.56 s時(shí)，大顆粒受編號(hào)為0007、0249、1440、1988的4個(gè)小顆粒的擠壓作用而上升；在模擬時(shí)間為1.58～1.66 s時(shí)，大顆粒主要受編號(hào)為1440、0249的前、后2個(gè)小顆粒的作用而維持高度不變；在模擬時(shí)間為1.68～1.94 s時(shí)，離心力的作用明顯，大顆粒在靠近外壁面一側(cè)編號(hào)為1440、5146、0007的3個(gè)小顆粒的共同作用下上升；在模擬時(shí)間為1.76～1.80 s時(shí)，大顆粒高度基本不變；在1.80 s以后，大顆粒在編號(hào)為1440的小顆粒的作用下上升，直到完成第1層的上升過(guò)程。

0007、1440、1988、0249、5146—模擬中5個(gè)小顆粒的編號(hào)。圖6 第1層5個(gè)小顆粒與大顆粒接觸的個(gè)數(shù)及結(jié)構(gòu)演變過(guò)程

0007、1440、1988、0249、5146—小顆粒編號(hào)。圖7 大顆粒與5個(gè)小顆粒的法向重疊量變化

第1層小顆粒層傳遞底板的剪切力是大顆粒受到的推動(dòng)力和剪切力的來(lái)源，第2層小顆粒層對(duì)大顆粒主要起固定作用。在慣性作用和尺寸差異的共同影響下，大顆粒只接觸下方1個(gè)后傳遞顆粒時(shí)的剪切力傳遞模式如圖8所示。大顆粒在填補(bǔ)空隙后，受到上阻擋顆粒的阻擋，同時(shí)，受到前阻滯顆粒的阻滯，推動(dòng)前方小顆粒前進(jìn)。在下層，大顆粒只接觸后傳遞顆粒，底板的剪切力通過(guò)旋轉(zhuǎn)的后傳遞顆粒傳遞給大顆粒。由于大顆粒的質(zhì)心高度大于小顆粒直徑，因此后傳遞顆粒的剪切力傳遞模式的效果更明顯。

圖8 大顆粒只接觸下方1個(gè)后傳遞顆粒時(shí)的剪切力傳遞模式

在慣性作用下，第1層后方小顆粒推動(dòng)大顆粒前進(jìn)，受到第2層前方小顆粒的阻滯，因此在第2層后方出現(xiàn)空隙，為大顆粒的上升提供了空間條件。第2層顆粒層對(duì)大顆粒的固定作用及空隙空間的演化如圖9所示，其中半透明顆粒為第1層小顆粒。從圖中可以看出，在模擬時(shí)間為1.50～1.56 s時(shí)，第2層后方出現(xiàn)明顯的空隙，對(duì)應(yīng)大顆粒的第1小段上升；在模擬時(shí)間為1.56～1.66 s時(shí)，大顆粒缺少上升空間而維持高度基本不變；在模擬時(shí)間為1.68～1.80 s時(shí)，大顆粒被第2層小顆粒固定，同樣缺少上升空隙，但是大顆粒的質(zhì)心高度大于0.004 m，超過(guò)了小顆粒的頂點(diǎn)高度，經(jīng)過(guò)前面時(shí)間段的弛豫，編號(hào)為1440的小顆粒對(duì)大顆粒的推動(dòng)作用逐漸增強(qiáng)，大顆粒實(shí)現(xiàn)了第2小段的上升；在1.80 s以后，大顆粒只受到編號(hào)為1440的小顆粒的作用，小顆粒對(duì)大顆粒的豎直力推動(dòng)明顯變大，如圖10所示，大顆粒完成最后的上升。

圖9 第2層顆粒層對(duì)大顆粒的固定作用及空隙空間的演化

圖10 編號(hào)為1440的小顆粒對(duì)大顆粒的豎直推動(dòng)力

大顆粒在第2層小顆粒層的固定作用下沿接觸的下方小顆粒滾動(dòng)表面上升過(guò)程如圖11所示。從圖中可以看出，在1.80 s以后，大顆粒被固定在第2層小顆粒中，只與下層編號(hào)為1440的小顆粒接觸。小顆粒在底板剪切作用下前進(jìn)，大顆粒則保持位置基本不變，因此大顆粒是沿著滾動(dòng)的小顆粒表面升高的。

3 結(jié)論

本文中基于DEM，采用開(kāi)源軟件LIGGGHTS，針對(duì)環(huán)形剪切流顆粒系統(tǒng)，研究了單一大顆粒在小顆粒群系統(tǒng)中分離的作用機(jī)理，探討了摩擦力大小、壁面、尺寸差異、慣性等因素對(duì)大顆?？臻g排列結(jié)構(gòu)的影響，得到如下結(jié)論：

圖11 大顆粒在第2層小顆粒層的固定作用下沿接觸的下方小顆粒滾動(dòng)表面上升過(guò)程

1)大顆粒軌跡呈現(xiàn)階梯狀爬升的過(guò)程，尺寸差異是大顆粒能夠沿著小顆粒群向上爬升的根本條件，大顆粒在上升過(guò)程中受到某接觸小顆粒的持續(xù)作用，并沿著此小顆粒表面上升，因此呈階梯狀向上爬升。

2)大顆粒受到空間條件、力學(xué)條件2種機(jī)制的共同作用而上升。①空間條件。顆粒群的運(yùn)動(dòng)速度在豎直方向呈非線性分布，底層顆粒的平均速度是第2層大顆粒速度的3倍多，動(dòng)量在由下到上的傳遞過(guò)程中發(fā)生損耗。大顆粒的尺寸效應(yīng)體現(xiàn)在對(duì)前上方小顆粒動(dòng)量傳遞模式的改變。大顆粒的尺寸入侵改變了前方上層小顆粒間的動(dòng)量傳遞模式，直接正面擠壓使得前方小顆粒具有與大顆粒相同的速度，而大顆粒后方的小顆粒以該層平均速度運(yùn)動(dòng)，因此在大顆粒后方出現(xiàn)空隙，成為大顆粒上升的通道。②力學(xué)條件。大顆粒的差異入侵引起小顆粒群的排斥，小顆粒群對(duì)小顆粒頂點(diǎn)高度范圍內(nèi)的大顆粒持續(xù)擠壓，空間條件的產(chǎn)生是該階段大顆粒上升的充分條件。當(dāng)大顆粒只接觸后下方1個(gè)后傳遞顆粒時(shí)，形成單一小顆粒提供傳動(dòng)的模式，大顆粒作為從動(dòng)顆粒沿著后傳遞顆粒表面上升。