李典哲, 劉 璐, 季順迎

(大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室, 遼寧 大連 116024)

0 引 言

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

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

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

本文通過(guò)組合擴(kuò)展多面體單元構(gòu)造了形態(tài)各異的顆粒材料．采用背景網(wǎng)格法計(jì)算該模型的質(zhì)量和轉(zhuǎn)動(dòng)慣量,消除了顆粒間重疊帶來(lái)的影響．將組合顆粒的接觸判斷轉(zhuǎn)化為其基本顆粒之間的接觸判斷．通過(guò)不同形態(tài)組合顆粒的堆積與卸料過(guò)程的離散元模擬,驗(yàn)證了該方法的可行性．

1 基于擴(kuò)展多面體組合單元的模型

基于組合離散元方法構(gòu)建擴(kuò)展多面體組合單元,采用背景網(wǎng)格法計(jì)算組合單元的質(zhì)量特性,并將組合單元間的接觸問(wèn)題簡(jiǎn)化為其基本單元間接觸判斷,最后對(duì)組合單元的平動(dòng)與轉(zhuǎn)動(dòng)進(jìn)行求解．

1.1 基于Minkowski和的擴(kuò)展多面體單元

Minkowski和理論最早由德國(guó)數(shù)學(xué)家Herman Minkowski提出,其在Euclid幾何空間中給定了兩個(gè)空間體A和B,將兩個(gè)空間點(diǎn)集相疊加得到新的空間體集合,其表達(dá)式為[38]

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

(1)

式中,x和y分別為空間體A和B對(duì)應(yīng)的空間坐標(biāo)．

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

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

1.2 擴(kuò)展多面體組合單元的構(gòu)造

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

圖2 不同形態(tài)的擴(kuò)展多面體組合單元Fig. 2 Multi-dilated polyhedron elements with various shapes

1.3 擴(kuò)展多面體組合單元計(jì)算參數(shù)的確定

考慮組合單元由可重疊的任意擴(kuò)展多面體顆粒組成,且基本顆粒單元之間可存在較大的重疊量．顆粒之間的重疊會(huì)導(dǎo)致質(zhì)量特性的重復(fù)計(jì)算,不宜準(zhǔn)確得到組合單元的運(yùn)動(dòng)信息．為此,本文采用背景網(wǎng)格法計(jì)算組合單元的質(zhì)量、質(zhì)心和轉(zhuǎn)動(dòng)慣量,以避免重復(fù)計(jì)算帶來(lái)的影響[37],如圖3所示．

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

通過(guò)組合單元所在空間位置確定網(wǎng)格的計(jì)算區(qū)域,并沿計(jì)算區(qū)域的3個(gè)方向進(jìn)行劃分,由此得到多個(gè)微單元．空間網(wǎng)格越多,微單元越小,計(jì)算精度越高．根據(jù)每個(gè)微單元網(wǎng)格質(zhì)心是否在組合單元內(nèi)來(lái)判斷有效微單元網(wǎng)格的數(shù)量,通過(guò)疊加有效微單元網(wǎng)格的質(zhì)量得到組合單元質(zhì)量,其表達(dá)式為

(2)

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

組合單元的質(zhì)心可寫(xiě)作

(3)

(4)

(5)

式中,Cx,Cy和Cz分別為組合單元質(zhì)心在x,y和z方向上的坐標(biāo)位置;xg,yg和zg分別為微單元網(wǎng)格在x,y和z方向上的中心坐標(biāo)位置．

組合單元的慣性張量可寫(xiě)作

(6)

(7)

(8)

(9)

(10)

(11)

1.4 組合單元的接觸判斷

接觸判斷是離散元數(shù)值模擬的一個(gè)重要環(huán)節(jié),對(duì)計(jì)算效率有很大的影響．這里采用二階多面體擴(kuò)展函數(shù)與球面函數(shù)加權(quán)求和的方法得到擴(kuò)展多面體的包絡(luò)函數(shù), 從而將擴(kuò)展多面體間的接觸問(wèn)題轉(zhuǎn)化為兩個(gè)包絡(luò)函數(shù)之間的優(yōu)化問(wèn)題．通過(guò)求解優(yōu)化問(wèn)題快速確定顆粒間的接觸法向和重疊量,有效地提高了擴(kuò)展多面體單元的接觸搜索效率[43]．該包絡(luò)函數(shù)的歸一化形式為

(12)

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

優(yōu)化模型的表達(dá)式為：

(13)

式中,fA和fB分別為兩個(gè)基本單元的包絡(luò)函數(shù)．

由包絡(luò)函數(shù)構(gòu)造的擴(kuò)展多面體顆粒需要滿足嚴(yán)格凸形的條件[44],而擴(kuò)展多面體組合單元形態(tài)往往是不規(guī)則的、凹形的．但事實(shí)上,已證明了兩個(gè)組合單元之間的接觸重疊,至少存在一對(duì)基本單元相互接觸[33]．因此,兩個(gè)組合單元之間的接觸問(wèn)題可轉(zhuǎn)化為其基本單元即擴(kuò)展多面體顆粒之間的接觸問(wèn)題．組合單元可采用與擴(kuò)展多面體顆粒相同的接觸判斷方法．圖4為兩個(gè)相互接觸的組合單元．

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

