基于擴展多面體組合單元的非規則顆粒材料離散元方法

李典哲, 劉 璐, 季順迎

(大連理工大學 工業裝備結構分析國家重點實驗室, 遼寧 大連 116024)

0 引 言

離散單元法(discrete element method)能全面反應顆粒材料的宏微觀力學特性,是預測顆粒材料力學行為的有效工具[1-3]．由于球形顆粒具有易于構造、接觸判斷簡單等優點,常被用于離散元數值模擬[4-8]．然而,自然界與工業生產中的碎礫石、泥沙和道砟等真實顆粒的形態往往都是非規則的、復雜的[9-10]．顆粒形狀對顆粒系統的動力學行為影響很大[11-14]．與球形顆粒相比,非規則顆粒表面不平整、不連續,其流動狀態也呈間歇性流動[15]．由于非規則顆粒之間容易形成互鎖效應,阻礙顆粒的相對運動,其往往具有較大的空隙率和較小的流動速度[16]．同時,在顆粒堆積過程中,非規則顆粒的體積分數較小而休止角較大[17]．從球形顆粒系統中得到的結論難以直接、準確地應用到非規則顆粒系統上．因此,在離散元數值模擬中發展非規則顆粒的理論模型是十分必要的．

為能更準確地模擬真實顆粒系統,各種非規則顆粒的構建方法不斷被提出．由二次曲面函數可構造得到橢球體單元,并且通過調整參數可得到二維圓盤、三維細長或扁平橢球體等不同形狀顆粒[18-19]．超二次曲面單元是橢球體單元的擴展,通過超二次曲面函數可構造出幾何對稱的顆粒單元,改變函數參數可獲得具有不同表面尖銳度和長寬比的圓柱體、橢球體等顆粒單元[20-22]．多面體單元是基于拓撲結構的顆粒單元,可表示尖銳的角點和棱邊,能很好地反映出自然界中顆粒材料的真實形態[23-24]．與多面體單元相比,擴展多面體單元是基于Minkowski和定義由任意多面體和擴展球體構建而成的,通過改變擴展半徑可得到具有不同表面尖銳度的顆粒單元[25-26]．其避免了純多面體難以直接通過幾何元素判斷接觸的缺點,提高了接觸搜索效率[27]．然而,這些顆粒模型只適用于構造凸形顆粒．近年來,針對凹形顆粒的離散元方法也在不斷發展．Li等用凹形函數確定心形顆粒形態,再通過網格法來確定顆粒之間的接觸點,其計算效率隨形狀函數復雜度的增加而降低[28]．王嗣強等采用球諧函數描述了任意幾何顆粒形態,運用基于水平集方法的任意形態接觸算法確定了顆粒間的接觸方向和重疊量[29]．Feng將任意顆粒形態離散成一系列三角單元,通過能量守恒接觸模型確定顆粒之間的接觸力、接觸法向[30-31]．盡管上述構造方法均能有效描述任意形態的顆粒單元,但由于接觸計算的復雜性,導致離散元數值模擬時間大大增加．

為描述任意形態顆粒材料的幾何構型,組合單元法被提出且不斷完善發展[32]．組合單元法可組合不同數目的任意基本單元,且允許顆粒之間重疊,進而構造出不同形態的顆粒．組合方法的優勢在于其接觸判斷可簡化為一系列簡單的基本顆粒之間的接觸判斷[33]．組合單元一直在不斷地完善與發展,其組成基本顆粒包括球體、橢球體、圓柱體和超二次曲面等．組合球體單元通過將一定數量的球體顆粒組合起來,構造出形狀復雜的非規則顆粒,因球體單元易于構造、接觸判斷簡單等優點而最先發展并廣泛應用[34-35]．組合橢球體單元基于二次曲面函數,將若干個橢球體顆粒組合而成[36]．與橢球體相比,組合單元更好地描述了顆粒形狀的不對稱性,更加接近真實顆粒形態．組合超二次曲面單元基于超二次曲面方程,將多個超二次曲面單元進行任意組合,且基本單元之間存在重疊量,可用于構造任意形態顆粒材料[37]．

