橢球顆粒材料在水平轉筒內混合特性的超二次曲面離散元分析*

王嗣強 季順迎

(大連理工大學,工業(yè)裝備結構分析國家重點實驗室,大連 116023)

轉筒中的顆粒流廣泛存在于工業(yè)生產(chǎn)中,顆粒形狀是影響顆粒流動的重要因素.本文基于超二次曲面方程描述球體和橢球顆粒的幾何形態(tài),采用離散元方法對水平轉筒中顆粒物質的流動和混合特性進行數(shù)值分析,并與橢球混合過程的實驗結果進行對比驗證.在此基礎之上,進一步研究了圓筒轉速、顆粒填充分數(shù)和顆粒長寬比對混合率的影響規(guī)律.計算結果表明：顆粒材料的混合率隨著轉速的增加或填充分數(shù)的減小而增加.在相同轉速和填充分數(shù)下,橢球顆粒的混合率高于球形顆粒.同時,長寬比為0.75和1.5時橢球顆粒具有最高的混合率.當長寬比小于0.75時,混合率隨著長寬比的增加而增加；當長寬比大于1.5時,混合率隨著長寬比的增加而減小.此外,橢球顆粒出現(xiàn)更明顯的速度分層現(xiàn)象.顆粒長寬比改變顆粒間的接觸模式和系統(tǒng)的密集度,增加了顆粒系統(tǒng)的平動而限制了轉動,在一定程度上提高了外部能量向顆粒系統(tǒng)轉化的效率.

1 引 言

顆粒流是由大量離散固體顆粒組成,區(qū)別于傳統(tǒng)流體介質的一種復雜流動[1,2].在外界激勵下,顆粒物質表現(xiàn)出類似固態(tài)或液態(tài)的物理力學特性,并在一定條件下相互轉化,是顆粒介質類固-液相變的宏觀表現(xiàn)[3,4].顆粒間摩擦和非彈性碰撞可有效地衰減系統(tǒng)能量,從而引起非平衡相變的特殊轉化規(guī)律[5,6].顆粒在流動過程中重新排列,毗鄰顆粒相互接觸形成強度迥異的力鏈網(wǎng)絡,進而改變系統(tǒng)的流動狀態(tài)[7,8].轉筒內顆粒流的數(shù)值模擬是研究顆粒流動的重要手段之一,其混合特性的機理研究對顆粒材料在工業(yè)生產(chǎn)和設計中的應用具有參考價值[9,10].

離散單元法(discrete element method,DEM)是由Cundall 和 Strack[11]提出,并已成為模擬顆粒流物理特性并解決相關工程問題的重要工具[12].目前,水平轉筒中顆粒材料的流動和混合特性研究主要通過對混合指數(shù)、持續(xù)時間和顆粒溫度等討論不同旋轉速度、方向、尺寸以及顆粒系統(tǒng)的填充分數(shù)、彈性模量、摩擦系數(shù)等因素的影響[13-16].實驗結果表明：顆粒摩擦和圓筒轉速顯著影響轉筒內顆粒流的動力特性,并且隨著摩擦系數(shù)和轉速的增加,顆粒系統(tǒng)的動態(tài)休止角隨之增加[17].進一步考慮容器形狀對混合率的影響,結果表明在不同填充分數(shù)下4葉混合器具有更好的混合效果[18].然而,大量實驗和數(shù)值模擬主要集中于球形顆粒,而對旋轉圓筒內非球形顆粒流的研究則相對較少[19].此外,自然界或工業(yè)中普遍是由非球形顆粒組成的復雜體系,其流動過程中的動力學行為與球形顆粒具有顯著差異[20].同時,球形顆粒系統(tǒng)得到的宏觀和細觀物理力學特性難以簡單地推廣到非球形顆粒系統(tǒng)[21].

