基于離散元方法的松散體滑動堆積特性及影響因素分析*

成浩 韓培鋒2)? 蘇有文

1) (西南科技大學土木工程與建筑學院, 綿陽 621010)

2) (水利部山洪地質災害防治工程技術研究中心, 武漢 430010)

1 引 言

松散體是由大量的松散顆粒組成, 在自然界廣泛存在, 具有結構松散、孔隙度大、粒間結合力差和易變形等特點, 還衍生著諸多的動力學行為和物理力學特性; 因結構構成的特殊性, 往往在自重或其他外界因素的影響下失穩破壞, 造成河流堵塞、道路阻斷及房屋掩埋等危害[1?5]. 當前, 對于松散體滑動后的堆積范圍很難科學預測及易損性評估,已成為國際上的研究熱點.

學者們在研究松散體滑動后的堆積特性時, 主要考慮初始松散體的體積、顆粒粒徑、坡高、啟動區坡度、坡角、摩擦系數和地震荷載等因素對松散顆粒變形過程和滑動后的堆積形態、沖程、平均粒徑分布、固相分數及分選系數等指標的影響規律[6?10].堆積形態與顆粒運動過程及速度密切相關, 通過物理模型試驗和理論分析提出了滑坡沖程與速度的預測方法[11]、落石運動速度與距離的數學公式[12]、松散體運動中的障礙物與沖程的動態關系模型[13].可見, 模型試驗和理論分析為其他方法研究松散顆粒滑動堆積特性提供了前提, 但缺乏滑動過程的實測數據去解釋松散顆粒滑動堆積特性的內在機理.此外, 在現場調研中發現坡度對松散體滑動堆積范圍分布各異, 松散體中的含石量大小與滑后堆積特性之間的內在關系仍需進一步定量研究.

離散元方法(discrete element method, DEM)最早由Cundall[14?16]教授1971 年提出, 當前正廣泛用于研究非連續介質顆粒運動等問題[6,17]. 由于現實的顆粒運動體量太大, 這會大幅度降低DEM 的計算效率; 學者們通常借鑒物理模型試驗中的縮尺滑槽試驗, 并采用DEM 研究松散顆粒的變形失穩過程及滑動后的顆粒分布特性[6]、坡高和坡角對滑后堰塞體的內部結構和堆積性質的影響規律[8]. 顆粒間和顆粒與滑動邊界間的摩擦系數對顆粒堆積特性的影響顯著, 摩擦系數對顆粒堆積角的影響較大[18,19]. 前述離散元研究都有一個共同點, 均是基于球形顆粒開展研究, 而自然界中的球形物質實際上很少; 近年來學者們開始研究非球形顆粒, 這就使得研究更加貼近現實. 滾動摩擦系數對非球形顆粒堆積角和顆粒堆內外邊界間隙z的影響較大[17].

松散體滑動后堆積特性的內在機理與運動中的顆粒能量及動力學關系密切. 根據牛頓第二定律及顆粒相互作用時的接觸力關系, 可從顆粒流分析法和能量交換方面分析顆粒整個運動過程[20,21]. 顆粒間的動勢能轉化與顆粒的摩擦聯系緊密, 可用土力學方法分析顆粒在運動中的能量衰減[22,23]. 從上述研究中可以看出, 顆粒間的能量及接觸力為解釋松散顆粒滑動堆積特性提供了切入點.

因此, 本文采用DEM 構建了松散體滑動堆積三維模型, 對非球形顆粒滑動堆積過程進行了數值模擬, 研究了含石量和坡度變化對松散體滑動堆積特性的影響規律, 分析了顆粒滑動堆積過程中的能量交換和接觸作用力, 為深入理解和進一步揭示松散顆粒堆積特性的內在機理提供參考.

2 顆粒滑動堆積的離散元模型

2.1 DEM 的計算原理[14?16,24]

DEM 認為顆粒單元間的運動是相互獨立的,只有顆粒單元間發生相互作用時, 才會在接觸點處產生作用力. 根據力與位移關系, 由位移得到顆粒受到的作用力, 再由牛頓第二定律得出顆粒i的運動方程

利用中心差分法對(1)式進行積分, 得到以兩次迭代時間步長的中間點表示的更新速度表達式為

