三維漏斗中顆粒物質(zhì)堵塞問題的數(shù)值實(shí)驗(yàn)研究?

麻禮東楊光輝張晟林平田園楊磊

1)(中國科學(xué)院近代物理研究所,蘭州 730000)

2)(中國科學(xué)院大學(xué)物理學(xué)院,北京 100049)

3)(蘭州大學(xué)信息科學(xué)與工程學(xué)院,蘭州 730000)

1 引 言

顆粒物質(zhì)是由眾多離散顆粒相互作用而形成的具有內(nèi)在有機(jī)聯(lián)系的復(fù)雜系統(tǒng).自然界中單個顆粒的典型尺度在10?6—10 m范圍內(nèi),其運(yùn)動規(guī)律服從牛頓定律;整個顆粒介質(zhì)在外力或內(nèi)部應(yīng)力狀況變化時發(fā)生流動,表現(xiàn)出流體的性質(zhì),從而構(gòu)成顆粒流.顆粒物質(zhì)看似非常簡單(相互作用通常只有摩擦、碰撞和重力),但隨著研究的不斷深入,人們開始發(fā)現(xiàn)顆粒物質(zhì)表現(xiàn)出的動力學(xué)行為其實(shí)異常復(fù)雜[1?3],傳統(tǒng)的流體動力學(xué)理論和熱力學(xué)理論均不能完整地用于描述顆粒介質(zhì)的流動行為,其中一個典型的挑戰(zhàn)就是如何描述顆粒物質(zhì)的堵塞現(xiàn)象.

近年來,顆粒物質(zhì)的堵塞行為是物理學(xué)和工程領(lǐng)域研究中的一個熱點(diǎn)和難點(diǎn)問題[4?6].顆粒體系的堵塞構(gòu)形、密度變化以及顆粒物質(zhì)從容器開口流出時的堵塞行為和臨界特征等,都得到廣泛而深入的研究和普遍關(guān)注.事實(shí)上,堵塞現(xiàn)象在很多系統(tǒng)都可以發(fā)生[7],如懸浮微粒通過狹窄通道,液氦表面電子通過納米通道,第二類超導(dǎo)中的旋渦結(jié)構(gòu)以及人或動物通過窄門等.這些系統(tǒng)的堵塞行為本質(zhì)上與顆粒系統(tǒng)類似,但是至今還沒有明確的解釋.對于顆粒系統(tǒng),當(dāng)開口尺寸與顆粒的大小相差不大時,在開口的上方很容易形成穩(wěn)定的拱結(jié)構(gòu),從而阻礙了顆粒的流動,導(dǎo)致堵塞的發(fā)生.

在各種顆粒流動的研究中,重力驅(qū)動下的漏斗流中的堵塞現(xiàn)象研究得最為深入,關(guān)于堵塞的解釋已有了許多理論模型[8,9],同時對于二維漏斗是否存在一個堵塞的臨界開口尺寸是一個有爭議的問題[9?11].在實(shí)驗(yàn)中為了得到盡可能多的堵塞事例,在堵塞形成后,通常會采取振動或噴入氣流的方法破壞拱結(jié)構(gòu),使得流動重新開始,直到下一次穩(wěn)定的拱結(jié)構(gòu)的出現(xiàn).Kondic[12]對于二維錐形漏斗的模擬顯示,從流動開始到堵塞的時間t的概率分布P(t)隨t呈e指數(shù)下降的趨勢.Guariguata等[13]構(gòu)建了水平放置的、由水流驅(qū)動的顆粒堵塞研究平臺,結(jié)果發(fā)現(xiàn)水流流速對于堵塞概率的影響很小.Lin和Fang[14]研究了橢圓顆粒流中的堵塞問題,發(fā)現(xiàn)橢圓的長軸短軸之比對于堵塞概率有著明顯的影響.Longjas等[15]的研究發(fā)現(xiàn)漏斗中顆粒粒徑不均一性會使得破壞拱結(jié)構(gòu)變得困難.Kunte等[16]對于多個開口的漏斗流中的堵塞問題進(jìn)行了研究,發(fā)現(xiàn)堵塞概率明顯小于只有一個開口時的情況.Hong等[17]利用油滴做實(shí)驗(yàn)的結(jié)果暗示了當(dāng)摩擦趨于零時堵塞概率也將趨于無窮小.