為合理地描述非球形顆粒形狀,基于粘結或鑲嵌的組合顆粒單元[22]、基于幾何拓撲的多面體單元[23]、基于Minknowski Sum理論的擴展多面體單元[24]、基于連續(xù)函數(shù)包絡的超二次曲面單元[25]、基于球諧函數(shù)擬合的星形顆粒單元[26]等不同的構造方法不斷發(fā)展和完善.通過SIPHPM方法模擬多面體顆粒在二維旋轉圓筒內的混合行為[27],結果表明：與球形顆粒相比,三角形、四邊形和六邊形顆粒具有更高的混合度.采用超二次曲面方程構造橢球和桿狀顆粒,研究不同長寬比對非球形顆粒系統(tǒng)混合過程的影響[28,29].結果表明：高長寬比顆粒的主軸方向基本平行于顆粒的流動方向.此外,采用球形鑲嵌單元和平面-柱面組合單元模擬圓柱體顆粒的流動過程并與實驗和理論結果進行對比,發(fā)現(xiàn)不同顆粒形狀的構造方法顯著影響顆粒流的動力學行為[30].因此,本文采用超二次曲面單元構造不同長寬比的橢球顆粒,并通過非線性迭代方法計算單元間的接觸力,進而模擬橢球顆粒在轉筒中的流動和混合行為.

本文基于連續(xù)函數(shù)包絡的超二次曲面方程構造不同長寬比的橢球顆粒,采用DEM建立不同顆粒形狀填充的顆粒床在水平轉筒作用下的數(shù)值模型,并通過橢球顆?；旌线^程的實驗結果[31]對數(shù)值模型進行驗證.在不同轉動速度、填充分數(shù)和顆粒長寬比下,通過對Lacey混合指數(shù)以及混合率的計算,探討顆粒形狀對顆粒材料混合特性的影響規(guī)律,并通過顆粒系統(tǒng)的平動和轉動動能揭示非球形顆粒混合特性的內在機理.

2 橢球顆粒的離散元方法

針對橢球顆粒單元的幾何特點,這里采用基于連續(xù)函數(shù)包絡的超二次曲面方程描述顆粒的幾何形態(tài),并采用相應的非線性接觸模型計算單元間的接觸作用,對旋轉圓筒內橢球顆粒流的混合過程進行數(shù)值模擬.

2.1 基于超二次曲面方程的橢球單元

超二次曲面方程是數(shù)學意義上描述非球形顆粒的普遍方法,并由二次曲面方程經(jīng)擴展得到的函數(shù)形式[32]：

式中,a ,b 和c 分別表示顆粒沿主軸方向的半軸長,n1和n2表示顆粒表面的尖銳度參數(shù).當參數(shù)滿足n1= n2= 2時,該方程可以構造球體或橢球顆粒.圖1顯示在超二次曲面方程中改變不同的半軸長得到不同的球體或橢球單元模型.

圖1 球體和橢球的超二次曲面離散元模型 (a) a = b =2c ；(b) a = b = c ；(c) a = b = 0.5cFig.1.Super-quadric discrete element model of spherical and ellipsoidal elements：(a) a = b = 2c ；(b) a = b = c ；(c) a = b = 0.5c.

2.2 橢球單元間的接觸搜索算法

由于橢球單元是通過超二次曲面方程進行構造,因此單元間的接觸判斷方法主要采用超二次曲面單元的求解方法[28].考慮將單元間最短距離的優(yōu)化問題轉化為求解四元非線性方程組進行求解,并通過牛頓迭代算法計算單元間的重疊量[33].因此,四元非線性方程組可表示為[34]

式中 X=(x,y,z)T,?F(X) 表示超二次曲面方程的一階導函數(shù),μ2表示兩個顆粒單元間的外法向平行且方向相反,下標i 和j 分別表示單元i 和j.這里,牛頓迭代公式可表示為[35]

式中X(k+1)=X(k)+dX,μ(k+1)=μ(k)+dμ,?2F(X)表示超二次曲面方程的二階導函數(shù).如果中間點X0滿足 Fi(X0)＜ 0 且 Fj(X0)＜ 0 ,則表明兩個顆粒單元發(fā)生接觸,接觸法向定義為n =?Fi(X0)/‖?Fi(X0)‖或n =-?Fj(X0)/‖?Fj(X0)‖,如圖2所示.顆粒表面點Xi和Xj滿足：Xi=X0+αn 和 Xj=X0+βn ,這里未知參數(shù)α和β可通過非線性牛頓迭代數(shù)值求解和因此,法向重疊量可表述為

圖2 超二次曲面單元間的接觸檢測Fig.2.Contact detection between super-quadric particles.

2.3 橢球單元間的非線性接觸模型