式中, ?t是時間步長;N對應時間t.

對(2)式進行積分, 可得到關于位移的等式表述為

進而得到顆粒單元新的位移值, 并將該新位移代入力與位移關系計算新的作用力, 循環反復就可實時跟蹤每個顆粒單元在任意時刻的運動, 如圖1所示.

圖1 顆粒計算迭代示意圖 (a) 力與位移關系; (b) 理論計算Fig. 1. Granular computing iteration diagram: (a) Relationship between force and displacement; (b) theoretical computing.

2.2 Hertz-Mindlin(no slip)接觸模型[14?16,25]

當顆粒與顆粒發生相互作用時, 會產生相應的接觸法向力Fn和切向力Ft, 如圖2 所示.

圖3 為顆粒法向力示意圖,Fn表達式為

其中,E*為等效彈性模量, 由顆粒的彈性模量E和泊松比v求得;R*為等效顆粒半徑, 由顆粒的半徑R求得, 其表達式為

顆粒間發生作用的法向重疊量a由顆粒半徑R和顆粒的球心位置矢量r求得, 即

圖2 顆粒相互作用時的受力Fig. 2. The forces by particles interacting.

圖3 顆粒法向力示意圖 (a) 法向重疊量; (b) 位置關系Fig. 3. Diagram of normal force of granular: (a) Normal overlap; (b) position relations.

式中,a為顆粒間接觸面為圓形時的接觸半徑.

其中,b為法向阻尼力系數, 由顆粒間相互作用時的彈性恢復系數e求出;Sn為法向剛度, 由顆粒的等效模量、等效顆粒半徑和顆粒的法向重疊量求出;m*為等效質量, 由顆粒的質量m求出; 即

在顆粒接觸的切線方向上, 其切向接觸力Ft, 切向阻尼力的表達式為

式中,G*為等效切向模量; 切向力Ft由切向重疊量d和切向剛度St確定;為切向相對速度.

2.3 數值模擬計算方案

為探討松散體滑動后的堆積特性, 以堆積體中的含石量s(%)和滑床坡度q(°)為變量參數開展數值模擬. 含石量s的取值分別為0, 30, 50, 70;坡度q的取值分別為30, 45, 65. 本研究共開展12 組模擬, 每組模擬三次取平均值作為該組的模擬結果.

本次數值模擬的主要目的是針對松散體滑動后堆積特性的規律及機理研究. 因此, 借鑒碎屑流滑槽物理模型試驗[26,27], 基于相似比原理開展數值模擬, 底板模擬堆積區, 初始堆積體由黃色粗顆粒和藍色細顆粒在模擬軟件EDEM 2.7[17]中自動生成并自然堆積, 如圖4 所示. 需要指出的是, 本文數值模型的邊界條件在處理時參照了文獻[6]和文獻[8]中的研究思路, 將料箱和底板沿滑床中軸線設置, 滑床和底板均為全約束并保持靜止狀態; 在松散體顆粒整個滑動堆積過程中, 滑床和底板的彈性形變忽略不計, 只考慮顆粒與滑床間以及顆粒與底板間的相互摩擦和碰撞作用. 此外, 顆粒流為固相運動, 僅在重力作用下滑動, 忽略外力(如風力等)的干擾.

圖4 松散體滑動堆積模型示意圖 (a) 三維數值模型;(b) 側視圖; (c) 俯視圖Fig. 4. Diagram of sliding accumulation model of loose accumulation body: (a) Three-dimensional numerical model;(b) side of view; (c) vertical view.

依據現場典型地質災害點的松散堆積體現場取樣, 通過室內的土工試驗, 對松散體中的塊石與碎石土粒徑及其體積范圍進行統計, 并參照文獻[28]中對松散體中顆粒粒徑及其分形理論的研究結論,最后結合模擬實際計算效率綜合確定粗顆粒體積約為細顆粒體積的46.9 倍, 滿足現場調查及相關研究的范圍要求. 其中, 模擬的大粒徑塊石和碎石土單個顆粒均由4 個球形基礎顆粒組成, 如圖5 所示. 通過野外調研和土工試驗, 確定了模擬材料的本征參數和彈性恢復系數, 顆粒與滑槽間的摩擦系數通過物理試驗確定[29], 主要計算參數見表1.

