基于隨機游走模型的筒倉內(nèi)顆粒流堵塞概率測量

王海濤，陳 泉，李 然，修文正，楊 暉

（1．上海理工大學(xué) 光電信息與計算機工程學(xué)院，上海 200093；2．上海理工大學(xué) 健康科學(xué)與工程學(xué)院，上海 200093）

引 言

筒倉是一種常見的顆粒物料存儲設(shè)備，廣泛應(yīng)用于工業(yè)生產(chǎn)和日常生活中[1-4]。理想的筒倉能夠在最小的占地面積上存儲最多的產(chǎn)品數(shù)量，并且以所需的流量正常流出。然而，當(dāng)顆粒流至筒倉出口附近時，由于筒倉壁間的距離急劇減小，顆粒之間的相互作用增強，顆粒間很容易形成較長的力鏈或力拱，導(dǎo)致筒倉顆粒流中斷，形成堵塞現(xiàn)象。該現(xiàn)象不僅會產(chǎn)生安全隱患，更會給工農(nóng)業(yè)生產(chǎn)帶來巨大損失。因此，預(yù)測筒倉內(nèi)顆粒流的堵塞概率，探究其變化規(guī)律，一直受到人們的廣泛關(guān)注[5-10]。

近年來，大量研究人員通過不同的實驗方法研究了筒倉內(nèi)顆粒流的堵塞概率[11-17]。Goldberg等[11]通過離散元法模擬了正多邊形和圓形顆粒在筒倉中的運動，發(fā)現(xiàn)正多邊形顆粒上頂點的數(shù)目越少，堵塞概率越大，并且堵塞概率與顆粒的頂點數(shù)目呈非線性關(guān)系。López-Rodríguez等[12]采用相機對不銹鋼球體在筒倉開口上方形成的堵塞拱進行拍攝，通過測量崩塌規(guī)模(〈s〉)，觀察到隨著漏斗壁傾斜角度的增大，堵塞概率在逐漸減小。尤其，當(dāng)漏斗壁傾斜角度從0°增加到80°時，堵塞的概率會減小3個數(shù)量級。Gella等[13]采用面陣CCD相機拍攝筒倉內(nèi)不銹鋼球體兩次連續(xù)堵塞中的流動過程，發(fā)現(xiàn)筒倉寬度增大到顆粒直徑的15～20倍過程中，筒倉中的堵塞概率在不斷增大，繼續(xù)增大筒倉的寬度后，堵塞概率將不再發(fā)生變化。To等[14]首次提出了隨機游走模型，并成功應(yīng)用于二維無摩擦漏斗系統(tǒng)中堵塞概率的預(yù)測。該模型根據(jù)堵塞拱的特征參數(shù)，計算出堵塞概率的預(yù)測值，將預(yù)測結(jié)果與實驗結(jié)果進行比較，發(fā)現(xiàn)兩者基本一致。在To等的研究中，該模型要求的條件過于嚴格，只能應(yīng)用在近似無摩擦的二維漏斗中。對于有摩擦的矩形二維筒倉系統(tǒng)，文中并沒有對該模型的適用范圍作出說明[18]。

本文在上述研究內(nèi)容的基礎(chǔ)上，針對矩形筒倉，將圖像法應(yīng)用于堵塞拱的特征參數(shù)檢測，同時結(jié)合隨機游走模型計算出理論堵塞概率，重點探究在不同流量下理論堵塞概率和實驗堵塞概率之間的關(guān)系。

1 實 驗

1.1 實驗裝置

圖1為實驗裝置示意圖。本實驗裝置系統(tǒng)主要由矩形筒倉、亞克力擋板、相機和LED光源組成。矩形筒倉由有機玻璃制作而成，筒倉壁厚10 mm，高550 mm，左右壁間距150 mm，前后壁間距8 mm，筒倉孔口是位于筒倉底部的一個矩形槽，其長度8 mm，寬度R可調(diào)節(jié)，R∈[0, 8 mm]。將相機（GE1050型號，分辨率為1 024×1 024）與筒倉放置在同一水平面上，通過對筒倉孔口附近50 mm×50 mm的正方形區(qū)域進行拍照,獲得顆粒堵塞后的正視圖[19]。LED光源型號為SL-200W，輸出功率為250 W，其光照正對孔口處，讓反光點盡量顯示在正視圖中顆粒圓形輪廓的中心處。另外，在實驗過程中，通過在筒倉內(nèi)插入不同寬度的亞克力擋板，從而調(diào)整矩形筒倉的寬度(12 mm、20 mm、28 mm、32 mm、38 mm、54 mm、86 mm、120 mm)。圖1中插圖顯示的是填充顆粒黑色電鍍玻璃珠，表1為實驗顆粒的特征參數(shù)。其中摩擦系數(shù)0.58是通過計算錐形顆粒堆休止角的正切值得到的，其代表的是顆粒與顆粒之間的摩擦力和作用在其表面上的垂直力之比值。實驗過程中采用分布式填充方式進行顆粒加載[20]，填充高度為550 mm[21]。