Zuriguel等[4,18?22]做了一系列實(shí)驗(yàn),研究了顆粒物質(zhì)從三維圓筒中流出的堵塞行為,他們得到了數(shù)千個堵塞的事例,測量了連續(xù)兩次堵塞之間的坍塌規(guī)模s(兩次堵塞之間從開口流出的顆粒數(shù)目).實(shí)驗(yàn)發(fā)現(xiàn):1)顆粒材料特性以及表面粗糙度對坍塌規(guī)模幾乎沒有影響,但顆粒形狀對坍塌規(guī)模有顯著影響;2)平均坍塌規(guī)?！磗〉和R?Rc的倒數(shù)成冪律關(guān)系,其中R為開口尺寸,Rc為臨界開口尺寸;3)臨界堵塞半徑隨形狀改變而變化,但臨界指數(shù)幾乎不變化;4)在開口上方放置障礙物可以明顯地減小堵塞的概率.Saraf和Franklin[23]用三維不對稱錐形漏斗做堵塞實(shí)驗(yàn),發(fā)現(xiàn)不同于Zuriguel等的結(jié)論,崩塌概率P(s)和崩塌規(guī)模s之間呈冪律關(guān)系.在以上實(shí)驗(yàn)過程中,為了觸發(fā)下一次坍塌,他們用壓縮氣體來沖擊開口附近行成的拱結(jié)構(gòu).然而氣體沖擊對實(shí)驗(yàn)結(jié)果的影響有多大并沒有進(jìn)行評估.為此,本文采用數(shù)值模擬的方法,研究當(dāng)漏斗開口打開后的第一次堵塞行為,避免了實(shí)驗(yàn)上為觸發(fā)下一次坍塌而引入的干擾,并討論了漏斗開口尺寸、漏斗錐角等對坍塌規(guī)模的影響.

2 數(shù)值模擬模型

目前,對于顆粒體運(yùn)動的數(shù)值分析方法主要分為:連續(xù)介質(zhì)方法和離散元方法.本文利用自主開發(fā)的、基于離散元方法和分子動力學(xué)的多GPU(graphics processing unit)并行程序[24,25].模擬中采用Hertz-Mindlin非線性接觸模型計(jì)算顆粒間的作用力,對于兩個相互接觸的球形顆粒i和j(半徑分別為Ri和Rj;位矢為ri和rj;速度為vi和vj;角速度為ωi和ωj),它們之間的作用力由法向力Fnij和切向力Ftij構(gòu)成.在法線方向上,顆粒間作用力包括Hertz彈性力和阻尼力,即

式中,δn=Ri+Rj?rij(其中rij=ri?rj,rij=|rij|)為顆粒間的法向重疊量;nij=rij/rij為顆粒間的法向單位向量;vnij=(vij·nij)nij(其中vij=vi?vj)為顆粒間的法向相對速度;

和

分別是法向彈性系數(shù)和法向阻尼系數(shù),其中等效彈性模量Y?,等效半徑R?以及等效質(zhì)量m?分別由公式

決定;e為顆粒間碰撞的恢復(fù)系數(shù);?為泊松比.在切向上,基于Mindlin理論,顆粒間的切向力可表述為

式中

為顆粒間的切向相對速度;utij為顆粒間的切向位移,它與顆粒的接觸時間有關(guān),剛接觸時它的值為0(即utij|t=0),隨后其值由公式

(此公式的第二項(xiàng)來源于剛體在接觸點(diǎn)的轉(zhuǎn)動,確保切向位移utij總是與接觸平面垂直)決定;和

分別是切向彈性系數(shù)和切向阻尼系數(shù),其中等效切向模量G?由公式

決定.此外,通過截?cái)嗲邢蛭灰苪tij的大小來確保Mohr-Coulomb屈服條件(即|Ftij|≤|μFnij|,μ為顆粒間的摩擦系數(shù)).如果顆粒與剛性邊界接觸,則可將邊界設(shè)為半徑無限大的球體.然后在重力的作用下,球形顆粒i的運(yùn)動方程可寫為

其中mi,vi,Ii,ωi分別為顆粒i的質(zhì)量、速度、轉(zhuǎn)動慣量和角速度;Ri是從顆粒的中心指向接觸點(diǎn)的向量,它的模就是顆粒的半徑.