圖5 單個顆粒大樣圖 (a) –Z 視角; (b) +Y 視角Fig. 5. Diagram large sample of single particle: (a) –Z of view; (b) +Y of view.

2.4 顆粒滑動堆積過程的試驗對比驗證

為了驗證本文數值計算的準確性與穩定性, 首先對參數為坡度為38°, 體積為0.015 m3, 顆粒粒徑為5—20 mm, 單個顆粒形狀如圖5 所示, 其他參數見表1 的松散體滑動過程進行了數值模擬計算, 并與文獻[30]中的物理模型試驗結果進行對比, 如圖6 所示. 可知, DEM 計算結果累積體積開始的時間略早, 最后達到穩定階段的累積體積較試驗結果稍小, 可能是因為試驗中滑槽的摩擦系數偏小, 但總體變化的規律與試驗結果吻合良好, 說明顆粒形狀對累積體積的影響較小, 同時反映了本文DEM 數值模型可以較好地模擬松散體的滑動過程.

表1 松散顆粒堆積離散元模擬的主要計算參數Table 1. Main computational parameters of discrete element simulation for loose granular accumulation.

圖6 顆粒滑動后累積體積試驗值[30]與DEM 計算值比較Fig. 6. Comparison of cumulative volume between experiment results[30] and DEM simulation results of granular after sliding.

3 數值模擬結果分析

3.1 堆積過程與速度分布

圖7 為坡度為65°、含石量為50% 的松散體滑動堆積過程與速度分布.t= 400 ms 時, 初始堆積體在重力下呈傾覆式變形, 前端和表層顆粒的啟動速度較內部和底層大.t= 600 ms 時, 因顆粒流頭部受約束較小, 所以顆粒速度較大, 少量顆粒脫離顆粒流運動.t= 860 ms 時, 顆粒流長度延伸最大,而厚度越扁平, 平均速度達到最大, 最早接觸底板時易發生彈跳.t= 1000 ms 時, 由于底板限制了顆粒流的運動方向, 部分顆粒到達坡底便開始減速堆積, 顆粒流長度變短的同時其厚度短暫增大; 速度出現明顯分層, 底層最小, 中間層次之, 頂層最大.t= 1 200 ms 時, 顆粒流長度再次延展, 相反的厚度逐漸減小; 頭部速度最大, 中部次之, 尾部最小. 隨著時間的推移, 顆粒漸漸準靜態堆積, 呈藍色狀態, 堆積結束.

綜上, 松散體滑動堆積大體可劃分為啟程變形加速, 近程高速滑動, 中程減速滑動, 遠程準靜態堆積四個階段; 這里只是定性的分析松散體堆積過程中長度、厚度、形狀及速度的變化, 接下來將對堆積結果進行定量分析.

3.2 含石量對堆積特性的影響

本文主要分析含石量和坡度變化對松散體滑動后堆積特性的影響規律, 最終得到12 組不同的顆粒堆. 通過提取顆粒堆中與底板接觸的顆粒位置坐標, 導入Origin 2018 軟件中, 采用網格劃分法,選取特殊交點可得到平面堆積區域的簡化計算示意圖, 進一步獲取堆積特征值(沖程l, 堆積寬度d,最大厚度h, 平面堆積面積S), 如圖8 所示. 其中,沖程l是指松散體滑動后沿坡腳線至主堆積區最前緣的距離, 是地質災害防治和預測的常用指標;堆積寬度d是指在堆積區中與沖程方向垂直的最大距離; 最大厚度h為堆積區豎直方向上底層至表層的最大距離; 平面堆積面積S則是松散體地質災害范圍預測、評估和保護區搬遷范圍的重要參考指標.

圖9 給出的是含石量對最終堆積特征值的影響.分析可知, 隨著含石量增大, 顆粒堆的沖程、堆積寬度、堆積面積均先增大后遞減, 而最大厚度呈先減小后遞增趨勢. 主要是因為運動中細顆粒的減少,使得粗顆粒間的空隙及粗顆粒與滑槽間的摩擦都相對變大, 粗顆粒相互碰撞作用增強, 能量耗散較大.