圖 1 筒倉內(nèi)顆粒流堵塞概率的實驗裝置示意圖Fig. 1 Schematic diagram of the experimental system used for calculating the jamming probability of granular flow in a silo

表 1 顆粒的特征參數(shù)Tab. 1 Properties of the granular materials used in this study

1.2 堵塞拱的特征提取

Fig. 2 Schematic diagram of feature extraction of jamming arch

圖2是顆粒堵塞拱的處理過程。當(dāng)筒倉內(nèi)的顆粒流堵塞后，采用面陣CCD相機拍攝形成顆粒堵塞拱的正視圖，將相機拍攝得到的圖像裁剪成適當(dāng)尺寸，使組成堵塞拱的顆粒都包含在圖像中，如圖2（a）所示。對于原始彩色圖像，首先將圖像轉(zhuǎn)換成灰度圖像，在灰度化處理過程中采用的是加權(quán)平均法，對彩色圖像中R、G、B分量進行加權(quán)平均。其次，對灰度圖像進行濾波處理，在濾波處理過程中采用的是雙邊濾波的方法，先用雙邊濾波器對灰度圖像進行濾波得到低頻部分，再用灰度圖和低頻部分作差后得到高頻分量，最后將高頻分量和低頻分量分別增強后再進行合成，得到濾波后的灰度圖像。該方法的特點是保邊去噪，相較于其他濾波，在平滑圖像的同時，增強了對圖像邊緣的保護。最后，對濾波后的圖像進行全局閾值分割，通過使用最大類間方差法找到輸入圖像合適的閾值，再根據(jù)這個閾值將灰度圖像轉(zhuǎn)換成二值化圖像。為了畫出堵塞拱的結(jié)構(gòu)圖，必須求出二值化圖像中每一個反光區(qū)域的中心像素坐標[22-23]。

首先，以圖像左下角的像素點作為坐標原點，以水平方向為橫軸、豎直方向為縱軸，相鄰兩個像素點間的距離為一個刻度，建立二維平面坐標系，如圖2（b）所示。分別將每個反光區(qū)域中所有像素點橫軸和縱軸方向最小值加上最大值的一半作為中心點坐標值（ c enterX[i] ， c enterY[i] ）：

式中：X[i] 是從左向右數(shù)第i個反光區(qū)域中存放所有像素點橫坐標的數(shù)組；Y[i] 是從左向右數(shù)第i個反光區(qū)域中存放所有像素點縱坐標的數(shù)組。將顆粒堵塞拱從左向右依次排序，用箭頭表示有向線段并依次連接反光區(qū)域的中心像素坐標，得到圖2（b）所示堵塞拱結(jié)構(gòu)示意圖。由于已知每個中心點的像素坐標值，再通過斜率公式和反切函數(shù)計算得到每一個有向線段與水平正方向夾角 θi：

式中： c enterY[i+1] 是第i+1個反光區(qū)域中心點縱軸坐標； c enterY[i] 是第i個反光區(qū)域中心點縱軸坐標； c enterX[i+1] 是第i+1個反光區(qū)域中心點橫軸坐標； c enterX[i] 是第i個反光區(qū)域中心點橫軸坐標。用ri表示第i個反光區(qū)域中心到第i+1個反光區(qū)域中心的方向向量， θi是堵塞拱上第i個顆粒的ri與水平正方向的夾角，X是堵塞拱上第一個顆粒到最后一個顆粒的方向向量的水平分量，D是顆粒直徑，R是筒倉孔口寬度。如圖2（c）所示，ri的長度始終等于D， θi∈[-180°,180°]，X + D總是大于孔口尺寸R[24-25]。

圖 3 堵塞概率J隨實驗次數(shù)N的變化散點圖Fig. 3 Scatter diagram of jamming probability J and number of experiments N

1.3 堵塞概率計算