本文通過組合擴展多面體單元構造了形態各異的顆粒材料．采用背景網格法計算該模型的質量和轉動慣量,消除了顆粒間重疊帶來的影響．將組合顆粒的接觸判斷轉化為其基本顆粒之間的接觸判斷．通過不同形態組合顆粒的堆積與卸料過程的離散元模擬,驗證了該方法的可行性．

1 基于擴展多面體組合單元的模型

基于組合離散元方法構建擴展多面體組合單元,采用背景網格法計算組合單元的質量特性,并將組合單元間的接觸問題簡化為其基本單元間接觸判斷,最后對組合單元的平動與轉動進行求解．

1.1 基于Minkowski和的擴展多面體單元

Minkowski和理論最早由德國數學家Herman Minkowski提出,其在Euclid幾何空間中給定了兩個空間體A和B,將兩個空間點集相疊加得到新的空間體集合,其表達式為[38]

A⊕B={x+y|x∈A,y∈B},

(1)

式中,x和y分別為空間體A和B對應的空間坐標．

依據Minkowski和的定義,將A和B兩個空間體設置為任意多面體和擴展球體,構建得到擴展多面體單元[39]．通過改變球體的擴展半徑可構建具有不同粒子尖銳度的擴展多面體單元,如圖1所示．

圖1 由不同擴展半徑球體與多面體構造的擴展多面體單元Fig. 1 Dilated polyhedrons composed of various dilated spheres and polyhedrons

1.2 擴展多面體組合單元的構造

以往離散元方法主要集中在構造凸形顆粒單元[40-42]．為更準確地構建形狀更為復雜的凹形顆粒,本文基于組合離散元方法將多個不同形態的擴展多面體單元組合起來,且顆粒之間可存在一定重疊量,進而構造任意形態的顆粒材料．不同顆粒形態的擴展多面體組合單元如圖2所示．

圖2 不同形態的擴展多面體組合單元Fig. 2 Multi-dilated polyhedron elements with various shapes

1.3 擴展多面體組合單元計算參數的確定

考慮組合單元由可重疊的任意擴展多面體顆粒組成,且基本顆粒單元之間可存在較大的重疊量．顆粒之間的重疊會導致質量特性的重復計算,不宜準確得到組合單元的運動信息．為此,本文采用背景網格法計算組合單元的質量、質心和轉動慣量,以避免重復計算帶來的影響[37],如圖3所示．

圖3 基于背景網格法的組合單元質量和慣性矩計算Fig. 3 Calculation of the mass and the moment of inertia for multi-dilated polyhedrons

通過組合單元所在空間位置確定網格的計算區域,并沿計算區域的3個方向進行劃分,由此得到多個微單元．空間網格越多,微單元越小,計算精度越高．根據每個微單元網格質心是否在組合單元內來判斷有效微單元網格的數量,通過疊加有效微單元網格的質量得到組合單元質量,其表達式為

(2)

式中,m為組合單元質量;ρ為密度;Nx,Ny和Nz分別為x,y和z方向上有效微單元網格的數量;lnx,lny和lnz分別為微單元在x,y和z方向上的邊長．

組合單元的質心可寫作

(3)

(4)

(5)

式中,Cx,Cy和Cz分別為組合單元質心在x,y和z方向上的坐標位置;xg,yg和zg分別為微單元網格在x,y和z方向上的中心坐標位置．

組合單元的慣性張量可寫作

(6)

(7)

(8)

(9)

(10)

(11)

1.4 組合單元的接觸判斷

接觸判斷是離散元數值模擬的一個重要環節,對計算效率有很大的影響．這里采用二階多面體擴展函數與球面函數加權求和的方法得到擴展多面體的包絡函數, 從而將擴展多面體間的接觸問題轉化為兩個包絡函數之間的優化問題．通過求解優化問題快速確定顆粒間的接觸法向和重疊量,有效地提高了擴展多面體單元的接觸搜索效率[43]．該包絡函數的歸一化形式為

(12)

式中,k為顆粒光滑度系數,R為球面函數的半徑．

優化模型的表達式為：

(13)

式中,fA和fB分別為兩個基本單元的包絡函數．