靜堆積角g是指底平面保持靜止時松散顆粒在重力下堆積后的自然坡度表面與水平面間的夾角, 它是研究散體顆粒材料基本物理堆積特性的常用測量參數[25]. 為了定量分析不同計算條件下靜堆積角的變化規律, 同時為了靜堆積角的測量精確, 對顆粒堆X軸正方向和X軸負方向采用Origin 2018 圖像數字化技術獲得顆粒堆的邊界輪廓線, 并對邊界輪廓線做線性擬合, 進一步獲取擬合方程及其斜率k. 如圖10 所示. 最終測定靜堆積角的表達式為

圖7 松散體滑動堆積全過程(+X 視角) (a) t = 400 ms; (b) t = 600 ms; (c) t = 860 ms; (d) t = 1 000 ms; (e) t = 1 200 ms;(f) t = 2 000 ms; (g) 堆積過程形狀變化Fig. 7. Loose materials the whole process of sliding accumulation(+X view): (a) t = 400 ms; (b) t = 600 ms; (c) t = 860 ms; (d) t =1 000 ms; (e) t = 1 200 ms; (f) t = 2 000 ms; (g) accumulation process changes shape.

圖8 最終平面堆積形態示意圖Fig. 8. The final diagram of plane accumulation form.

由表2 分析可知, 隨著含石量增大, 在一定坡度時, 靜堆積角會出現小幅度減小, 而超過一定坡度時, 靜堆積角會小幅度增大. 由于靜堆積角與松散顆粒的流動性密切相關, 從而反映了松散體中含石量值的大小對顆粒流動性有不同程度的影響, 坡度較小時顆粒流動性有增強的趨勢, 而坡度較大時顆粒流動性則有所減弱; 主要是因為含石量變化會使得顆粒滑動過程中的碰撞強度、摩擦及接觸形式在一定程度上改變顆粒的流動性, 除此之外影響顆粒流動性的因素還有很多, 如黏聚力、內摩擦力和初始堆積體密度等, 所以未來的工作會針對這一點進一步研究. 另外, 相比于含石量, 坡度對靜堆積角值的影響顯著; 換句話說, 坡度對顆粒流動性的影響關系更加密切, 至于坡度是如何影響靜堆積角和顆粒流動性的, 以及二者的關系模型在后文會進一步分析.

為了進一步建立含石量與松散體滑動后靜堆積角的關系模型, 同時考慮到坡度在一定程度上對含石量的影響效應, 在其他參數不變的條件下, 對表2 中的靜堆積角g均值與含石量s值進行關系式擬合, 得到擬合關系模型方程式為

在建立的靜堆積角g均值與含石量s值的關系模型((20)式)可以看出, 二者之間的關系近似二次拋物線, 相關性系數R2趨近于1, 說明靜堆積角與含石量的相關性擬合很好; 也即是說, 在一定條件下, 這種關系模型可以從散體顆粒流動性的基本物理特性角度出發為該類松散體滑動后的靜堆積角對含石量的響應提供預測模型.

圖9 含石量對堆積形態的影響 (a) 沖程; (b) 堆積寬度; (c) 最大厚度; (d) 堆積面積Fig. 9. Influence of stone content on the accumulation form: (a) Stroke; (b) accumulation width; (c) maximum thickness; (d) accumulation area.

圖10 靜堆積角邊界輪廓提取Fig. 10. Static accumulation angle boundary contour acquire.

表2 不同計算條件下靜堆積角測量值Table 2. Measured value of static accumulation angle under different computing conditions.

為研究含石量對堆積區顆粒稀疏程度的影響,這里以坡度65°、含石量50%的松散體滑動堆積區為例. 在模擬軟件中, 以坡腳線為基線沿顆粒主流運動方向每隔150 mm 劃分一個區域, 并自動獲取每個區域的顆粒體積, 進而分析A1～A10 的顆粒稀疏程度, 體積計算平面示意如圖11 所示.

由圖12 分析可知, 距離坡腳線越遠, 顆粒堆體積逐漸減小, 顆粒越稀疏. 隨著含石量增大, 堆積區體積有所下降, 但差異不明顯; 其中含石量70%的體積最小, 下降速度最快. 主要原因是粗顆粒間的相互碰撞及粗顆粒與滑面的摩擦變化.