圖1所示為本文所用的模型.模擬過程中只考慮顆粒所受的重力、顆粒之間及顆粒與漏斗之間的摩擦力和碰撞力,不考慮空氣阻力.重力方向沿漏斗軸向豎直向下.初始時刻起,5000個顆粒隨機(jī)從一個固定水平面產(chǎn)生,然后在重力的作用下落入漏斗中.一段時間后,顆粒在漏斗中處于穩(wěn)定靜止堆積狀態(tài),去除漏斗出口的擋板后,顆粒在重力作用下開始運(yùn)動.然后統(tǒng)計(jì)從去除擋板到第一次堵塞之間流出的顆粒數(shù)目,整個過程模擬1.2×107個時間步長(即6 s,其中堆積用時1 s,時間步長見表1).為了使模擬結(jié)果具有統(tǒng)計(jì)學(xué)意義,對于每個固定的漏斗開口尺寸和漏斗錐角,重復(fù)10000次模擬,而且每次模擬中用到的初始堆積構(gòu)型都不一樣(在模擬過程中通過改變顆粒初始產(chǎn)生速度獲得).對于5000個顆粒構(gòu)成的系統(tǒng),使用一塊GPU模擬是足夠的,且計(jì)算采用雙精度浮點(diǎn)數(shù),一個時間步長將耗費(fèi)約0.1 ms計(jì)算時間,因此一個模擬過程將花費(fèi)約20 min.模型中的參數(shù)取值見表1.

圖1 數(shù)值模型 (a)三維錐形漏斗;(b)三維平底漏斗(θ=0°)Fig.1.The simulation models:(a)Three-dimensional conical hopper;(b)three-dimensional fl at hopper(θ =0°).

表1 顆粒、漏斗的模擬參數(shù)Table 1.Simulation parameters of grains and hopper.

3 結(jié)果分析與討論

3.1 坍塌規(guī)模的概率分布

根據(jù)從去除擋板到第一次堵塞之間流出的顆粒數(shù),統(tǒng)計(jì)得到平底漏斗(θ=0°)當(dāng)開口D=3d時在10000次模擬中的坍塌規(guī)模,統(tǒng)計(jì)結(jié)果(坍塌規(guī)模的概率分布)如圖2所示,橫坐標(biāo)為坍塌規(guī)模(流出的顆粒數(shù)),縱坐標(biāo)為坍塌規(guī)模出現(xiàn)的概率.從圖2中可以看出,坍塌規(guī)模大于40的概率非常小.而且在坍塌規(guī)模等于6時,概率分布出現(xiàn)了1個峰值.峰值的左邊數(shù)據(jù)似乎為冪函數(shù)增長(圖3中的插圖),但由于數(shù)據(jù)點(diǎn)太少,這個現(xiàn)象還沒有得到很好的理解,在Zuriguel等[18]的實(shí)驗(yàn)中,氣體沖擊的持續(xù)時間、氣壓以及漏斗尺寸對它都有影響.

圖2 坍塌規(guī)模的概率分布(平底漏斗開口D=3d)Fig.2.Probability distribution of collapse scale(hopper outlet D=3d).

圖3 坍塌規(guī)模的概率分布(平底漏斗開口D=3d,主圖表示峰值右邊概率分布(半對數(shù)圖);插圖表示峰值左邊概率分布(雙對數(shù)圖))Fig.3.Probability distribution of collapse scale(hopper outlet D=3d;main fi gure,log plot of probability on the right of the peak in Fig.2;inset,log-log plot of probability on the left of the peak in Fig.2).

峰值右邊的概率分布呈指數(shù)衰減形式(圖3主圖),意味著存在一個特征參數(shù)控制這個系統(tǒng).為此圖4顯示了所有開口尺寸下歸一化坍塌規(guī)模(s/〈s〉)的概率分布.從圖中可以看出,所有概率分布曲線基本上塌縮到一條單一曲線.由于峰值左邊的分布呈冪函數(shù)(不是指數(shù)形式),所以這部分沒有很好地塌縮至單一曲線.但是峰值右邊幾乎可以坍縮至相同斜率的曲線.峰值右邊數(shù)據(jù)呈指數(shù)衰減可以用如下簡單概率模型[11]解釋:假定每個顆粒從開口流出的概率為p,并與其他顆粒無關(guān)(獨(dú)立事件),那么一個顆粒發(fā)生堵塞的概率為(1?p),則坍塌規(guī)模為s的概率為有趣的是,坍塌規(guī)模的概率分布與顆粒間的相互作用力[26]非常相似,在某種程度上表明,顆粒的堵塞行為與顆粒介觀結(jié)構(gòu)(顆粒間的力鏈和力網(wǎng))的靜力學(xué)穩(wěn)定性密切有關(guān).