橢球顆粒間的接觸力F主要是由法向力Fn和切向力Fs組成,并考慮法向力和切向力不通過單元形心引起的力矩Mn和Ms以及由單元間發(fā)生相對轉動時滾動摩擦引起的力矩Mr[36,37].顆粒間的法向力主要由彈性力Fcn和黏滯力Fdn組成,可分別表示為：

單元間的切向接觸力主要有彈性力Fcs和黏滯力Fds組成,同時考慮基于Mohr-Coulomb摩擦定律.則切向彈性力Fcs為

式中μs為滑動摩擦系數(shù),δt為切向重疊量,δt,max=μs(2-ν)/2(1-ν)·δn.

單元間的切向黏滯力Fds可表示為

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

當單元間發(fā)生相對轉動時,由滾動摩擦引起的力矩Mr可表示為

式中,μr為滾動摩擦系數(shù),為單位相對轉速,可表示為

3 橢球顆粒材料混合過程的實驗驗證

為驗證橢球離散元計算模型的可靠性,這里分別對轉速ωr為20 r/min和40 r/min的扁長形顆?；旌线^程進行數(shù)值模擬,并與文獻[31]的實驗結果進行對比.離散元模擬中水平轉筒的內徑D0為200 mm,長L0為20 mm.橢球顆粒的函數(shù)參數(shù)滿足：2a = 2b = 5 mm,2c = 10 mm,n1= n2= 2,主要的離散元計算參數(shù)列于表1中.同時,Lacey混合指數(shù)被用于計算轉筒中顆粒材料的混合度,可表示為[39]

總計1000個橢球顆粒在圓筒容器中考慮單元的隨機位置和角度并在重力作用下實現(xiàn)堆積.當容器內顆粒無相對運動時,圓筒開始旋轉并將離散元模擬的數(shù)值結果與實驗結果[31]進行對比,如圖3(a)所示.同時,將不同轉速ωr下的Lacey混合指數(shù)進行對比,如圖3(b)和圖3(c)所示.可以發(fā)現(xiàn),隨著轉筒圈數(shù)的增加,紅色和藍色顆粒呈螺旋形進行混合,且混合指數(shù)逐漸增加并最終達到穩(wěn)定.盡管離散元計算的數(shù)值結果與實驗結果存在一定的偏差,但在一定程度上可以很好地反映橢球顆粒的流動過程和混合特性,同時表明基于超二次曲面DEM模型的有效性.

表1 橢球顆粒離散元模擬的主要計算參數(shù)Table 1.DEM simulation parameters of ellipsoids.

圖3 不同轉動圈數(shù)下顆?；旌线^程的實驗結果[31]和離散元數(shù)值結果的對比 (a) ωr= 20 r/min的流動圖案；(b) ωr=20 r/min的Lacey混合指數(shù)；(c) ωr= 40 r/min的Lacey混合指數(shù)Fig.3.Comparison of mixing process between experiment results[31]and DEM simulation results at different rotating speeds：(a) Mixing pattern at 20 r/min；(b) Lacey mixing index at 20 r/min；(c) Lacey mixing index at 40 r/min.

4 轉動速度與長寬比對混合率的影響

顆粒形狀、填充分數(shù)和圓筒轉速是影響顆粒材料混合特性的重要因素.這里,水平轉筒的內徑D0為200 mm,長L0為50 mm,并且在z 方向施加周期性邊界條件,如圖4(a)所示.采用超二次曲面方程構造不同長寬比的橢球顆粒,函數(shù)參數(shù)滿足：σ= c /a (=b)且σ∈[0.25,3.0],如圖4(b)所示.不同形態(tài)的顆粒具有相同的質量,且等體積球體半徑為2 mm,顆粒數(shù)量為10000個.顆粒間的摩擦系數(shù)μs= 0.65,顆粒與筒壁間的摩擦系數(shù)μws=0.85,阻尼系數(shù)Cn= Ct= 0.1,其余計算參數(shù)列于表1中.

4.1 顆粒材料混合過程的離散元模擬