圖11 堆積區體積計算平面示意圖Fig. 11. Schematic diagram of volume calculation of accumulation area.

圖12 坡度65°時不同含石量對堆積體積的影響Fig. 12. Influence of different stone contents on the accumulation volume at a slope of 65°.

3.3 坡度對堆積特性的影響

圖13 給出的是坡度對最終堆積特征值的影響. 由圖13 分析可知, 雖然坡長相同代表滑槽的滑動距離相同, 但顆粒的勢能隨坡度的增大而增大, 重力勢能轉化為更大的動能. 這就使得相應的沖程, 堆積寬度, 堆積面積都有不同程度的增大,而最大厚度呈減小趨勢; 這一結論與文獻[7]中的相關結論基本相符, 同時也說明了顆粒形狀相較于坡度對堆積特性的影響較小. 其中, 沖程和堆積面積呈線性增大; 堆積寬度隨坡度增大, 增長趨勢變緩. 坡度影響堆積特性的主要原因是: 坡度越小時,坡高越小, 顆粒作用在滑面上的自重力及受到滑面的摩阻力相對增大, 導致顆粒能量損失增大; 此外,坡度越小的顆粒運動速度和沖擊能力較小, 即是說后緣顆粒無法推動或躍過前緣顆粒向前運動.

在以上關于含石量變化對靜堆積角以及顆粒流動性的分析時, 發現坡度變化對靜堆積角的影響顯著, 且隨著坡度增大, 靜堆積角呈遞減趨勢; 這就反映了坡度越大的散體顆粒流動性越好. 為了深入研究坡度因素與靜堆積角的相關關系, 在其他參數不變的情況下, 對表2 中的靜堆積角g均值與坡度q進行多項式擬合, 得到如下關系式

從建立的靜堆積角g均值與坡度q值的關系模型((21)式)可以看出, 采用二次函數擬合時, 得到的相關性系數R2均為1, 說明靜堆積角與坡度的相關性關系程度極強, 幾乎不受含石量變化的影響; 同樣, 在一定條件下, 這種關系模型可以從散體顆粒流動性的基本物理特性角度出發為該類松散體滑動后的靜堆積角對坡度的響應提供預測模型.

圖14 表明, 隨著坡度增大, 堆積輪廓平面形狀由近似圓狀趨于近似橢圓狀, 主要是因為顆粒堆沖程的延伸性較堆積寬度顯著. 同時也說明坡度比含石量對堆積輪廓平面形狀影響明顯.

由圖15 可以觀察到, 無論何種坡度下, 堆積區顆粒的稀疏程度由后緣往前緣逐漸增大; 隨著坡度增大, 顆粒的稀疏程度更加明顯. 其中, 細、粗顆粒在堆積區中的分布各異, 坡度較小時, 粗顆粒大部分均勻分布在細顆粒表層; 隨著坡度增大, 粗顆粒的分布路徑趨于沖積扇, 這一現象與現實中的溝谷碎屑流堆積相似. 說明粗顆粒較細顆粒的運動更加活躍, 平均沖擊與攜帶能力更強.

圖13 坡度對堆積形態的影響 (a) 沖程; (b) 堆積寬度; (c) 最大厚度; (d) 堆積面積Fig. 13. Influence of slope on the accumulation form: (a) Stroke; (b) accumulation width; (c) maximum thickness; (d) accumulation area.

圖14 不同坡度下的堆積輪廓平面形狀 (a) 含石量0%; (b) 含石量30%; (c) 含石量50%; (d) 含石量70%Fig. 14. Plane accumulation morphology under different slope: (a) Stone content 0%; (b) stone content 30%; (c) stone content 50%;(d) stone content 70%.

圖16 為含石量50%的松散體在不同坡度下滑動后的體積計算結果. 分析可知, 當坡度為30°和65°時, 距離坡腳線越遠, 顆粒堆體積逐漸減小, 顆粒越稀疏; 坡度45°時, 顆粒堆體積在坡腳線附近形成隆起區, 當距離持續增大, 顆粒堆體積依然逐漸減小.

