基于離散元法的二元顆粒流剪切分離機理分析

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

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

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

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

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

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

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

1 物理模型與模擬方法

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

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

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

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

(1)

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

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

(2)

(3)

(4)

(5)

其中

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

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

2 結果與討論

2.1 軌跡

摩擦力越大，大顆粒上升越快，但在施加滾動摩擦阻力抑制顆粒滾動的同時，也會抑制顆粒分離的發(fā)生，可見摩擦力和滾動對大顆粒的上升十分重要。在前期工作中，已模擬了11種不同的摩擦系數(shù)與3種不同初始化位置的組合共計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。當摩擦系數(shù)為0.6時，初始化位置處于內壁和中間的大顆粒才能在模擬設置的時間40 s內爬升到頂部，因此選用0.6作為摩擦系數(shù)參數(shù)。

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

(a)大顆粒上升軌跡

2.2 速度和空間排列結構

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

(a)水平方向

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

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

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

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

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

2.3 上升分析

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

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

0007、1440、1988、0249、5146—小顆粒編號。圖7 大顆粒與5個小顆粒的法向重疊量變化

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

圖8 大顆粒只接觸下方1個后傳遞顆粒時的剪切力傳遞模式

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

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

圖10 編號為1440的小顆粒對大顆粒的豎直推動力

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

3 結論

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

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

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

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