圖5(a)和圖5(b)顯示長寬比σ= 3.0時橢球顆粒在ωr= 20 r/min和40 r/min下的混合過程.可以發(fā)現(xiàn),紅色和藍色顆粒的流動呈螺旋狀并出現(xiàn)明顯的分層現(xiàn)象.隨著轉動圈數(shù)的增加,這個現(xiàn)象逐漸消失并且紅色和藍色顆?；揪鶆虻胤植荚趫A筒內.同時,靠近圓筒的顆粒被連續(xù)提升并流動至自由表面,并在重力作用下呈現(xiàn)S形流動圖案.隨著轉動速度的增加,顆?；旌纤俣茸兛烨铱拷鼒A筒的顆粒被提升的高度增加,并最終引起更顯著的S形圖案.圖5(a),圖5(c)和圖5(d)顯示當ωr=20 r/min時長寬比σ= 3.0,1.0和0.5的橢球和球體顆粒在不同時刻下的混合過程.可以發(fā)現(xiàn),在相同轉速下,橢球顆粒相比球形顆粒更快地達到均勻混合狀態(tài).下面將進一步通過改變轉動速度和顆粒長寬比來定量分析顆粒流動性能的差異.

圖4 三維水平圓筒的離散元數(shù)值模型(a)和不同長寬比的橢球模型(b)Fig.4.Schematic diagram of three-dimensional horizontal drum simulated by DEM model (a) and examples of ellipsoids with different aspect ratios (b).

圖5 旋轉速度和顆粒長寬比對混合過程的影響 (a) ωr= 20 r/min的橢球顆粒(σ= 3.0)；(b) ωr= 40 r/min的橢球顆粒(σ=3.0)；(c) ωr= 20 r/min的球體顆粒(σ= 1.0)；(d) ωr= 20 r/min的橢球顆粒(σ= 0.5)Fig.5.The influence of rotating speed and aspect ratio on the mixing process：(a) Ellipsoids (σ= 3.0) at 20 r/min；(b) ellipsoids(σ= 3.0) at 40 r/min；(c) sphere (σ= 1.0) at 20 r/min；(d) ellipsoids (σ= 0.5) at 20 r/min.

4.2 轉動速度對混合率的影響

初始時刻下紅色和藍色顆粒完全分離,即初始混合指數(shù)I0為0.當圓筒轉動充足時間后,顆粒系統(tǒng)最終達到完全混合狀態(tài),即最終混合指數(shù)If為1.對于球體或橢球顆粒,Lacey混合指數(shù)隨時間的變化過程可以通過積分函數(shù)進行擬合,該函數(shù)可表示為[41]

式中,G(Rmt) 是積分函數(shù),可表示為G(x)=這里,Rm是混合率,可通過擬合DEM模擬數(shù)據(jù)得到.圖6(a)和圖6(b)所示為長寬比σ= 0.5和3.0的橢球顆粒的Lacey混合指數(shù)隨時間的變化過程.其中,實線表示采用 (11) 式對模擬數(shù)據(jù)進行擬合的結果.可以發(fā)現(xiàn),對于不同長寬比的橢球顆粒,Lacey混合指數(shù)最終保持穩(wěn)定并趨近于1.同時,(11) 式可以較好地反映顆粒系統(tǒng)混合指數(shù)的變化.因此,這里對球體和長寬比σ=0.5和3.0的橢球顆粒在不同轉速下的混合率進行比較,如圖6(c)所示.可以發(fā)現(xiàn),隨著轉速的增加,顆?；旌纤俾孰S之增加.在相同轉速下,橢球顆粒的混合率高于球形顆粒.同時,橢球顆粒的長寬比顯著影響顆粒系統(tǒng)的混合率.

圖6 旋轉速度和顆粒形狀對Lacey混合指數(shù)和混合率的影響 (a) 橢球顆粒(σ= 0.5)；(b) 橢球顆粒(σ= 3.0)；(c) 混合率Fig.6.The influence of rotating speed and particle shape on the Lacey mixing index and mixing rate：(a) Ellipsoids (σ= 0.5)；(b) ellipsoids (σ= 3.0)；(c) mixing rate.

4.3 顆粒填充分數(shù)對混合率的影響