由包絡函數構造的擴展多面體顆粒需要滿足嚴格凸形的條件[44],而擴展多面體組合單元形態往往是不規則的、凹形的．但事實上,已證明了兩個組合單元之間的接觸重疊,至少存在一對基本單元相互接觸[33]．因此,兩個組合單元之間的接觸問題可轉化為其基本單元即擴展多面體顆粒之間的接觸問題．組合單元可采用與擴展多面體顆粒相同的接觸判斷方法．圖4為兩個相互接觸的組合單元．

圖4 擴展多面體組合單元間的接觸判斷Fig. 4 The contact detection between multi-dilated polyhedrons

1.5 擴展多面體組合單元的運動

1.5.1 組合單元合力、合力矩計算

組合單元的運動可分為平動和轉動兩部分,根據Newton第二定律列出其運動方程如下：

(14)

(15)

式中,m,I分別為組合單元的質量和慣性張量;dv/dt,dω/dt分別為組合單元的加速度和角加速度;∑F,∑M分別為組合單元的合力和合力矩;g為重力加速度．

(16)

(17)

對于組合單元來說,合力∑F和合力矩∑M的計算是將所有基本單元的接觸力、接觸力矩進行疊加．將作用在基本單元上的接觸力相對于組合單元質心進行累加．

1.5.2 組合單元的平動與轉動

組合單元所在空間坐標系可分為整體坐標系(G)和以組合單元質心為原點的全局坐標系(GL)和局部坐標系(B)．組合單元平動時,局部坐標系(B)隨之發生平動,其動力學方程可寫作

ak+1=Fk+1/m,

(18)

(19)

(20)

式中,F,a和v分別為組合單元的合力、加速度和速度;x為組合單元的位置信息;k,k+1分別表示當前時刻和下一時刻,其間隔用Δt表示．

組合單元是一個剛體,所有基本單元與組合單元之間的相對位置都是恒定的．因此,通過更新組合單元的坐標位置與所有基本單元和組合單元的相對位置,就可以得到所有基本單元的位置信息,可表示為

(21)

式中,i為基本單元的編號,xi為下一時刻第i個基本單元位置信息;Δxi為第i個基本單元與其組合單元質心的相對位置．

組合單元轉動時其運動可視為局部坐標系(B)在全局坐標系(GL)下的運動．局部坐標系下的坐標eB與全局坐標系下的坐標eGL,其轉換關系可表示為

eB=R·eGL,

(22)

式中,R為轉換矩陣,可由四元數q(q0,q1,q2,q3)求得,其可寫作

(23)

作用在組合單元上的合力矩M在局部坐標系下可表示為MB=R·M．根據組合單元轉動的動力學方程可求得下一時刻局部坐標系下組合單元的角速度,即

(24)

(25)

(26)

由式(24)—(26)求得下一時刻局部坐標系下組合單元轉速,然后再對組合單元四元數進行更新,即

(27)

隨后,將更新后的四元數再次代入方程 (23)中,計算得到下一時刻的轉換矩陣．通過更新后的轉化矩陣將局部坐標系下組合單元的角速度轉化為全局坐標系下的轉速．

2 擴展多面體組合單元運動過程的離散元分析

通過與已有文獻的結果對比了兩類組合單元的卸料過程,驗證了該模型的有效性．其次,模擬不同形態單元的堆積、卸料過程,并分析顆粒形狀對體積分數、卸料流量和配位角的影響．

2.1 顆粒材料在平底漏斗中卸料過程的離散元分析

為驗證擴展多面體組合單元的有效性,并與已有文獻結果進行對比[46]．本文選取了文獻[46]中的兩類單元,分別模擬其在長方形平底漏斗中的卸料過程．圖5中的這兩類單元分別是由一個基本顆粒構成的凸形三棱柱單元和由兩個基本顆粒構成的凹形正倒錐體單元,其可看作是在凸形三棱柱單元基礎上將其兩個矩形面旋轉45°得到．

圖5 凸形三棱柱單元和凹形正倒錐體單元Fig. 5 The convex triangular prism element and the concave upward-downward conical element

