沙堆休止角與含水量關系探究

王涵斌 史鼎元 王文玲

(北京航空航天大學 1航空科學與工程學院； 2高等理工學院；3物理科學與核能工程學院，北京 100083)

顆粒休止角是在重力場中，顆粒在粉體堆積層的自由斜面上滑動時所受重力和粒子之間摩擦力達到平衡而處于靜止狀態下測得的最大角。它對泥石流預測與土壤結構穩定等問題的探索具有重要意義，日本學者Yoshioka N. 在研究地震中的土壤穩定性問題時就使用休止角作為土壤穩定性的力學指標[1]。

研究表明，影響休止角大小最為明顯的因素是含水量[2]。自然世界中土壤顆粒的半徑、形狀、密度等性質相對穩定(不會受環境影響而急速改變)，含水量則不同，它會受降水的影響在短時間顯著變化，極大地改變土壤的力學性質。故含水量對土壤結構的影響既是很主要的，也具有很強的實際意義。因此本文將重點研究沙堆休止角隨含水量的變化關系。但目前對于含水量與休止角的關系研究并不成熟。在《顆粒物質力學導論》[3]一書中，建立了含水量與顆粒間液橋力的大小關系，但是沒有將含水量與休止角兩個宏觀量建立聯系。 Cundall P.A.等人發表文章提出了離散元方法(discrete element method, DEM)，這種方法的基本思想是把每個顆粒看作獨立的微元，通過給定的相互作用力和邊界條件，使用計算機蠻力迭代求解[4],這可以建立顆粒粘聚能量密度與休止角的關系，依舊沒有將兩個宏觀量建立起聯系。

在實驗方面，本文搭建了合理限制自由度的休止角生成裝置，并使用計算機視覺的處理方法，利用OpenCV這款開源的計算機視覺庫，把圖像的顏色轉化成灰度值，將灰度值平均運算，進行平滑處理，剔除圖像中不必要的信息。最后計算灰度變化梯度，因為沙堆的顏色和背景差異很大，梯度較大的地方可以推測為沙堆的輪廓。有了輪廓上每個點的坐標，就可以使用數學工具進行直線擬合，進而通過直線斜率得到角度大小。

我們的視覺處理方法可以準確地進行了休止角測量，由此可以得到可靠的含水量與休止角數據。在理論方面，初步構建了休止角-含水量數學模型，并給出了合理的物理解釋。

1 實驗設計

1.1 休止角生成的方法

常見的休止角生成方法有3種：注入法、排出法和傾斜法。注入法是在漏斗中添加粉體使之下落在承接臺上形成錐體的方法，休止角為錐體側面與地面夾角。排出法是使粉體從容器流出部分之后剩余粉體形成的凹陷的方法，休止角是凹陷面與水平面的夾角。傾斜法是在繞水平軸慢速回轉的圓筒容器內加入1/2～1/3容積的粉體，當粉體的表面產生滑動時，測定其表面的傾角的方法。具體模式如圖1所示。

3種方法中，排出法形成的休止角位于裝置內部，不便于測量；傾斜法是動態過程實驗裝置較為復雜，休止角的生成角度誤差較大。特別是不便于使用本實驗的測量方法。因此，我們選擇注入法。

1.2 實驗儀器的搭建

傳統的休止角生成和測量如圖2所示[2、5]。由圖2可以看到，此裝置(圖2(a)[2])有一些不足之處：第一，支架和漏斗的固定處可能存在轉動自由度，使漏斗的出口方向和下面的承接臺不完全垂直，造成誤差。第二，漏斗的固定臂會因為漏斗重力的作用產生撓度與支端轉角。第三，承接臺不平也會影響到實驗準確度。

圖2 傳統休止角測定儀模式圖

針對上述3個問題，本文分別加以改進。

圖3 實驗儀器

1.3 測量方法設計

沙堆是一個三維錐體，不便于測量，因此可以通過拍照將其轉化為二維。但是沙堆并不是標準錐體，圖像邊緣不是直線，人手工測量會有一定的主觀性和不確定性，因此采用計算機視覺技術進行測量。OpenCV是一款開源成熟度的計算機視覺庫，應用這個庫，可以在照片上找出沙堆的輪廓，并把輪廓上點的坐標進行直線擬合，直線的斜率即對應休止角的正切值，這樣就規避了人在測量非標準直線構成的夾角時的主觀性。

在實際拍攝的時候，使用三腳架固定了鏡頭的位置，保持鏡頭和水平面垂直，且距離被測沙堆20cm左右。這樣每次拍攝的位置固定，可以避免因為鏡頭傾斜或俯仰造成的誤差。另外通過保持被測物體在照片中間等手段大大減少了手機鏡頭引起的鏡頭畸變或感光元件間隔不均勻造成的測量誤差。

顆粒在下落過程中，使用漏斗上方固定的攪拌裝置緩慢攪拌，這讓顆粒流出漏斗后具有微小且隨機的水平速度，保證了沙堆的形狀是近似中心對稱的，加上多次測量求平均值的處理，就可以消除因為拍攝位置固定而可能帶來的系統誤差。