顆粒填充分數(shù)Ψ是指顆粒材料總體積與轉筒體積的百分比,這里分別取圓筒內徑D0= 200,220,240,280和320 mm,進而對球體和橢球顆粒的混合率隨填充分數(shù)在8%—22%范圍內的變化關系進行離散元分析,如圖7所示.可以發(fā)現(xiàn),隨著填充分數(shù)的增加,球體和橢球顆粒的混合率都逐漸降低.這主要是由于靠近轉筒的顆粒數(shù)目隨著轉筒內徑的增加而增加,從而提高了外部能量通過顆粒與轉筒間的摩擦向顆粒系統(tǒng)轉化的效率,這在一定程度上加快了顆粒的運動和混合行為.此外,在相同填充分數(shù)下,橢球顆粒的混合率高于球形顆粒.顆粒長寬比對混合率的影響為主要因素,而填充分數(shù)對混合率的影響為次要因素.下面將進一步通過改變顆粒長寬比來定量分析顆粒材料混合性能的差異.

圖7 ωr= 30 r/min時橢球和球形顆粒的混合率隨填充分數(shù)的變化Fig.7.Mixing rate of the ellipsoid and spherical particles varies with the fill level at ωr= 30 r/min.

4.4 顆粒長寬比對混合率的影響

圖8(a)顯示ωr= 30 r/min和50 r/min時橢球顆粒的混合率隨長寬比的變化.可以發(fā)現(xiàn),球體具有最低的混合率,而橢球在長寬比σ= 0.75和1.5時具有最高的混合率.當長寬比小于0.75時,混合率隨著長寬比的增加而增加；當長寬比大于1.5時,混合率隨著長寬比的增加而減小.這主要是由于顆粒系統(tǒng)的混合過程緊密聯(lián)系于顆粒間堆積和接觸狀態(tài),緊密的面-面接觸和有序的堆積結構具有更好的混合效果[27].然而,顆粒長寬比顯著影響顆粒間的接觸模式和緊密堆積程度,如圖8(b)所示.其中,C0為橢球顆粒系統(tǒng)的初始密集度.可以發(fā)現(xiàn),球形顆粒的體積分數(shù)最低,而橢球在長寬比σ= 0.75和1.5時具有最高的體積分數(shù),并且體積分數(shù)隨長寬比的變化規(guī)律與文獻[42,43]中基本一致.顆粒長寬比減小了橢球間的空隙,并形成排列緊密的顆粒系統(tǒng),這增加了非球形顆粒間的接觸時間,并在一定程度上提高了外部能量通過顆粒間的摩擦和非彈性碰撞向顆粒系統(tǒng)轉化的效率.

圖8 ωr= 30 r/min時橢球顆粒的混合率和初始體積分數(shù)的變化 (a)混合率；(b)初始密集度Fig.8.Mixing rate and initial packing fraction under various aspect ratios at 30 r/min：(a) Mixing rate；(b) initial packing fraction.

5 顆粒材料在混合過程中的動能分布特性

5.1 顆粒系統(tǒng)的速度分布

水平轉筒通過顆粒間的摩擦和互鎖效應將外部能量轉化為顆粒系統(tǒng)自身的動能,進而驅動或加速顆粒的運動和混合.圖9顯示球體和長寬比σ=0.5和3.0的橢球在不同轉速下顆粒系統(tǒng)的速度分布.可以發(fā)現(xiàn),顆粒系統(tǒng)被分為三部分：靠近圓筒的流動層、顆粒床中間的靜態(tài)層和顆粒床表面的滑動層.同時,隨著轉速的增加或球形向非球形的轉變,這種顆粒速度的分層現(xiàn)象變得更加顯著.一般而言,與球形顆粒相比,非球形顆粒更容易被圓筒提升,這使得顆粒床呈現(xiàn)更明顯的S形狀.同時,非球形顆粒床表面的顆粒具有更高的滑動速度,其流動層和滑動層的相對面積增加,而處于中間的靜態(tài)層面積減小,這表明整個顆粒系統(tǒng)具有更高的動能.

5.2 顆粒系統(tǒng)的平動和轉動動能

