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

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

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

引 言

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

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

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

1 實 驗

1.1 實驗裝置

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

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

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

1.2 堵塞拱的特征提取

Fig. 2 Schematic diagram of feature extraction of jamming arch

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

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

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

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

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

1.3 堵塞概率計算

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

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

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

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

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

2 結果與討論

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

2.1 堵塞概率的統計

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

2.2 堵塞拱的特征值

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

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

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

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

2.3 流量對堵塞概率的影響

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

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

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

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

3 結 論

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