圖15 含石量50%時不同坡度松散顆粒滑動堆積模擬結果 (a) 30°; (b) 45°; (c) 65°Fig. 15. The results of the sliding accumulation simulation of stone content 50% loose granular with different slopes:(a) 30°; (b) 45°; (c) 65°.

圖16 含石量50%時不同坡度對堆積體積的影響Fig. 16. Influence of different slope of 50% stone content on the accumulation volume.

進一步分析初始松散體中相同比重的粗、細顆粒分別在堆積區體積中的占額, 如圖17 所示. 可以觀察到, 粗、細顆粒的體積占額均出現一個臨界距離Lc, 當距坡腳線距離L < Lc時, 細顆粒大于粗顆粒所占體積; 當L > Lc時, 粗顆粒反而較細顆粒體積大; 坡度越大, Lc值越大, 即距離坡腳線越遠. 說明粗顆粒的平均位移能力較細顆粒大, 換句話說, 粗顆粒較細顆粒平均運動距離更遠. 此外,由于實際的松散體滑動堆積形式和運動過程十分復雜, 與顆粒的能量耗散和接觸力等因素, 以及顆粒材質、粒徑和堆積體密集度等因素聯系緊密; 因此, 不同邊界計算條件下的臨界距離存在差異.

圖17 含石 量50%時不同坡度顆粒 體積 對比 (a) 30°;(b) 45°; (c) 65°Fig. 17. The volume comparison of granular with different slopes with 50% stone content: (a) 30°; (b) 45°; (c) 65°.

3.4 含石量和坡度對堆積區顆粒質量補給的影響

累積質量或累積體積常用于表征松散體滑動后對河流堵塞程度或公路掩埋等的影響[6,30]. 圖18給出的是同一坡度下, 不同含石量的松散體對堆積區顆粒累積質量供應時程曲線. 分析可知, 含石量為30%, 50%和70%時, 松散體幾乎同一時間開始對堆積區進行顆粒質量供給, 累積質量曲線增長速度基本一致; 在達到穩定值時, 含石量越大, 堆積區的顆粒累積質量越小.

圖19 給出的是同一含石量下, 不同坡度的松散體對堆積區顆粒累積質量供應時程曲線. 分析可知, 坡度越小的松散體對堆積區的顆粒質量供給開始時間越晚, 且供給持續時間越長; 隨著坡度增大,累積質量曲線的增長率越大. 在達到穩定值時, 坡度越大, 堆積區的顆粒累積質量越大; 其中, 坡度為30°—45°的顆粒累積質量增幅十分顯著, 當坡度持續增大時, 增幅趨于平緩.

圖18 坡度65°時不同含石量對累積質量的影響Fig. 18. Influence of different stone contents on cumulative mass at slope of 65°.

圖19 含石量50%時不同坡度對累積質量的影響Fig. 19. Influence of different slope on cumulative mass at stone content of 50%.

圖20 含石量50%時各坡度顆粒累積質量對比 (a) 30°;(b) 45°; (c) 65°Fig. 20. The cumulative mass comparison of granular with different slopes with 50% stone content: (a) 30°; (b) 45°;(c) 65°.

進一步分析松散體中相同比重的粗、細顆粒分別在堆積區顆粒累積質量供給中的占額, 如圖20所示. 這里定義粗、細顆粒的累積質量差值最早大于0.5 kg 的時間點作為二者累積質量占額的分異點. 圖20 表明, 當坡度為30°時, 粗、細顆粒在堆積區的顆粒累積質量供給差異很小, 無分異點; 隨著坡度增大, 粗、細顆粒的分異點出現時間越早; 達到穩定值時, 細顆粒的累積質量供給大于粗顆粒,但隨坡度增大, 二者間的差距有縮小趨勢. 主要是由于細顆粒摩擦力較小, 整體滑動性好; 而粗顆粒運動十分活躍, 除了消耗自身能量, 也會使得細顆粒獲得一部分能量.

3.5 顆粒堆積過程的動能與接觸力分布特性