圖10顯示球體和橢球顆粒在不同轉速下的平動動能Et和轉動動能Er隨時間的變化關系.可以發(fā)現(xiàn),當圓筒開始旋轉時,顆粒系統(tǒng)從初始堆積狀態(tài)到混合狀態(tài)需要吸收更多的動能,并且橢球顆粒系統(tǒng)吸收的動能高于球形顆粒系統(tǒng).這主要是由于橢球顆粒間具有更緊密的接觸模式和互鎖效應,從而形成更穩(wěn)定的堆積結構.因此,促使橢球顆粒系統(tǒng)開始流動需要更多的動能.當轉動時間大于3 s時,橢球顆粒系統(tǒng)基本進入穩(wěn)定流動和混合狀態(tài).然而,球形顆粒系統(tǒng)的平動動能隨時間的變化具有較大的波動性.這主要是由于球形顆粒間更容易發(fā)生相對滑動或滾動,這不利于外部能量向顆粒系統(tǒng)的連續(xù)轉化.此外,隨著轉速的增加,顆粒系統(tǒng)的平動和轉動動能都隨之增加.

圖9 不同轉速下球體和橢球顆粒的速度分布 (a) 橢球(σ= 0.5)；(b) 球體(σ= 1.0)；(c) 橢球(σ= 3.0)Fig.9.The velocity profiles of granular bed for differently shaped particles：(a) Ellipsoids (σ= 0.5)；(b) spheres (σ= 1.0)；(c) ellipsoids (σ= 3.0).

為進一步比較顆粒形狀對系統(tǒng)動能的影響,這里將3—10 s內顆粒系統(tǒng)的平動和轉動動能取平均值,并分別統(tǒng)計轉速為30 r/min和50 r/min下球體和橢球顆粒系統(tǒng)的平均動能Et*和Er*,如圖11所示.可以發(fā)現(xiàn),在相同轉速下,橢球顆粒系統(tǒng)的平動動能高于球形顆粒系統(tǒng),而轉動動能低于球形顆粒系統(tǒng).同時,隨著長寬比的增加或減小,系統(tǒng)的平動動能隨之增加而轉動動能隨之減小.這主要是由于橢球顆粒在混合過程中其主軸方向基本平行于流動方向[28],有利于外部能量通過顆粒間的摩擦轉化為系統(tǒng)自身的平動動能.然而,顆粒長寬比引起的互鎖效應不利于顆粒間的相對轉動,這在一定程度上降低了外部能量轉化為系統(tǒng)的轉動動能.

圖10 在不同轉速下球體和橢球顆粒的平動動能和轉動動能隨時間的變化 (a),(b) 橢球(σ= 0.5)；(c),(d) 球體(σ= 1.0)；(e),(f) 橢球(σ= 3.0)Fig.10.Translational and rotational kinetic energy at different rotating speeds for differently shaped particles：(a),(b) Ellipsoids (σ=0.5)；(c),(d) spheres (σ= 1.0)；(e),(f) ellipsoids (σ= 3.0).

圖11 ωr= 30和50 r/min時球體和橢球顆粒的平動和轉動動能 (a)平均平動動能；(b)平均轉動動能Fig.11.Translational and rotational kinetic energy at 30 and 50 r/min for spheres and ellipsoids：(a) Average translational kinetic energy；(b) average rotational kinetic energy.

6 結 論

基于超二次曲面方程構造球體和橢球顆粒,采用離散單元法對顆粒材料在水平轉筒內的流動和混合特性進行了數(shù)值分析.通過與橢球顆?；旌线^程的實驗結果進行對比驗證,表明基于超二次曲面模型的離散元方法適用于分析橢球顆粒流的混合特性.在此基礎之上,進一步研究了轉動速度、填充分數(shù)和顆粒長寬比對混合率的影響.研究結果表明：顆粒材料的混合率隨著圓筒轉速的增加或填充分數(shù)的減小而增加.同時,相比表面光滑的球形顆粒,橢球顆粒具有更高的混合率.長寬比σ= 0.75和1.5時橢球顆粒具有最高的混合率.當長寬比小于0.75時,混合率隨著長寬比的增加而增加；當長寬比大于1.5時,混合率隨著長寬比的增加而減小.顆粒長寬比顯著影響系統(tǒng)的密集度,緊密的接觸模式和有序的堆積結構導致更高的混合率.此外,外部能量通過顆粒間的摩擦和互鎖效應轉變?yōu)轭w粒系統(tǒng)自身的動能,同時橢球顆粒出現(xiàn)更明顯的速度分層現(xiàn)象.顆粒長寬比提高了外部能量向顆粒系統(tǒng)平動動能的轉化效率,而其引起的互鎖效應限制了外部能量向系統(tǒng)轉動動能的轉化.