王涵斌 史鼎元 王文玲

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

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

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

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

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

1 實(shí)驗(yàn)設(shè)計(jì)

1.1 休止角生成的方法

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

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

1.2 實(shí)驗(yàn)儀器的搭建

傳統(tǒng)的休止角生成和測(cè)量如圖2所示[2、5]。由圖2可以看到，此裝置(圖2(a)[2])有一些不足之處：第一，支架和漏斗的固定處可能存在轉(zhuǎn)動(dòng)自由度，使漏斗的出口方向和下面的承接臺(tái)不完全垂直，造成誤差。第二，漏斗的固定臂會(huì)因?yàn)槁┒分亓Φ淖饔卯a(chǎn)生撓度與支端轉(zhuǎn)角。第三，承接臺(tái)不平也會(huì)影響到實(shí)驗(yàn)準(zhǔn)確度。

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

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

圖3 實(shí)驗(yàn)儀器

1.3 測(cè)量方法設(shè)計(jì)

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

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

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

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

圖4 測(cè)量流程

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

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

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

圖5 實(shí)驗(yàn)數(shù)據(jù)

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

2 理論分析

2.1 液橋力的計(jì)算

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

(1)

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

圖6 顆粒液橋模型

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

(2)

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

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

計(jì)算式如下：

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

(3)

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

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

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

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

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

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

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

θ=0.0000793381k+35.356

(4)

為了確定修正系數(shù)的具體函數(shù)表達(dá)式，需要使用加入修正系數(shù)的公式和上述關(guān)系，使用實(shí)驗(yà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)計(jì)算顆粒液橋體積，再將顆粒液橋體積代入式(3)計(jì)算系數(shù)F*，隨后將系數(shù)F*與式(5)中的系數(shù)代入式(6)來計(jì)算黏聚能量密度，最后將黏聚能量密度代入式(4)即可得到休止角。

(6)

圖7 引入δ的擬合曲線

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

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

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

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

2.2.1 線性階段

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

2.2.2 水平階段

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

2.2.3 躍升階段

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

2.2.4 下降階段

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

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

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

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

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

2.3.1 線性階段

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

2.3.2 水平階段

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

2.3.3 躍升階段

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

2.3.4 下降階段

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

3 結(jié)語

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