為了在照片上找到輪廓，要在灰化的基礎上采用Sobel算子進行梯度檢測，這樣有助于定位邊界的位置和檢測邊緣，之后再使用模糊化處理通過灰度梯度的變化找到邊界，就可以計算休止角的大小了。其完整的測量流程和各階段效果如圖4。

圖4 測量流程

該程序能夠十分準確地找到沙堆的邊界，并且通過直線擬合推測休止角的大小，完全規避了人在測量中的主觀性，測量十分準確。

另外，由于拍照通常是在手機上進行，但是程序的運行和調試使用電腦更加方便，為了彌合兩者之間的矛盾，搭建了一個小服務器實現前后端架構。手機瀏覽器作為前端上傳圖片并顯示結果，服務器作為后端進行測量。

此測量方式跟傳統測量方式(圖2(b)[5])比較，有了很大改進。首先，它在測量時保持了沙堆的完整性，從而能夠保證數據的可靠性，避免了傳統測量方式只生成了半個沙堆而一定程度上破壞了沙堆的力學性質的缺憾。其次，傳統測量方式進行測量時需要使用標尺，這個操作極易對沙堆造成破壞導致測量不準確。而拍照的方式能夠最大程度地對被測對象進行保護。此外，在得到沙堆輪廓之后采用數學方法擬合直線計算角度，也要好于人的手工測量。因為沙堆的輪廓往往不是嚴格的直線，這使得手工測量極易因有主觀性造成誤差，而采用計算機數值計算的方法，可以消除這種主觀性造成的影響。

圖5 實驗數據

通過此實驗裝置，我們實驗得到了一系列沙堆含水量和休止角的關系，列在圖5中，橫坐標是含水量(水的質量與顆粒質量之比)，變化范圍從0到0.12，每變化0.01進行一組實驗，縱坐標是休止角(大小用角度表示)。為得到充分的實驗數據，同一含水量重復實驗3次，由圖顯示，同一含水量下3次實驗數據很相近，進一步驗證了實驗結果的可靠性。

2 理論分析

2.1 液橋力的計算

我們利用液橋力的相關理論，進一步理論分析沙堆休止角與含水量的關系。計算液橋力的前提是確立液橋模型(圖6)，并根據宏觀含水量確定液橋含水體積。為簡化模型，假設顆粒為直徑相同的球體，根據含水量計算液橋體積如下：

(1)

式(1)中，ρp是顆粒密度；Rp是顆粒半徑；nc°是配位數；ρ1是液體密度；ω0是含水量。

圖6 顆粒液橋模型

液橋力是表面張力與靜水壓力之和

(2)

式(2)中，γ是填隙液體的表面張力；Δp液橋汽液界面壓強差；rN是液橋頸部曲率半徑。

為簡化式(2)的計算，本文采取其回歸形式.

計算式如下：

lnF*=f1°°-f2°exp(f3°lnS++f4°ln2S+)

(3)

2.2 修正系數的引入

《顆粒物質力學導論》[3]中使用液橋理論和水的表面張力理論，認為液體在顆粒間以小液柱(液橋)的形式存在，建立了含水量與顆粒間液橋力的大小關系，但是其公式中液橋與顆粒之間的接觸角、液橋的距離都不能確定。因此公式使用起來有困難。

微小尺度上成千上萬的顆粒、液橋相互作用最終決定宏觀尺度上休止角大小，其物理過程需要的巨量計算在計算機技術發展后得以解決。Cundall P.A.等人發表文章提出了離散元方法(discrete element method, DEM)，通過給定的相互作用力和邊界條件，使用計算機蠻力迭代求解[4],這可以建立顆粒粘聚能量密度與休止角的關系，但是粘聚能量密度與液橋力之間的關系取決于目前無法確定的液橋與顆粒的接觸面積。

在現有理論中，存在許多未知量，比如我們很難確定液橋距離S(參考文獻[3]給出了液橋距離的計算公式，但實際距離與極限距離的關系難以確定)，顆粒與液橋的接觸角，以及在從液橋力向顆粒粘聚能量密度過渡時我們難以確定軟球模型中顆粒的接觸面積A。

諸多不確定因素使該公式無法計算出準確解，本文將不確定因素造成的影響集中起來引入修正系數修正。這些不確定因素都和含水量有關，因此修正系數應為含水量的函數設為δ(ω)。

EDEM是大型商用離散元仿真軟件，它的原理是從某些在外力(如重力)條件下發生微小位移的離散元(如顆粒)開始計算，通過位移引發的受力變化計算其他顆粒發生的位移，不斷迭代計算，直至所有的顆粒達到穩定狀態，不再發生位移，就得到了受力平衡時的沙堆形狀。

