沙堆休止角與含水量關(guān)系探究

王涵斌 史鼎元 王文玲

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

顆粒休止角是在重力場中，顆粒在粉體堆積層的自由斜面上滑動時所受重力和粒子之間摩擦力達(dá)到平衡而處于靜止?fàn)顟B(tài)下測得的最大角。它對泥石流預(yù)測與土壤結(jié)構(gòu)穩(wěn)定等問題的探索具有重要意義，日本學(xué)者Yoshioka N. 在研究地震中的土壤穩(wěn)定性問題時就使用休止角作為土壤穩(wěn)定性的力學(xué)指標(biāo)[1]。

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

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

我們的視覺處理方法可以準(zhǔn)確地進(jìn)行了休止角測量，由此可以得到可靠的含水量與休止角數(shù)據(jù)。在理論方面，初步構(gòu)建了休止角-含水量數(shù)學(xué)模型，并給出了合理的物理解釋。

1 實驗設(shè)計

1.1 休止角生成的方法

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

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

1.2 實驗儀器的搭建

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

圖2 傳統(tǒng)休止角測定儀模式圖

針對上述3個問題，本文分別加以改進(jìn)。

圖3 實驗儀器

1.3 測量方法設(shè)計

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

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

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

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

圖4 測量流程

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

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

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

圖5 實驗數(shù)據(jù)

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

2 理論分析

2.1 液橋力的計算

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

(1)

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

圖6 顆粒液橋模型

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

(2)

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

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

計算式如下：

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

(3)

2.2 修正系數(shù)的引入

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

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

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

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

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

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

θ=0.0000793381k+35.356

(4)

為了確定修正系數(shù)的具體函數(shù)表達(dá)式，需要使用加入修正系數(shù)的公式和上述關(guān)系，使用實驗數(shù)據(jù)進(jìn)行最小二乘法擬合。發(fā)現(xiàn)修正系數(shù)的表達(dá)式采用如下形式效果最好。

(5)

式中參數(shù)如下：

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。

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

(6)

圖7 引入δ的擬合曲線

2.2 休止角階段性變化的性質(zhì)描述

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

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

進(jìn)一步，我們分別對這4個階段依次分析。

2.2.1 線性階段

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

2.2.2 水平階段

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

2.2.3 躍升階段

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

2.2.4 下降階段

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

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

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

圖9 隨含水量遞增導(dǎo)致的液橋變化

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

2.3.1 線性階段

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

2.3.2 水平階段

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

2.3.3 躍升階段

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

2.3.4 下降階段

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

3 結(jié)語

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