長方形平底漏斗幾何形狀如圖6所示,其長、寬、高分別為0.14 m,0.05 m和0.4 m,漏斗擋板尺寸為0.04 m×0.05 m．這里選取600個藍色顆粒和400個紅色顆粒,并將其劃分為5層,以便更好地觀察顆粒的堆積和卸料過程．表1給出了數值模擬中的主要計算參數,包括顆粒的特定參數和一般模擬參數[46]．初始時漏斗中放置擋板,所有顆粒具有隨機位置和空間方位,在重力作用下自由下落堆積．當顆粒堆積穩定后撤去擋板,完成漏斗卸料過程．

表1 兩種擴展多面體組合單元的主要計算參數

圖6 平底漏斗幾何形狀(單位： m)Fig. 6 The geometry shape of the flat bottom hopper (unit： m)

將兩種組合單元模擬平底卸料過程所得結果與相關試驗以及離散元模擬結果[46]進行比較分析,如圖7、圖8所示．在三棱柱單元的模擬中,組合單元的流動從中間向邊緣發展,組合單元逐漸呈V型流動．這主要是由于漏斗中心的顆粒速度要高于兩側邊壁的顆粒速度．已有的試驗和模擬中分別在1.6 s和2.4 s時形成準穩定拱門,但在大多數情況下這些準穩定拱門最終會坍塌[46]．而在本文的模擬中,首先在1.6 s處觀察到準穩定拱門的形成,而后拱門坍塌,漏斗繼續卸料,2.4 s時再次形成準穩定拱門．這主要是因為凸形顆粒的流動是連續的,也更容易滑動和轉動,不能像凹形顆粒一樣互鎖形成穩定的拱形結構．而對于凹形正倒錐體單元,在已有的試驗、模擬以及本文的數值模擬中均在0.8s處形成穩定的拱形結構,且在整個模擬過程中拱形結構保持穩定．這主要是因為凹形顆粒之間的互鎖結構阻礙了組合顆粒的連續流動．其次,靠近邊壁的凹形顆粒更容易形成互鎖效應,從而形成更為穩定的拱形結構．

圖7 凸形三棱柱單元卸料過程的離散元模擬與試驗對比Fig. 7 The convex triangular prism element’s discharge process simulated with the DEM and compared with the physical experimental results

為進一步驗證本文離散元方法的可靠性,下面將卸料過程的流量進行對比分析．文獻[46]共進行了20次試驗和5次數值模擬,其中20次試驗的結果存在一定差異,這主要是由于顆粒材料的初始排列狀態具有很強的隨機分布規律,導致最終堆積狀態不同,這會對卸料結果產生一定影響[46]．圖9為平底漏斗卸料過程中剩余顆粒比例隨時間的變化情況,在此給出了文獻[46]20次試驗結果的上下限和一次數值模擬的結果．可以看到,對于凸形三棱柱單元,本文模擬結果在文獻[46]試驗結果的區域范圍內,且與離散元模擬結果在趨勢上是一致的．此外,本文模擬結果與文獻[46]的模擬結果最終剩余顆粒比例相近,不同的是開始時本文模擬漏斗卸料速度更快．這主要是因為在本次模擬中我們選擇的基本單元為擴展多面體單元,即使給定的擴展半徑很小,但其較多面體而言,顆粒表面更為光滑,顆粒也更易流動．

對于凹形正倒錐體單元,本文模擬所得結果與文獻[46]的試驗、離散元模擬結果在趨勢上一致,且其位于試驗結果的區域范圍內．顯然,在最終穩定狀態下平底漏斗中剩余顆粒比例與已有試驗結果上限相近．從整體上看,本文模擬結果與文獻[46]的試驗和離散元模擬結果是一致的．

2.2 顆粒材料堆積過程的離散元模擬

通過組合擴展多面體顆粒構造6種不同形態的組合單元,如圖10所示,其中L,W和H分別表示顆粒的長、寬和高,這里取L=100 mm,W=30 mm,H=180 mm．采用4 000個顆粒材料進行堆積過程的離散元模擬,其主要計算參數列于表2中．

圖10 不同形態的擴展多面體組合單元Fig. 10 Multi-dilated polyhedrons with various shapes

