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

成浩 韓培鋒2)? 蘇有文

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

2) (水利部山洪地質(zhì)災(zāi)害防治工程技術(shù)研究中心, 武漢 430010)

1 引 言

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

學(xué)者們?cè)谘芯克缮Ⅲw滑動(dòng)后的堆積特性時(shí), 主要考慮初始松散體的體積、顆粒粒徑、坡高、啟動(dòng)區(qū)坡度、坡角、摩擦系數(shù)和地震荷載等因素對(duì)松散顆粒變形過(guò)程和滑動(dòng)后的堆積形態(tài)、沖程、平均粒徑分布、固相分?jǐn)?shù)及分選系數(shù)等指標(biāo)的影響規(guī)律[6?10].堆積形態(tài)與顆粒運(yùn)動(dòng)過(guò)程及速度密切相關(guān), 通過(guò)物理模型試驗(yàn)和理論分析提出了滑坡沖程與速度的預(yù)測(cè)方法[11]、落石運(yùn)動(dòng)速度與距離的數(shù)學(xué)公式[12]、松散體運(yùn)動(dòng)中的障礙物與沖程的動(dòng)態(tài)關(guān)系模型[13].可見, 模型試驗(yàn)和理論分析為其他方法研究松散顆?；瑒?dòng)堆積特性提供了前提, 但缺乏滑動(dòng)過(guò)程的實(shí)測(cè)數(shù)據(jù)去解釋松散顆粒滑動(dòng)堆積特性的內(nèi)在機(jī)理.此外, 在現(xiàn)場(chǎng)調(diào)研中發(fā)現(xiàn)坡度對(duì)松散體滑動(dòng)堆積范圍分布各異, 松散體中的含石量大小與滑后堆積特性之間的內(nèi)在關(guān)系仍需進(jìn)一步定量研究.

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

松散體滑動(dòng)后堆積特性的內(nèi)在機(jī)理與運(yùn)動(dòng)中的顆粒能量及動(dòng)力學(xué)關(guān)系密切. 根據(jù)牛頓第二定律及顆粒相互作用時(shí)的接觸力關(guān)系, 可從顆粒流分析法和能量交換方面分析顆粒整個(gè)運(yùn)動(dòng)過(guò)程[20,21]. 顆粒間的動(dòng)勢(shì)能轉(zhuǎn)化與顆粒的摩擦聯(lián)系緊密, 可用土力學(xué)方法分析顆粒在運(yùn)動(dòng)中的能量衰減[22,23]. 從上述研究中可以看出, 顆粒間的能量及接觸力為解釋松散顆粒滑動(dòng)堆積特性提供了切入點(diǎn).

因此, 本文采用DEM 構(gòu)建了松散體滑動(dòng)堆積三維模型, 對(duì)非球形顆?；瑒?dòng)堆積過(guò)程進(jìn)行了數(shù)值模擬, 研究了含石量和坡度變化對(duì)松散體滑動(dòng)堆積特性的影響規(guī)律, 分析了顆粒滑動(dòng)堆積過(guò)程中的能量交換和接觸作用力, 為深入理解和進(jìn)一步揭示松散顆粒堆積特性的內(nèi)在機(jī)理提供參考.

2 顆?；瑒?dòng)堆積的離散元模型

2.1 DEM 的計(jì)算原理[14?16,24]

DEM 認(rèn)為顆粒單元間的運(yùn)動(dòng)是相互獨(dú)立的,只有顆粒單元間發(fā)生相互作用時(shí), 才會(huì)在接觸點(diǎn)處產(chǎn)生作用力. 根據(jù)力與位移關(guān)系, 由位移得到顆粒受到的作用力, 再由牛頓第二定律得出顆粒i的運(yùn)動(dòng)方程

利用中心差分法對(duì)(1)式進(jìn)行積分, 得到以兩次迭代時(shí)間步長(zhǎng)的中間點(diǎn)表示的更新速度表達(dá)式為

式中, ?t是時(shí)間步長(zhǎng);N對(duì)應(yīng)時(shí)間t.

對(duì)(2)式進(jìn)行積分, 可得到關(guān)于位移的等式表述為

進(jìn)而得到顆粒單元新的位移值, 并將該新位移代入力與位移關(guān)系計(jì)算新的作用力, 循環(huán)反復(fù)就可實(shí)時(shí)跟蹤每個(gè)顆粒單元在任意時(shí)刻的運(yùn)動(dòng), 如圖1所示.

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

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

當(dāng)顆粒與顆粒發(fā)生相互作用時(shí), 會(huì)產(chǎn)生相應(yīng)的接觸法向力Fn和切向力Ft, 如圖2 所示.

圖3 為顆粒法向力示意圖,Fn表達(dá)式為

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

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

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

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

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

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

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

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

2.3 數(shù)值模擬計(jì)算方案

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

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

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

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

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

2.4 顆?；瑒?dòng)堆積過(guò)程的試驗(yàn)對(duì)比驗(yàn)證