圖4 坍塌規(guī)模的概率分布(不同平底漏斗開口尺寸,半對數(shù)圖)Fig.4.Probability distribution of collapse scale(with various hopper outlet size,log plot).

3.2 臨界開口尺寸

值得特別注意的一個問題是:漏斗流中是否存在一個臨界開口尺寸Dc,當(dāng)開口尺寸大于臨界尺寸時不會出現(xiàn)堵塞現(xiàn)象?直觀上,這個問題很難回答.從原理上來說,即使當(dāng)漏斗開口尺寸很大時,還是可能形成非常大的拱結(jié)構(gòu)進(jìn)而發(fā)生堵塞現(xiàn)象,只是這種情況下一般要等待足夠長的時間才能觀察到.為此,圖5顯示了在平底漏斗中,平均堵塞規(guī)?！磗〉與漏斗開口尺寸D之間的變化規(guī)律.從圖中可以看出,隨著開口尺寸的增加,平均堵塞規(guī)模開始發(fā)散.為了說明這是一個相變(低于臨界開口尺寸,會出現(xiàn)堵塞現(xiàn)象;高于臨界開口尺寸,不會出現(xiàn)堵塞現(xiàn)象),需要觀察平均堵塞規(guī)?！磗〉在臨界開口附近是否滿足冪函數(shù)分布.由于不清楚確切的臨界開口尺寸的大小,模擬上會比較困難.不過,參考文獻(xiàn)[18],利用用如下冪函數(shù)去擬合數(shù)據(jù):

其中A和γ是常數(shù).擬合結(jié)果顯示在圖5中,其中A=123.8,Dc=4.75,γ=1.5.通過擬合,在本文的系統(tǒng)中,臨界開口尺寸在4.75d左右.

圖5 平底漏斗中平均坍塌規(guī)?！磗〉與開口尺寸D的關(guān)系Fig.5.Relation of averaged collapse scale 〈s〉and hopper outlet size.

3.3 漏斗錐角的影響

進(jìn)一步地,研究了漏斗錐角θ對于堵塞問題的影響,其中漏斗直徑(D0=15d)以及漏斗開口尺寸(D=3d)都固定不變.研究發(fā)現(xiàn)當(dāng)θ為15°時(圖6),坍塌規(guī)模的分布和平底漏斗(即θ=0°)區(qū)別不大,而當(dāng)θ大于45°時,坍塌規(guī)模的分布為兩個峰,一個集中在s＜3的區(qū)域,另一個峰值大約在s=9,這個峰值略大于平底漏斗時s的峰值(s=6).關(guān)于第一個峰值的存在很容易理解,錐形漏斗開口容易使得在堆積過程中開口附近就形成比較穩(wěn)定的拱結(jié)構(gòu),從而使得體系未經(jīng)流動而直接進(jìn)入了堵塞的狀態(tài).而從s＞10的概率分布來看,漏斗錐角的增大會使得s＞10的事件出現(xiàn)的頻率變高.

同樣地,在漏斗流中是否存在一個臨界漏斗錐角θ,當(dāng)錐角大于臨界值時堵塞不會發(fā)生.考慮到當(dāng)漏斗錐角接近臨界值時坍塌規(guī)模將非常大,為此我們模擬了一個顆粒體系更大的漏斗流(顆粒數(shù)為40000).研究結(jié)果發(fā)現(xiàn),當(dāng)θ≥77°時,系統(tǒng)幾乎不會堵塞(坍塌規(guī)模大于40000).顯而易見的是,當(dāng)θ=90°時(即漏斗變成一根直管),不應(yīng)該發(fā)生堵塞,我們的數(shù)值結(jié)果表明在三維漏斗中這個臨界值大約為77°.這個臨界值和文獻(xiàn)[8]中提到的二維漏斗中的臨界角度75°很接近.

圖6 不同漏斗錐角θ(漏斗開口D=3d)時,坍塌規(guī)模的概率分布Fig.6.Probability distribution of collapse scale with various hopper angle(hopper outlet D=3d).

4 結(jié) 論

利用自主開發(fā)的GPU程序,我們通過建立三維漏斗流顆粒模型,模擬了顆粒通過圓形開口過程中顆粒堵塞現(xiàn)象,研究了坍塌規(guī)模與開口尺寸之間的關(guān)系,臨界開口尺寸是否存在,漏斗錐角對于堵塞概率影響等問題.得到如下結(jié)論:

1)對于固定的漏斗開口尺寸,坍塌規(guī)模概率分布在峰值右邊滿足指數(shù)衰減分布;

2)通過歸一化,發(fā)現(xiàn)所有開口尺寸的峰值右邊的坍塌規(guī)模分布可以基本上坍縮到單一曲線上;

3)坍塌規(guī)模概率分布在峰值右邊呈指數(shù)分布的現(xiàn)象,可以用一個簡單的概率模型進(jìn)行解釋;

4)通過擬合平均坍塌規(guī)模與開口尺寸之間的關(guān)系,發(fā)現(xiàn)存在一個臨界開口尺寸,當(dāng)開口尺寸大于這個臨界尺寸時,堵塞將不再發(fā)生;

5)當(dāng)漏斗錐角θ≤60°時,θ對于平均坍塌規(guī)模影響不大,但存在一個臨界角度,當(dāng)θ超過這個臨界角度后系統(tǒng)將不再發(fā)生堵塞.

[1]Peng Y J,Zhang Z,Wang Y,et al. 2012Acta Phys.Sin.61 134501(in Chinese)[彭亞晶,張卓,王勇,等2012物理學(xué)報(bào)61 134501]

[2]Xie X M,Jiang Y M,Wang H Y,et al. 2003Acta Phys.Sin.52 2194(in Chinese)[謝曉明,蔣亦民,王煥友,等2003物理學(xué)報(bào)52 2194]

[3]Lu K Q,Hou M Y,Jiang Z H,et al. 2012Acta Phys.Sin.61 119103(in Chinese)[陸坤權(quán),厚美瑛,姜澤輝,等2012物理學(xué)報(bào)61 119103]

[4]Zuriguel I,Pugnaloni L A,Garcimartin A,Maza D 2003Phys.Rev.E68 3

[5]Thomas C C,Durian D J 2015Phys.Rev.Lett.114 17

[6]Lu K,Liu J 2004Physics33 10(in Chinese)[陸坤權(quán),劉寄星2004物理33 10]

[7]Zuriguel I,Parisi D R,Hidalgo R C,Lozano C,Janda A,Gago P A,Peralta J P,Ferrer L M,Pugnaloni L A,Clement E,Maza D,Pagonabarraga I,Garcimartin A 2014Sci.Reports4 7324

[8]To K,Lai P Y,Pak H K 2001Phys.Rev.Lett.86 1

[9]Masuda T,Nishinari K,Schadschneider A 2014Phys.Rev.Lett.112 13

[10]To K W 2005Phys.Rev.E71 6

[11]Janda A,Zuriguel I,Garcimartin A,Pugnaloni L A,Maza D 2008EPL84 4

[12]Kondic L 2014Granular Matter16 2

[13]Guariguata A,Pascall M A,Gilmer M W,Sum A K,Sloan E D,Koh C A,Wu D T 2012Phys.Rev.E86 6

[14]Lin Y J,Fang C 2016J.Mech.32 6

[15]Longjas A,Monterola C,Saloma C 2009J.Statist.Mech.Theory and Experiment2009 05006

[16]Kunte A,Doshi P,Orpe A V 2014Phys.Rev.E90 2

[17]Hong X,Kohne M,Weeks E R 2015arXivpreprint

[18]Zuriguel I,Garcimartin A,Maza D 2005Phys.Rev.E71 5

[19]Mankoc C,Garcimartin A,Zuriguel I,Maza D 2009Phys.Rev.E80 1

[20]Zuriguel I 2014Papers in Physics6 060014

[21]Zuriguel I,Janda A,Garcimartin A,Lozano C,Arevalo R,Maza D 2011Phys.Rev.Lett.107 27

[22]Lozano C,Janda A,Garcimartin A,Maza D,Zuriguel I 2012Phys.Rev.E86 3

[23]Saraf S,Franklin S V 2011Phys.Rev.E83 3

[24]Tian Y,Qi J,Lai J,Zhou Q,Yang L 2013Proceedings of the Awareness Science and Technology and Ubi-Media ComputingAizu-Wakamatsu,Japan,November 2–4,2013 p547

[25]Tian Y,Zhang S,Lin P,Yang Q,Yang G,Yang L 2017Comput.Chem.Engineer.104 231

[26]Snoeijer J H,van Hecke M,Somfai E,van Saarloos W 2003Phys.Rev.E67 3