為探究松散體滑動堆積特性的內在機理, 以坡度為65°、含石量為50%的松散體為例, 分析顆粒滑動過程中的動能與接觸力分布特性, 如圖21—圖25 所示. 顆粒平均動能與平均接觸力的均值見表3. 需要指出的是, 本文模擬的非球形(復合幾何體)顆粒是由四個球形顆粒組成, 在文獻[31]和文獻[32]中的研究結論中, 認為基于球形的接觸模型及其動能和接觸力的統計方法可以用于對復合幾何體顆粒的研究.

圖21 顆粒平均動能分布特性 (a) 平動動能; (b) 轉動動能Fig. 21. Granular average kinetic energy distribution characteristics: (a)Translational kinetic energy; (b) rotational kinetic energy.

結合圖21 與表3 可以看出, 粗顆粒的平均動能波動幅度最大, 說明粗顆粒的運動十分活躍, 從而為細顆粒的下滲提供足夠的間隙, 當細顆粒相對位移后, 粗顆粒就變得更容易發生滾動或滑動. 顆粒的平動動能遠遠大于轉動動能, 并且粗顆粒的平均動能及變化幅度均大于細顆粒的平均動能. 其中, 粗顆粒的平均平動動能均值最大, 達到了細顆粒的70.2 倍, 粗顆粒的平均轉動動能均值達到了細顆粒的5.6 倍. 這就從能量角度解釋了粗顆粒在堆積運動中較細顆粒更活躍的原因, 粗顆粒較細顆粒具有更大的能量, 動能衰減時間更長.

圖22、圖23 給出的分別是顆粒相互作用時的平均法向接觸力與平均切向接觸力時程曲線. 結合表3 分析可知, 顆粒間的法向力均為切向力的兩倍多, z 方向的法向力與切向力最大, y 方向次之,x 方向最小; 這符合顆粒運動主要沿水平和豎直運動, 同時也說明顆粒的橫向運動不顯著. 可以觀察到, 粗顆粒間的接觸力波動最劇烈, 其平均接觸力均值約為細顆粒與粗顆粒間平均接觸力均值的4.8 倍, 約為細顆粒間平均接觸力均值的8.1 倍. 說明粗顆粒間的相互碰撞作用更明顯, 滾動和轉動運動十分活躍, 使得接觸力及能量消耗較大. 細顆粒與粗顆粒的碰撞接觸力主要是來自二者的相對位移運動, 而它們的接觸力變化及力的傳遞形式較為復雜, 還有待進一步研究. 細顆粒間的接觸力最穩定, 這就反映了細顆粒在整個運動中以滑動為主.

圖24 給出的是顆粒間平均法向接觸力重疊量與平均切向接觸力重疊量時程曲線. 可以觀察到,粗顆粒間的重疊量最大, 細顆粒與粗顆粒次之, 細顆粒間最小; 根據(4)式, 顆粒間的接觸力隨接觸力重疊量a 的增大而增大. 當大部分顆粒到達底板時, 接觸重疊量達到最大, 法向重疊量約為切向重疊量的三倍多; 粗顆粒間的接觸重疊量均為細顆粒間、細顆粒與粗顆粒間的兩倍多; 說明顆粒的位移運動主要來自顆粒間接觸的法向力. 最后穩定階段粗顆粒間的接觸重疊量最大, 驗證了圖17 的結論, 粗顆粒較細顆粒的平均位移能力更強.

圖22 顆粒間平均法向接觸力時程曲線 (a) x 方 向;(b) y 方向; (c) z 方向Fig. 22. Time-history curve of average normal contact force between granulars: (a) x direction; (b) y direction; (c) z direction.

圖23 顆粒間平均切向接觸力時程曲線 (a) x 方向;(b) y 方向; (c) z 方向Fig. 23. Time-history curve of average tangential contact force between granulars: (a) x direction; (b) y direction;(c) z direction.