為了驗(yàn)證本文數(shù)值計(jì)算的準(zhǔn)確性與穩(wěn)定性, 首先對(duì)參數(shù)為坡度為38°, 體積為0.015 m3, 顆粒粒徑為5—20 mm, 單個(gè)顆粒形狀如圖5 所示, 其他參數(shù)見表1 的松散體滑動(dòng)過(guò)程進(jìn)行了數(shù)值模擬計(jì)算, 并與文獻(xiàn)[30]中的物理模型試驗(yàn)結(jié)果進(jìn)行對(duì)比, 如圖6 所示. 可知, DEM 計(jì)算結(jié)果累積體積開始的時(shí)間略早, 最后達(dá)到穩(wěn)定階段的累積體積較試驗(yàn)結(jié)果稍小, 可能是因?yàn)樵囼?yàn)中滑槽的摩擦系數(shù)偏小, 但總體變化的規(guī)律與試驗(yàn)結(jié)果吻合良好, 說(shuō)明顆粒形狀對(duì)累積體積的影響較小, 同時(shí)反映了本文DEM 數(shù)值模型可以較好地模擬松散體的滑動(dòng)過(guò)程.

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

圖6 顆?；瑒?dòng)后累積體積試驗(yàn)值[30]與DEM 計(jì)算值比較Fig. 6. Comparison of cumulative volume between experiment results[30] and DEM simulation results of granular after sliding.

3 數(shù)值模擬結(jié)果分析

3.1 堆積過(guò)程與速度分布

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

綜上, 松散體滑動(dòng)堆積大體可劃分為啟程變形加速, 近程高速滑動(dòng), 中程減速滑動(dòng), 遠(yuǎn)程準(zhǔn)靜態(tài)堆積四個(gè)階段; 這里只是定性的分析松散體堆積過(guò)程中長(zhǎng)度、厚度、形狀及速度的變化, 接下來(lái)將對(duì)堆積結(jié)果進(jìn)行定量分析.

3.2 含石量對(duì)堆積特性的影響

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

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

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

圖7 松散體滑動(dòng)堆積全過(guò)程(+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) 堆積過(guò)程形狀變化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 最終平面堆積形態(tài)示意圖Fig. 8. The final diagram of plane accumulation form.

由表2 分析可知, 隨著含石量增大, 在一定坡度時(shí), 靜堆積角會(huì)出現(xiàn)小幅度減小, 而超過(guò)一定坡度時(shí), 靜堆積角會(huì)小幅度增大. 由于靜堆積角與松散顆粒的流動(dòng)性密切相關(guān), 從而反映了松散體中含石量值的大小對(duì)顆粒流動(dòng)性有不同程度的影響, 坡度較小時(shí)顆粒流動(dòng)性有增強(qiáng)的趨勢(shì), 而坡度較大時(shí)顆粒流動(dòng)性則有所減弱; 主要是因?yàn)楹孔兓瘯?huì)使得顆?；瑒?dòng)過(guò)程中的碰撞強(qiáng)度、摩擦及接觸形式在一定程度上改變顆粒的流動(dòng)性, 除此之外影響顆粒流動(dòng)性的因素還有很多, 如黏聚力、內(nèi)摩擦力和初始堆積體密度等, 所以未來(lái)的工作會(huì)針對(duì)這一點(diǎn)進(jìn)一步研究. 另外, 相比于含石量, 坡度對(duì)靜堆積角值的影響顯著; 換句話說(shuō), 坡度對(duì)顆粒流動(dòng)性的影響關(guān)系更加密切, 至于坡度是如何影響靜堆積角和顆粒流動(dòng)性的, 以及二者的關(guān)系模型在后文會(huì)進(jìn)一步分析.

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

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

圖9 含石量對(duì)堆積形態(tài)的影響 (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 不同計(jì)算條件下靜堆積角測(cè)量值Table 2. Measured value of static accumulation angle under different computing conditions.

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

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

圖11 堆積區(qū)體積計(jì)算平面示意圖Fig. 11. Schematic diagram of volume calculation of accumulation area.

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

3.3 坡度對(duì)堆積特性的影響

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

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

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

圖14 表明, 隨著坡度增大, 堆積輪廓平面形狀由近似圓狀趨于近似橢圓狀, 主要是因?yàn)轭w粒堆沖程的延伸性較堆積寬度顯著. 同時(shí)也說(shuō)明坡度比含石量對(duì)堆積輪廓平面形狀影響明顯.

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

