成浩 韓培鋒2)? 蘇有文

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

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

1 引 言

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

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

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

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

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

2 顆?；瑒佣逊e的離散元模型

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

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

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

式中, ?t是時間步長;N對應(yīng)時間t.

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

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

圖1 顆粒計算迭代示意圖 (a) 力與位移關(guān)系; (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]

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

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

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

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

圖2 顆粒相互作用時的受力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為顆粒間接觸面為圓形時的接觸半徑.

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

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

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

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

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

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

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

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

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

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

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

圖6 顆?；瑒雍罄鄯e體積試驗值[30]與DEM 計算值比較Fig. 6. Comparison of cumulative volume between experiment results[30] and DEM simulation results of granular after sliding.

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

3.1 堆積過程與速度分布

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

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

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

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

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

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

圖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 最終平面堆積形態(tài)示意圖Fig. 8. The final diagram of plane accumulation form.

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

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

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

圖9 含石量對堆積形態(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 不同計算條件下靜堆積角測量值Table 2. Measured value of static accumulation angle under different computing conditions.

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

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

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

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

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

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

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

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

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

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

圖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 含石量和坡度對堆積區(qū)顆粒質(zhì)量補給的影響

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

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

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

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

圖20 含石量50%時各坡度顆粒累積質(zhì)量對比 (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 的時間點作為二者累積質(zhì)量占額的分異點. 圖20 表明, 當(dāng)坡度為30°時, 粗、細(xì)顆粒在堆積區(qū)的顆粒累積質(zhì)量供給差異很小, 無分異點; 隨著坡度增大, 粗、細(xì)顆粒的分異點出現(xiàn)時間越早; 達(dá)到穩(wěn)定值時, 細(xì)顆粒的累積質(zhì)量供給大于粗顆粒,但隨坡度增大, 二者間的差距有縮小趨勢. 主要是由于細(xì)顆粒摩擦力較小, 整體滑動性好; 而粗顆粒運動十分活躍, 除了消耗自身能量, 也會使得細(xì)顆粒獲得一部分能量.

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

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

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

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

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

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

圖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.

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

表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 小 結(jié)

將上述含石量和坡度對松散體滑動后堆積特性的影響規(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, 探究了不同含石量和坡度對松散體滑動后堆積特性的影響, 主要結(jié)論如下:

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

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

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

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

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

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