筒倉內(nèi)顆粒流形成穩(wěn)定堵塞拱的過程中，離散顆粒組成的堵塞拱結(jié)構(gòu)軌跡是滿足隨機游走模型的。該模型最早由卡爾-皮爾遜于1905年在《自然》雜志中提出[26]，2001年To將該模型應(yīng)用到二維漏斗的堵塞中[14]。根據(jù)堵塞拱的幾何特征參數(shù)，堵塞拱的概率分布函數(shù)（an(x)）可以表示為

式中：n為堵塞拱上的顆粒個數(shù)；x=X/D；βn-1=max{-π/2,θn-2-2π/3} ；An是歸一化系數(shù)；δ (X) 為狄拉克函數(shù)。對式（3）進一步處理可以將an(x)簡化成如下公式：

當(dāng)筒倉孔口寬度為R時，對式（4）進行積分運算獲得堵塞拱上顆粒個數(shù)為n時的概率（jn(R) ），其計算公式為

式中D為顆粒直徑。因此，卸料過程中，筒倉內(nèi)顆粒流發(fā)生堵塞事件的概率為

式中g(shù)R(n) 是由n個顆粒組成的堵塞拱發(fā)生的次數(shù)占總堵塞事件次數(shù)的比值。式(6)的結(jié)果就是理論堵塞概率J1的計算結(jié)果。

2 結(jié)果與討論

本實驗中，改變筒倉寬度后，加載滿整個筒倉的顆?？倲?shù)量發(fā)生了改變。為了保證不同筒倉寬度條件下實驗結(jié)果的可對比性，必須選擇流出相同質(zhì)量的顆粒作為一次實驗。同時，也要滿足在相同條件下，不同實驗次數(shù)中筒倉孔口處流出的流量保持不變。因此，選擇筒倉寬度為最小值（W=12 mm）時，在顆粒從滿載流到距筒倉底部為兩倍筒倉寬度位置處時，記錄流出的顆粒質(zhì)量為80 g±2 g。因此，在所有的實驗中，以流出80 g顆粒質(zhì)量作為一次實驗。筒倉內(nèi)顆粒的應(yīng)力重定向到側(cè)壁而引起筒倉底部壓力迅速飽和，從而導(dǎo)致顆粒從孔口流出的流量幾乎不變，這種現(xiàn)象必須滿足筒倉內(nèi)顆粒的高度是筒倉寬度的兩倍以上[27]。

2.1 堵塞概率的統(tǒng)計

圖3為堵塞概率J與實驗次數(shù)N的變化散點圖。在本實驗中，使用圖1所示裝置，平底筒倉寬度為38 mm，孔口寬度為7 mm，實驗的總次數(shù)N=1 325次，顆粒流發(fā)生堵塞的實驗次數(shù)Na=106次。根據(jù)以下公式計算出實驗堵塞概率J2=8%。從圖中可以看出，當(dāng)實驗次數(shù)較少時，堵塞概率的變化幅值較大。這是由于筒倉內(nèi)顆粒流的堵塞事件是隨機發(fā)生的。當(dāng)實驗次數(shù)大于400次時，堵塞概率J的值趨于穩(wěn)定。因此，在接下來的不同實驗中，堵塞實驗的總次數(shù)都設(shè)定為400次。

2.2 堵塞拱的特征值

圖4（a）為2.1節(jié)中106次堵塞事件的堵塞拱上顆粒個數(shù)n的統(tǒng)計分布圖?？梢钥闯?，堵塞拱上顆粒個數(shù)n分布在[4,16]范圍內(nèi)。當(dāng)11＜n＜16時，Na隨n的增大而減小。而當(dāng)n＞16時，Na=0。另一方面，對具有相同顆粒個數(shù)n的堵塞拱上 θi的值進行算數(shù)平均計算，得到 μθi與顆粒位置i的分布，如圖4（b）所示。在誤差小于5%的條件下，采用最小二乘法對測量數(shù)據(jù)進行擬合處理，擬合曲線如圖中虛線所示。可以看出，當(dāng)n∈[4,11] 范圍內(nèi)時， μθi與i呈負線性關(guān)系。這說明了堵塞拱的形狀為凸形結(jié)構(gòu)。而當(dāng)n∈[12,16] 范圍內(nèi)時， μθi與i不再滿足線性關(guān)系，如圖4（c）所示。

圖 4 堵塞拱的特征圖Fig. 4 Characteristic diagram of the jamming arch