圖13 坡度對(duì)堆積形態(tài)的影響 (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%的松散體在不同坡度下滑動(dòng)后的體積計(jì)算結(jié)果. 分析可知, 當(dāng)坡度為30°和65°時(shí), 距離坡腳線越遠(yuǎn), 顆粒堆體積逐漸減小, 顆粒越稀疏; 坡度45°時(shí), 顆粒堆體積在坡腳線附近形成隆起區(qū), 當(dāng)距離持續(xù)增大, 顆粒堆體積依然逐漸減小.

圖15 含石量50%時(shí)不同坡度松散顆粒滑動(dòng)堆積模擬結(jié)果 (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%時(shí)不同坡度對(duì)堆積體積的影響Fig. 16. Influence of different slope of 50% stone content on the accumulation volume.

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

圖17 含石 量50%時(shí)不同坡度顆粒 體積 對(duì)比 (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 含石量和坡度對(duì)堆積區(qū)顆粒質(zhì)量補(bǔ)給的影響

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

圖19 給出的是同一含石量下, 不同坡度的松散體對(duì)堆積區(qū)顆粒累積質(zhì)量供應(yīng)時(shí)程曲線. 分析可知, 坡度越小的松散體對(duì)堆積區(qū)的顆粒質(zhì)量供給開始時(shí)間越晚, 且供給持續(xù)時(shí)間越長(zhǎng); 隨著坡度增大,累積質(zhì)量曲線的增長(zhǎng)率越大. 在達(dá)到穩(wěn)定值時(shí), 坡度越大, 堆積區(qū)的顆粒累積質(zhì)量越大; 其中, 坡度為30°—45°的顆粒累積質(zhì)量增幅十分顯著, 當(dāng)坡度持續(xù)增大時(shí), 增幅趨于平緩.

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

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

圖20 含石量50%時(shí)各坡度顆粒累積質(zhì)量對(duì)比 (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°.

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

3.5 顆粒堆積過(guò)程的動(dòng)能與接觸力分布特性

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

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

結(jié)合圖21 與表3 可以看出, 粗顆粒的平均動(dòng)能波動(dòng)幅度最大, 說(shuō)明粗顆粒的運(yùn)動(dòng)十分活躍, 從而為細(xì)顆粒的下滲提供足夠的間隙, 當(dāng)細(xì)顆粒相對(duì)位移后, 粗顆粒就變得更容易發(fā)生滾動(dòng)或滑動(dòng). 顆粒的平動(dòng)動(dòng)能遠(yuǎn)遠(yuǎn)大于轉(zhuǎn)動(dòng)動(dòng)能, 并且粗顆粒的平均動(dòng)能及變化幅度均大于細(xì)顆粒的平均動(dòng)能. 其中, 粗顆粒的平均平動(dòng)動(dòng)能均值最大, 達(dá)到了細(xì)顆粒的70.2 倍, 粗顆粒的平均轉(zhuǎn)動(dòng)動(dòng)能均值達(dá)到了細(xì)顆粒的5.6 倍. 這就從能量角度解釋了粗顆粒在堆積運(yùn)動(dòng)中較細(xì)顆粒更活躍的原因, 粗顆粒較細(xì)顆粒具有更大的能量, 動(dòng)能衰減時(shí)間更長(zhǎng).

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

圖24 給出的是顆粒間平均法向接觸力重疊量與平均切向接觸力重疊量時(shí)程曲線. 可以觀察到,粗顆粒間的重疊量最大, 細(xì)顆粒與粗顆粒次之, 細(xì)顆粒間最小; 根據(jù)(4)式, 顆粒間的接觸力隨接觸力重疊量a 的增大而增大. 當(dāng)大部分顆粒到達(dá)底板時(shí), 接觸重疊量達(dá)到最大, 法向重疊量約為切向重疊量的三倍多; 粗顆粒間的接觸重疊量均為細(xì)顆粒間、細(xì)顆粒與粗顆粒間的兩倍多; 說(shuō)明顆粒的位移運(yùn)動(dòng)主要來(lái)自顆粒間接觸的法向力. 最后穩(wěn)定階段粗顆粒間的接觸重疊量最大, 驗(yàn)證了圖17 的結(jié)論, 粗顆粒較細(xì)顆粒的平均位移能力更強(qiáng).

圖22 顆粒間平均法向接觸力時(shí)程曲線 (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 顆粒間平均切向接觸力時(shí)程曲線 (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.

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

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

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

3.6 小 結(jié)

將上述含石量和坡度對(duì)松散體滑動(dòng)后堆積特性的影響規(guī)律匯總, 見表4.

圖25 顆粒間平均法向接觸力概率密度函數(shù)(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 模擬結(jié)果匯總表Table 4. Summary table of simulation results.

4 結(jié) 論

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

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

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

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

4)堆積區(qū)體積計(jì)算中, 粗、細(xì)顆粒的體積占額存在一個(gè)臨界距離Lc, 當(dāng)距坡腳線距離LLc時(shí), 粗顆粒反而較細(xì)顆粒體積大; 坡度越大,Lc值越大. 堆積區(qū)累積質(zhì)量計(jì)算中, 坡度越大, 粗、細(xì)顆粒累積質(zhì)量占額分異點(diǎn)出現(xiàn)越早.

5)最后, 對(duì)顆?；瑒?dòng)與堆積過(guò)程的平動(dòng)動(dòng)能和轉(zhuǎn)動(dòng)動(dòng)能進(jìn)行了討論, 詳細(xì)分析了顆粒間不同方向、不同顆粒間的接觸力及其概率密度分布特性;從而在細(xì)觀尺度上探討了松散顆?；瑒?dòng)與堆積運(yùn)動(dòng)的內(nèi)在機(jī)理.

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