為了進一步從顆粒不同方向相互作用的細觀接觸力分布特性角度出發, 深入分析散體顆粒的滑動堆積特性, 故采用概率密度函數(probability density function, PDF)進行定量研究. 由于顆粒的切向接觸力與法向接觸力變化規律基本一致, 因此這里以分析法向接觸力為主. 由圖25 分析可知,法向接觸力越大或越小時, 概率密度越小, 而在0 附近的接觸力概率密度越大; 可以觀察到, 相比于粗顆粒間的作用力, 細顆粒與粗顆粒間、細顆粒間的接觸力概率密度逐漸增大, 但力的分布寬度逐漸減小. 此外, 從接觸力概率分布的復雜程度來看,細顆粒間的概率密度分布比較穩定, x 方向和y 方向的粗顆粒間、細顆粒與粗顆粒間的概率密度分別在0 附近對稱分布; 而z 方向則表現出粗顆粒間的概率密度均分布在0 的右側, 且隨著接觸力值的增大, 概率密度越小, 細顆粒與粗顆粒間的概率密度均分布在0 的左側, 隨著接觸力的減小, 概率密度越小. 上述數據分析表明了一個有趣的可能, 在更大的松散顆粒體系滑動運動中, 細顆粒間的接觸力長度較短或趨于點狀并位于顆粒流底層; 細顆粒與粗顆粒間的接觸力形式變化復雜, 可能存在耦合作用, 使得一定的細顆粒作用在粗顆粒周圍; 在z 方向上的粗顆粒間接觸力或許存在一個臨界值, 接觸力小于臨界值的概率密度呈指數函數分布, 接觸力大于臨界值的概率密度則呈冪函數分布.

表3 滑動堆積過程中顆粒的平均動能和接觸力均值Table 3. Average kinetic energy and contact force of granular in the process of sliding accumulation.

圖24 顆粒間平均接觸力重疊量時程曲線 (a)法向;(b) 切向Fig. 24. Time-history curve of average contact force overlap between granulars: (a) Normal; (b) tangential.

3.6 小 結

將上述含石量和坡度對松散體滑動后堆積特性的影響規律匯總, 見表4.

圖25 顆粒間平均法向接觸力概率密度函數(PDF)分布 (a) x 方向; (b) y 方向; (c) z 方向Fig. 25. Probability density functions (PDF) of average normal contact force between granulars: (a) x direction; (b) y direction; (c) z direction.

表4 模擬結果匯總表Table 4. Summary table of simulation results.

4 結 論

本文基于DEM, 探究了不同含石量和坡度對松散體滑動后堆積特性的影響, 主要結論如下:

1)相同坡度, 隨著含石量增加, 松散體最終堆積區的沖程、堆積寬度和堆積面積均表現為先增大后減小, 而最大厚度則是先減小后增大. 含石量對堆積區輪廓平面形狀的影響較小, 最終累積質量隨含石量的增加而減小.

2)相同含石量, 隨著坡度增大, 松散體最終堆積區的沖程、堆積寬度、堆積面積和累積質量均會變大, 而最大厚度近似線性減小. 坡度對沖程、堆積面積和平面堆積輪廓形狀的影響顯著, 平面堆積輪廓由近似圓狀趨于近似橢圓狀.

3)相比于含石量, 坡度對靜堆積角的影響更加顯著, 且隨坡度的增大而減小, 顆粒的流動性越佳. 含石量、坡度與靜堆積角的擬合關系模型表明,在一定條件下, 可為靜堆積角對含石量、坡度的響應提供預測模型.

4)堆積區體積計算中, 粗、細顆粒的體積占額存在一個臨界距離Lc, 當距坡腳線距離LLc時, 粗顆粒反而較細顆粒體積大; 坡度越大,Lc值越大. 堆積區累積質量計算中, 坡度越大, 粗、細顆粒累積質量占額分異點出現越早.

5)最后, 對顆粒滑動與堆積過程的平動動能和轉動動能進行了討論, 詳細分析了顆粒間不同方向、不同顆粒間的接觸力及其概率密度分布特性;從而在細觀尺度上探討了松散顆粒滑動與堆積運動的內在機理.

本文的研究結論是在組合球模擬非球形顆粒的基礎上所得出的. 由于研究重點主要是含石量和坡度變化對松散顆粒物理堆積特性的影響, 因此并未深入考慮顆粒形狀因素對研究結果的影響. 顆粒形狀特別是研究復雜顆粒對其物理堆積特性、動力學行為及接觸力學的影響, 一直是作者研究的重點方向之一, 更多的研究工作正在進行中, 以期研究結論更具普適性, 進而為顆粒的物理機制揭示工程地質行為提供更為現實的參考價值.