進一步，將堵塞拱的特征參數(shù)n與 θi代入式（3）～（6）中，分別計算n∈[4,11] 、n∈[12,16] 和n∈[4,16]的理論堵塞概率J1，結(jié)果如表2所示，表2中還列出了通過實驗統(tǒng)計得到的實驗堵塞概率J2?？梢钥闯?，當(dāng)n∈[4,11] 時，理論堵塞概率J1=7.7%接近于實驗堵塞概率8%。而當(dāng)n∈[12,16] 和n∈[4,16] 時，理論堵塞概率J1遠大于實驗堵塞概率J2。因此，結(jié)合圖4（b）和圖4（c）可以認為，在有摩擦力的筒倉中，相同筒倉寬度和相同孔口寬度下，選擇堵塞拱的特征參數(shù)μθi與i呈線性關(guān)系，即選擇堵塞拱為凸性拱時計算得到的理論堵塞概率要比選擇非凸性拱計算得到的理論堵塞概率更接近實驗堵塞概率。

表 2 筒倉內(nèi)顆粒流的堵塞概率Tab. 2 Jamming probability of granular flow in silo

2.3 流量對堵塞概率的影響

由2.2節(jié)結(jié)果可知，在相同筒倉寬度和相同孔口寬度條件下，選擇凸性拱計算得到的理論堵塞概率更接近實驗堵塞概率。本節(jié)中，將探究通過改變筒倉寬度和孔口寬度，選擇凸性拱后計算得到的理論堵塞概率與實驗堵塞概率是如何變化的。同時，再進一步討論筒倉寬度和孔口寬度發(fā)生改變后，顆粒流出孔口的流量與堵塞概率之間的關(guān)系。

圖5（a）為流量與堵塞概率隨筒倉孔口寬度的變化散點圖?？梢钥闯?，流量Q與孔口寬度R呈正相關(guān)關(guān)系，這與Beverloo方程所描述的流量變化趨勢是一致的[28-29]。然而，堵塞概率J隨R的增大而減小，即堵塞概率J與流量Q呈負相關(guān)關(guān)系，Gella在研究筒倉的堵塞與顆粒材料的運動學(xué)關(guān)系時也得出過相似的結(jié)論[13]。另一方面，在流量較小時，理論堵塞概率與實驗堵塞概率存在較大偏差；而當(dāng)流量Q＞8.3 g/s（孔口寬度R＞7 mm）時，理論堵塞概率與實驗堵塞概率的差值δ＜0.3%。此時，隨機游走模型計算得到的理論堵塞概率能夠準確預(yù)測出實驗堵塞概率。

圖 5 流量與堵塞概率的散點圖Fig. 5 Scatter plot of mass discharge rate and jamming probability

圖5（b）為流量和堵塞概率隨筒倉寬度的變化散點圖。可以看出，當(dāng)12 mm＜W＜38 mm時，流量Q隨筒倉寬度W的減小而迅速增大，最大流量為11.2 g/s；而當(dāng)W＞54 mm時，流量基本保持不變，Q= 8.3 g/s±0.1 g/s，這與Hirshfeld的研究結(jié)論相一致[30]。進一步，通過改變筒倉寬度增大流量，堵塞概率J與流量Q仍然滿足負相關(guān)關(guān)系。值得注意的是，當(dāng)流量Q＞8.3 g/s（筒倉寬度W＜38 mm）時，理論堵塞概率與實驗堵塞概率之間的差值δ＜0. 3%。此時，隨機游走模型計算得到的理論堵塞概率能夠準確預(yù)測出實驗堵塞概率。

3 結(jié) 論

本文在顆粒間存在摩擦的實驗裝置內(nèi)，通過改變筒倉孔口寬度和筒倉寬度，選用均值粒徑為1.85 mm的球形玻璃顆粒，利用圖像法測量了堵塞拱結(jié)構(gòu)特征參數(shù)，根據(jù)隨機游走模型，計算出了顆粒的理論堵塞概率，并分析了筒倉中堵塞概率的變化情況，得到如下結(jié)論：在有摩擦力的筒倉中，筒倉寬度和孔口寬度不變，堵塞拱的特征參數(shù) μθi與i呈線性關(guān)系時，堵塞拱結(jié)構(gòu)為凸性結(jié)構(gòu)；無論是改變筒倉寬度還是孔口寬度，當(dāng)顆粒流過筒倉孔口的流量達到8.3 g/s以上，且使用凸性拱的特征參數(shù)計算隨機游走模型時，理論堵塞概率可以更加準確地預(yù)測出筒倉的實際堵塞概率。