長方體容器的長、寬和高分別為L0=1.8 m,W0=1.8 m和H0=5.0 m．初始時刻,所有組合單元具有隨機位置和空間方位,且組合單元之間無重疊,在重力作用下下落堆積．下落一段時間后,組合單元保持靜止,形成穩定的顆粒床,如圖11所示,其中顆粒顏色表示其距離底部的距離．顆粒形狀對堆積分數的影響如圖12所示．不同顆粒之間堆積密度大不相同,其中O形顆粒較其他顆粒具有更高的堆積密度和更低的孔隙率．這主要是由于O形構造使得其空心部分位于顆粒中間,難以與其他顆粒勾結互鎖,顆粒更容易發生相對滑動和轉動．其次,T形、Z形和K形相比于A形和C形有更高的堆積密度和更低的孔隙率．這主要是由于A形、C形相較于T形、Z形和K形的空心部分更大,其更容易與多個顆粒產生顯著的互鎖效應,使得顆粒之間難以發生相對運動．由此可見,顆粒形狀對顆粒床的堆積密度有顯著影響,且隨顆粒形狀復雜度的增加,顆粒的體積分數不斷降低．

圖11 擴展多面體組合單元形成的穩定顆粒床Fig. 11 Stable granular beds composed of multi-dilated polyhedrons

圖12 顆粒形狀對堆積分數的影響Fig. 12 Effects of particle shapes on the piling fraction

2.3 顆粒材料卸料過程的離散元模擬

通過不同擴展多面體顆粒可構造形態各異的組合單元,其顆粒形態將影響流動行為．這里采用擴展多面體的不同組合形式構造4種組合單元,組合單元總數為4 000,其長、寬和高均相同,分別為L=1 m,W=1 m,H=1 m．初始時刻,組合顆粒在上部漏斗中進行堆積,形成穩定床后,采用3種顏色按照高度對顆粒床均勻地劃分為3層(藍色、綠色和紅色)．顆粒床堆積穩定后,將上部漏斗的孔口打開,顆粒開始卸料．圖13顯示了不同時刻顆粒的卸料過程．顆粒尺寸較大且凹形顆粒之間易形成互鎖結構阻礙顆粒流動,因此,將上部漏斗的孔口尺寸設置偏大,L1=W1=8 m．其中,L0∶W0∶H0=2∶2∶1,如圖14所示．漏斗孔口打開后底層藍色顆粒全部流出,顆粒逐漸呈V型流動,這主要是由于漏斗中心的顆粒速度高于靠近邊壁的顆粒速度．凹形顆粒之間易形成互鎖結構,顆粒不易滑動和滾動,因此在卸料結束后顆粒出現明顯的分層圖案．

圖13 不同時刻下組合單元的卸料過程Fig. 13 Discharging processes of multi-dilated polyhedrons at different moments

圖14 卸料漏斗的幾何尺寸Fig. 14 The geometry shape of the hopper model

不同顆粒形狀所對應的卸料流量和休止角影響如圖15所示．在卸料過程中,空心體組合單元具有最快的流動速率且最先完成卸料過程,而扭王字組合單元具有最慢的流動速率．卸料完成后,組合單元在漏斗底部保持靜止,形成穩定的堆積床．其中,空心體組合單元具有最小的休止角,而扭王字組合單元具有最大的休止角．這主要是因為顆粒材料形狀的復雜度顯著影響顆粒的堆積和卸料過程,凹形顆粒形狀越復雜,顆粒之間的互鎖效應越明顯,顆粒之間越不容易滑動和相對轉動．

(a) 卸料流量 (b) 休止角(a) The mass flow rate (b) The angle of repose

3 結 論

本文以Minkowski和方法構造的擴展多面體顆粒為基本單元,基于組合單元方法提出了一種任意形態擴展多面體組合模型．為驗證該模型的準確性,分別模擬了凸形三棱柱單元、凹形正倒錐體單元在平底漏斗中的卸料過程,并與已有試驗結果進行了對比驗證．基于組合方法構造了形狀各異的顆粒材料,并進行了試驗驗證,包括下落堆積和動態漏斗卸料．此外,本文還研究了顆粒形狀對堆積密度、卸料流量和休止角的影響．計算結果表明,組合單元形狀的復雜度顯著影響顆粒材料的堆積與卸料,并且隨著組合單元形狀復雜度的增加,顆粒的堆積密度與流動速率隨之降低,而休止角隨之增加．組合單元模型的有效應用,為任意形態顆粒材料的離散元數值模擬提供了一種新思路．