三維漏斗中顆粒物質堵塞問題的數值實驗研究?

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

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

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

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

1 引 言

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

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

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

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

2 數值模擬模型

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

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

和

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

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

式中

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

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

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

決定.此外,通過截斷切向位移utij的大小來確保Mohr-Coulomb屈服條件(即|Ftij|≤|μFnij|,μ為顆粒間的摩擦系數).如果顆粒與剛性邊界接觸,則可將邊界設為半徑無限大的球體.然后在重力的作用下,球形顆粒i的運動方程可寫為

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

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

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

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

3 結果分析與討論

3.1 坍塌規模的概率分布

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

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

圖3 坍塌規模的概率分布(平底漏斗開口D=3d,主圖表示峰值右邊概率分布(半對數圖);插圖表示峰值左邊概率分布(雙對數圖))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).

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

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

3.2 臨界開口尺寸

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

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

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

3.3 漏斗錐角的影響

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

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

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

4 結 論

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

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

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

3)坍塌規模概率分布在峰值右邊呈指數分布的現象,可以用一個簡單的概率模型進行解釋;

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

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

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

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

[3]Lu K Q,Hou M Y,Jiang Z H,et al. 2012Acta Phys.Sin.61 119103(in Chinese)[陸坤權,厚美瑛,姜澤輝,等2012物理學報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)[陸坤權,劉寄星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