仿真時，需要控制的參數主要有4個：顆粒的形狀、顆粒的半徑、顆粒的密度以及粘聚能量密度。顆粒的形狀、半徑和密度應該與實驗中實際使用的顆粒相近，所以顆粒形狀近似為球形，平均半徑取0.5mm，顆粒密度使用二氧化硅的密度2600kg/m3。而粘聚能量密度則是仿真的自變量，在其他3個參數確定不變的情況下，每指定一個粘聚能量密度，就進行一次仿真得到一個對應的休止角。仿真時的粘聚能量密度的范圍是50000～140000J/m3，每隔10000J/m3做一次，一共得到了10組粘聚能量密度k和對應休止角θ的大小，發現有如式(4)的線性關系：

θ=0.0000793381k+35.356

(4)

為了確定修正系數的具體函數表達式，需要使用加入修正系數的公式和上述關系，使用實驗數據進行最小二乘法擬合。發現修正系數的表達式采用如下形式效果最好。

(5)

式中參數如下：

p1=-51.01,p2=14.69,p3=-1.241,

p4=0.03918,p5=-0.0003692,

p6=6.177×10-7,q1=-0.09502,

q2=0.002622,q3=-6.745×10-5,

q4=3.11×10-6。

進一步，在理論模型下的擬合曲線如圖7所示。至此我們找到了休止角與含水量的單值關系式。首先將含水量ω代入式(1)計算顆粒液橋體積，再將顆粒液橋體積代入式(3)計算系數F*，隨后將系數F*與式(5)中的系數代入式(6)來計算黏聚能量密度，最后將黏聚能量密度代入式(4)即可得到休止角。

(6)

圖7 引入δ的擬合曲線

2.2 休止角階段性變化的性質描述

圖8 休止角隨含水量變化的四個階段

由圖8的擬合曲線，可以看到，當沙堆含水量在逐漸增大的時候，其休止角一直處于增大階段，這說明沙堆在含水量較小時，隨著含水量的增加，休止角依舊可以穩定增加，但當含水量達到一定區域約8%時，沙堆休止角就將隨著其含水量迅速下降，此時沙堆的穩定性越來越差，直至坍塌。

進一步，我們分別對這4個階段依次分析。

2.2.1 線性階段

含水量在0%～4%左右，休止角-含水量曲線，近似為一條直線，休止角與含水量成正比，在該階段過程中隨著含水量的增加休止角也在增加，說明土壤的蓄水能力，沒有達到飽和，還可以繼續儲水，在降水量較少時，土壤穩定性高。

2.2.2 水平階段

含水量在4%～5%左右，休止角與含水量之間不再保持正比關系，當含水量增加的時候，休止角-含水量曲線出現水平線段，休止角不再隨含水量的增加而增加。此時為休止角的穩定狀態階段，含水量微量變化時，休止角可以維持在一個固定的角度。

2.2.3 躍升階段

含水量在5%～8%左右，先是含水量的小范圍增加導致的休止角突變，之后隨著含水量的增加休止角繼續增加，且增長速度較快，該階段沙堆依然保持一定的儲水能力，但相比于近似線性階段，儲水能力下降，休止角隨含水量的變化的反饋程度變得劇烈，8%為該階段的極點，為非穩定的狀態。土壤變得粘稠因而可以更穩定。

2.2.4 下降階段

含水量大于8%，此時土壤的儲水能力將遭到破壞，休止角隨含水量的增加而減少，系統穩定性喪失，這表明降水量很大的時候土壤會從粘稠變滑，穩定性下降。

2.3 休止角階段性變化的理論解釋

由圖9，可以看到修正系數大小變化的特征為，在不同的含水量區域修正系數隨含水量的變化趨勢不同，可以把含水量的區間分為0～4%，4%～5%，5%～9%，9%～12% 4個區間。

圖9 隨含水量遞增導致的液橋變化

參考液橋變化圖9[3]，我們進一步給出理論解釋。

2.3.1 線性階段

含水量在0%～4%階段，休止角與含水量成正比，該階段顆粒間液橋通過接觸點上的液橋發生黏連，該階段含水量增加過程，導致休止角增加的原因可以理解為液體逐漸填隙在顆粒之間導致顆粒間的接觸點數目增多，所以該階段是線性階段。

2.3.2 水平階段

含水量在4%～5%階段，接觸點數目已經達到飽和，顆粒間液橋力的成型液橋形狀基本不變，依舊可以等效為鐘擺狀，該階段顆粒液橋達到該形狀的最大容納量。

2.3.3 躍升階段

含水量在5%～8%階段，此時液橋的形狀正在由鐘擺狀轉型為環索狀和毛細管狀，部分顆粒間隙充滿液體黏連力增加。導致了休止角的突變和后續休止角的增加。

2.3.4 下降階段

含水量大于8%，此時顆粒間有填隙液體，液體壓力開始增加，液橋發生破壞，所以呈現微小的下降階段。

3 結語

本文搭建了合理限制自由度的休止角生成裝置，提出了以計算機視覺方式進行休止角測量的方法。進一步，通過引入有效修正系數理論，構建了休止角-含水量數學模型，并給出了與實驗相吻合的物理解釋。這有利于人們在已知含水量的情況下對休止角進行合理預測，對于泥石流預測等土力學應用問題，將具有重大意義。