1.5 擴(kuò)展多面體組合單元的運(yùn)動(dòng)

1.5.1 組合單元合力、合力矩計(jì)算

組合單元的運(yùn)動(dòng)可分為平動(dòng)和轉(zhuǎn)動(dòng)兩部分,根據(jù)Newton第二定律列出其運(yùn)動(dòng)方程如下：

(14)

(15)

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

(16)

(17)

對(duì)于組合單元來(lái)說(shuō),合力∑F和合力矩∑M的計(jì)算是將所有基本單元的接觸力、接觸力矩進(jìn)行疊加．將作用在基本單元上的接觸力相對(duì)于組合單元質(zhì)心進(jìn)行累加．

1.5.2 組合單元的平動(dòng)與轉(zhuǎn)動(dòng)

組合單元所在空間坐標(biāo)系可分為整體坐標(biāo)系(G)和以組合單元質(zhì)心為原點(diǎn)的全局坐標(biāo)系(GL)和局部坐標(biāo)系(B)．組合單元平動(dòng)時(shí),局部坐標(biāo)系(B)隨之發(fā)生平動(dòng),其動(dòng)力學(xué)方程可寫(xiě)作

ak+1=Fk+1/m,

(18)

(19)

(20)

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

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

(21)

式中,i為基本單元的編號(hào),xi為下一時(shí)刻第i個(gè)基本單元位置信息;Δxi為第i個(gè)基本單元與其組合單元質(zhì)心的相對(duì)位置．

組合單元轉(zhuǎn)動(dòng)時(shí)其運(yùn)動(dòng)可視為局部坐標(biāo)系(B)在全局坐標(biāo)系(GL)下的運(yùn)動(dòng)．局部坐標(biāo)系下的坐標(biāo)eB與全局坐標(biāo)系下的坐標(biāo)eGL,其轉(zhuǎn)換關(guān)系可表示為

eB=R·eGL,

(22)

式中,R為轉(zhuǎn)換矩陣,可由四元數(shù)q(q0,q1,q2,q3)求得,其可寫(xiě)作

(23)

作用在組合單元上的合力矩M在局部坐標(biāo)系下可表示為MB=R·M．根據(jù)組合單元轉(zhuǎn)動(dòng)的動(dòng)力學(xué)方程可求得下一時(shí)刻局部坐標(biāo)系下組合單元的角速度,即

(24)

(25)

(26)

由式(24)—(26)求得下一時(shí)刻局部坐標(biāo)系下組合單元轉(zhuǎn)速,然后再對(duì)組合單元四元數(shù)進(jìn)行更新,即

(27)

隨后,將更新后的四元數(shù)再次代入方程 (23)中,計(jì)算得到下一時(shí)刻的轉(zhuǎn)換矩陣．通過(guò)更新后的轉(zhuǎn)化矩陣將局部坐標(biāo)系下組合單元的角速度轉(zhuǎn)化為全局坐標(biāo)系下的轉(zhuǎn)速．

2 擴(kuò)展多面體組合單元運(yùn)動(dòng)過(guò)程的離散元分析

通過(guò)與已有文獻(xiàn)的結(jié)果對(duì)比了兩類組合單元的卸料過(guò)程,驗(yàn)證了該模型的有效性．其次,模擬不同形態(tài)單元的堆積、卸料過(guò)程,并分析顆粒形狀對(duì)體積分?jǐn)?shù)、卸料流量和配位角的影響．

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

為驗(yàn)證擴(kuò)展多面體組合單元的有效性,并與已有文獻(xiàn)結(jié)果進(jìn)行對(duì)比[46]．本文選取了文獻(xiàn)[46]中的兩類單元,分別模擬其在長(zhǎng)方形平底漏斗中的卸料過(guò)程．圖5中的這兩類單元分別是由一個(gè)基本顆粒構(gòu)成的凸形三棱柱單元和由兩個(gè)基本顆粒構(gòu)成的凹形正倒錐體單元,其可看作是在凸形三棱柱單元基礎(chǔ)上將其兩個(gè)矩形面旋轉(zhuǎn)45°得到．

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

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

表1 兩種擴(kuò)展多面體組合單元的主要計(jì)算參數(shù)

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

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

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

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

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

2.2 顆粒材料堆積過(guò)程的離散元模擬

通過(guò)組合擴(kuò)展多面體顆粒構(gòu)造6種不同形態(tài)的組合單元,如圖10所示,其中L,W和H分別表示顆粒的長(zhǎng)、寬和高,這里取L=100 mm,W=30 mm,H=180 mm．采用4 000個(gè)顆粒材料進(jìn)行堆積過(guò)程的離散元模擬,其主要計(jì)算參數(shù)列于表2中．

圖10 不同形態(tài)的擴(kuò)展多面體組合單元Fig. 10 Multi-dilated polyhedrons with various shapes

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

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

圖12 顆粒形狀對(duì)堆積分?jǐn)?shù)的影響Fig. 12 Effects of particle shapes on the piling fraction

2.3 顆粒材料卸料過(guò)程的離散元模擬

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

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

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

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

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

3 結(jié) 論

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