基于DEM-CFD耦合方法的煤系土邊坡失穩(wěn)機(jī)理宏細(xì)觀分析

張 鴻，張 榜，豐浩然，吳 燦

(1.南昌工程學(xué)院 江西省水利土木工程基礎(chǔ)設(shè)施安全重點(diǎn)實(shí)驗(yàn)室，江西 南昌 330099；2.南昌工程學(xué)院 土木與建筑工程學(xué)院，江西 南昌 330099)

煤系地層邊坡開挖后，風(fēng)化作用加強(qiáng)，當(dāng)雨水滲入使得煤系土快速崩解軟化，力學(xué)強(qiáng)度急速降低，邊坡極易發(fā)生坍塌滑坡現(xiàn)象。目前，對(duì)煤系土的研究較少，尚處于起步階段，而且現(xiàn)有的研究主要集中在對(duì)煤系土不良地質(zhì)邊坡的加固處理方法[1-3]、邊坡穩(wěn)定性分析[4-6]、以及煤系土強(qiáng)度特性等宏觀層面上[7-9]，從細(xì)觀角度對(duì)煤系土邊坡失穩(wěn)的細(xì)觀機(jī)理進(jìn)行分析與解釋的文獻(xiàn)未有報(bào)道。煤系土邊坡的失穩(wěn)破壞運(yùn)動(dòng)是一個(gè)存在巖土體滑動(dòng)、平移、轉(zhuǎn)動(dòng)的復(fù)雜過(guò)程，具有宏觀上的不連續(xù)性和單個(gè)塊體運(yùn)動(dòng)的隨機(jī)性，用傳統(tǒng)的宏觀連續(xù)介質(zhì)理論已經(jīng)不能合理地分析散體過(guò)程，上述煤系土邊坡風(fēng)化之后的散、動(dòng)特征都與傳統(tǒng)的均勻、連續(xù)等假定沖突，將導(dǎo)致理論與實(shí)際的偏離。因此，要想揭示煤系土邊坡的失穩(wěn)機(jī)理，有必要從微細(xì)觀角度進(jìn)行研究[10]。

降雨入滲是邊坡發(fā)生失穩(wěn)破壞的主要誘發(fā)因素之一，目前考慮流固耦合的計(jì)算方法主要有3類，第1類是Euler-Euler方法，該方法是將流體和固體都視為連續(xù)介質(zhì)；第2類是Lagrange-Lagrange方法，該方法是把流體和固體都假設(shè)為離散介質(zhì)；第3類是Euler-Lagrange計(jì)算方法，該方法將流體視為連續(xù)介質(zhì)，采用流體動(dòng)力學(xué)方法（computational fluid dynamic，CFD）進(jìn)行模擬計(jì)算，固體視為離散介質(zhì)，采用離散單元方法（discrete element method，DEM）進(jìn)行模擬計(jì)算。這種DEM-CFD流固耦合方法早期主要應(yīng)用于工業(yè)領(lǐng)域，Tsuji等[11]首次將CFD方法和DEM方法進(jìn)行耦合，對(duì)固體顆粒中的氣體流動(dòng)行為進(jìn)行模擬，仿真計(jì)算結(jié)果很好地展現(xiàn)了DEM-CFD耦合方法模擬流體-固體相互運(yùn)動(dòng)問題的能力。Xu等[12]提出了動(dòng)力碰撞模型，將CFD和DEM的時(shí)間與空間尺度設(shè)置為一致，采用牛頓第三定律來(lái)描述流體與顆粒的相互作用。在此基礎(chǔ)上，Xu[13]、Yu[14]等進(jìn)一步將他們的方法進(jìn)行了改進(jìn)，對(duì)介于固體和流體之間的運(yùn)動(dòng)進(jìn)行了模擬。

Shamy等[15-16]首次將DEM-CFD耦合方法引入到巖土工程問題分析中，將Navier-Stokes方程進(jìn)行了簡(jiǎn)化處理，采用3維DEM-CFD耦合方法計(jì)算分析了邊坡滲流問題，對(duì)土體邊坡滲流的細(xì)觀機(jī)理進(jìn)行了分析。國(guó)內(nèi)一些研究人員在土體滲流和土體液化研究中也采用了DEM-CFD流固耦合的計(jì)算方法，取得了較理想的研究效果[17-21]。蔣明鏡等[22]在2維離散單元軟件（PFC2D）的基礎(chǔ)上，將Tait狀態(tài)方程引入到CFD，建立了弱可壓縮流體DEM-CFD流固耦合模型，并通過(guò)模擬計(jì)算單個(gè)顆粒水下自由沉降和單面固結(jié)排水試驗(yàn)來(lái)驗(yàn)證所建耦合模型的可行性。Khalili等[23]將離散化的DEM-CFD控制方程嵌入PFC2D軟件中，對(duì)雙軸不排水壓縮試驗(yàn)進(jìn)行了數(shù)值模擬計(jì)算。王胤等[24]利用開源程序CFDEM建立了DEM-CFD流固耦合數(shù)值計(jì)算模型，并將顆粒滾動(dòng)阻力機(jī)制引入到CFD-DEM耦合模型中，計(jì)算分析了滾動(dòng)阻力模型的有效性以及所建立耦合模型的可行性。從以上文獻(xiàn)可以發(fā)現(xiàn)，DEM-CFD耦合方法是一種非常有發(fā)展前景的數(shù)值模擬方法，特別是在土體與流體相互耦合作用的微細(xì)觀分析方面具有明顯的優(yōu)勢(shì)。

作者以江西省萬(wàn)載至宜春高速公路項(xiàng)目的煤系土邊坡為研究背景，將離散元方法（DEM）與計(jì)算流體動(dòng)力學(xué)方法（CFD）進(jìn)行耦合，建立了煤系土邊坡3維DEM-CFD流固耦合細(xì)觀作用計(jì)算模型，分析降雨作用下煤系土邊坡失穩(wěn)破壞的細(xì)觀機(jī)理，計(jì)算結(jié)果與模型試驗(yàn)結(jié)果比較吻合，研究成果可以為該區(qū)域煤系土邊坡的防護(hù)設(shè)計(jì)與施工提供理論依據(jù)。

1 DEM-CFD流固耦合計(jì)算基本原理

1.1 流體-顆粒相互作用

采用流-固耦合理論（DEM-CFD）分析流體對(duì)巖土顆粒的作用，是將流體視為連續(xù)體，將巖土體視為不連續(xù)的顆粒（Euler-Lagrange法）[25]。在CFD中流體的速度v是宏觀的速度，即單位橫截面積上的體積流量，假定流體在整個(gè)橫截面上產(chǎn)生，而實(shí)際流體流速只發(fā)生在顆?？紫犊臻g中。作者將流體與顆粒之間的相互作用力進(jìn)行了簡(jiǎn)化，僅考慮流體對(duì)顆粒的拖曳力和流體壓力梯度力，流體與顆粒間作用方程如下[26]：

式（1）～（2）中，u為顆粒的速度，fmech為作用于顆粒的外力之和，ffluid為流體對(duì)顆粒的作用力（包括流體對(duì)顆粒的拖曳力和流體壓力梯度力），m為顆粒的質(zhì)量，g為重力加速度， ω為顆粒旋轉(zhuǎn)角速度，M為作用于顆粒上的力矩，I為慣性矩。

假設(shè)流體施加于顆粒的力總是作用于形心，不考慮彎矩，則拖曳力fdrag為：

式中：f0為 單個(gè)顆粒所受的拖曳力；ε為流體所在單元的孔隙度； ε-χ項(xiàng)為考慮局部孔隙度的經(jīng)驗(yàn)系數(shù)，它可以使拖拽力同時(shí)適用于高孔隙度和低孔隙度的情況。

式（3）中單個(gè)顆粒所受的拖曳力為：

式中， ρf為流體密度，r為顆粒半徑，u為 顆粒的速度，v為流體速度， ρf為流體的動(dòng)力黏滯系數(shù)，Cd為拖拽力系數(shù)，其表示為：

式（3）中的經(jīng)驗(yàn)系數(shù) χ為：

式中，Rep為顆粒的雷諾數(shù)，定義為：

施加在流體單位體積上的體力為：

式中，V為流體單元體積， ?p為流體梯度，分子上求和對(duì)象為流體單元疊加的顆粒。

因此，流體對(duì)顆粒的相互作用力為：

1.2 流體-流體相互作用

1）連續(xù)方程

平均的Navier-Stokes方程可以定量地描述孔隙流體特性，可用于流體的細(xì)觀分析。流體的連續(xù)方程是根據(jù)質(zhì)量守恒定律進(jìn)行推導(dǎo)得出，如式（10）所示，其表達(dá)的物理意義為：?jiǎn)挝粫r(shí)間內(nèi)流入體積元的質(zhì)量與流出體積元的質(zhì)量差值等于體積元內(nèi)質(zhì)量的增加或減少。

式中，?為拉普拉斯算子，ε為流體所在單元的孔隙度。

2）動(dòng)量方程

牛頓第二定律流體流動(dòng)模型所獲得的控制方程推導(dǎo)得到動(dòng)量方程。Anderson等[27]將外力在體積元內(nèi)取平均，推導(dǎo)出流體的動(dòng)量方程為：

式中，p為流體壓力。

1.3 顆粒運(yùn)動(dòng)控制方程

顆粒的運(yùn)動(dòng)方程是在牛頓第二定律的基礎(chǔ)上推導(dǎo)出來(lái)的，為了考慮流體對(duì)顆粒的作用，在控制方程中引入顆粒與流體間的相互作用力ffluid，得到：

式（12）～（13）中，ui為第i個(gè)顆粒的速度，fij為作用于顆粒i的接觸力，fg(i)為第i個(gè)顆粒的重力，Ii為顆粒的轉(zhuǎn)動(dòng)慣量，Mij為 第i個(gè)顆粒受到第j個(gè)顆粒的轉(zhuǎn)動(dòng)力矩。

1.4 CFD-DEM耦合流程與程序?qū)崿F(xiàn)

本文以離散單元法商業(yè)軟件PFC3D為基礎(chǔ)，采用Fish語(yǔ)言編寫相關(guān)程序，使PFC3D與其內(nèi)嵌的CFD求解器進(jìn)行數(shù)據(jù)交換，從而實(shí)現(xiàn)流-固耦合計(jì)算。具體計(jì)算過(guò)程為：首先，根據(jù)邊界條件，利用有限體積法將CFD模塊中的流體控制方程式（10）、（11）進(jìn)行離散，并采用PISO應(yīng)力速度耦合算法進(jìn)行求解[28]，與此同時(shí)，將網(wǎng)格數(shù)據(jù)發(fā)送至DEM模塊，通過(guò)DEM模塊計(jì)算孔隙度和拖曳力并將數(shù)據(jù)發(fā)回CFD模塊；然后，再將每一時(shí)間步控制單元內(nèi)的流體與顆粒之間作用力發(fā)送至DEM計(jì)算模塊，再運(yùn)用離散元方法將流體作用力施加在土顆粒上，同時(shí)進(jìn)行顆粒之間的力學(xué)計(jì)算；最后，將流體作用力及孔隙度傳回CFD模型，以此循環(huán)計(jì)算直至結(jié)束。DEM-CFD耦合計(jì)算流程如圖1所示，當(dāng)兩者計(jì)算時(shí)間相等時(shí)觸發(fā)數(shù)據(jù)交互，完成流體與顆粒相互作用的計(jì)算過(guò)程。

圖1 DEM-CFD耦合計(jì)算流程Fig.1 DEM-CFD coupling calculation flow

2 煤系土邊坡失穩(wěn)DEM-CFD耦合數(shù)值模擬

2.1 煤系土邊坡室外模型試驗(yàn)

為了研究降雨入滲條件下煤系土邊坡失穩(wěn)破壞機(jī)理，項(xiàng)目組進(jìn)行了人工降雨室外邊坡模型試驗(yàn)。試驗(yàn)邊坡模型箱尺寸為4.5 m×3.0 m×2.4 m（長(zhǎng)×寬×高），試驗(yàn)土樣取自江西省萬(wàn)載至宜春高速公路項(xiàng)目A5標(biāo)段K30+120處煤系土邊坡開挖出露的原狀土，邊坡模型每40 cm填筑一層土體，坡率為1.0∶1.5。本次人工降雨模擬的降雨強(qiáng)度為江西宜春地區(qū)7月和8月份最大降雨強(qiáng)度為（9.84×10-7m/s）。試驗(yàn)從2018年11月19日上午8點(diǎn)開始，一直持續(xù)至晚上18點(diǎn)結(jié)束，現(xiàn)場(chǎng)觀察邊坡的失穩(wěn)破壞模式主要是水土流失和雨水沖蝕，形成了不同深度的沖蝕溝，如圖2（a）所示。

本次試驗(yàn)在邊坡土體埋設(shè)了4行29個(gè)側(cè)面示蹤點(diǎn)，通過(guò)模型箱有機(jī)玻璃側(cè)面可以清楚地觀察示蹤點(diǎn)的運(yùn)動(dòng)情況，如圖2（b）所示。通過(guò)分析量測(cè)結(jié)果可知，在持續(xù)降雨過(guò)程中，邊坡土體的主要破壞形式為雨水沖刷，在坡面形成深淺不一的沖蝕溝，最大沖蝕溝達(dá)0.36 m，發(fā)生在靠近坡中偏上的位置。詳細(xì)的邊坡模型試驗(yàn)結(jié)果分析參見文獻(xiàn)[29]。

圖2 邊坡失穩(wěn)破壞試驗(yàn)現(xiàn)場(chǎng)Fig.2 Site of slope failure test

2.2 煤系土邊坡數(shù)值計(jì)算模型建立

2.2.1 模型參數(shù)標(biāo)定

作者采用數(shù)值模擬室內(nèi)三軸試驗(yàn)的方法進(jìn)行顆粒細(xì)觀參數(shù)的標(biāo)定。利用PFC3D軟件作為計(jì)算平臺(tái)，自編程序建立三軸試驗(yàn)數(shù)值計(jì)算模型，試驗(yàn)?zāi)Ｐ蛪w尺寸為高4.0 m，直徑2.0 m的三軸墻體模型。為了節(jié)約計(jì)算時(shí)間和提高分析效率，將實(shí)際的土顆粒尺寸進(jìn)行放大和集中處理。顆粒粒徑在50～90 mm間等概率隨機(jī)分布，生成的計(jì)算模型如圖3所示。

圖3 三軸試驗(yàn)數(shù)值模擬計(jì)算模型Fig.3 Numerical simulation model of triaxial test

根據(jù)煤系土的性質(zhì)，標(biāo)定與邊坡數(shù)值建模均選用接觸黏結(jié)模型，該黏結(jié)作用于顆粒間很小的接觸點(diǎn)，可承受法向和切向應(yīng)力，但不能傳遞彎矩[30]，這與煤系土體性質(zhì)比較接近。模型內(nèi)生成顆粒之后，通過(guò)改變顆粒的細(xì)觀參數(shù)反復(fù)試算，使得數(shù)值模擬的偏應(yīng)力-軸應(yīng)變曲線逼近室內(nèi)試驗(yàn)結(jié)果，如圖4所示。從圖4中可以看出，模擬曲線的變化趨勢(shì)及大小與室內(nèi)試驗(yàn)結(jié)果相比都較為接近，說(shuō)明該模型參數(shù)取值可以較好地用于模擬煤系土細(xì)觀分析。顆粒參數(shù)細(xì)觀標(biāo)定結(jié)果如表1所示，其結(jié)果將直接用于邊坡的數(shù)值模擬試驗(yàn)。

圖4 三軸試驗(yàn)細(xì)觀參數(shù)標(biāo)定Fig.4 Meso parameters calibrated by simulated triaxial test

表1 煤系土細(xì)觀參數(shù)Tab. 1 Microscopic parameters of coal bearing soil

2.2.2 邊坡計(jì)算模型建立

以室外人工降雨邊坡模型試驗(yàn)為研究對(duì)象，通過(guò)fish語(yǔ)言編程，建立了3維邊坡流固耦合（DEM-CFD）相互作用細(xì)觀計(jì)算模型，如圖5所示。圖5（a）為雨水作用下邊坡3維流固耦合計(jì)算模型，通過(guò)顆粒分組可以看出，整個(gè)模型顆粒分成兩組，上部藍(lán)色顆粒為雨水作用的飽和區(qū)域，生成2 646個(gè)顆粒；下部淺綠色顆粒為邊坡土體的非飽和區(qū)域，生成4 572個(gè)顆粒。

本次模擬的狀態(tài)為模型邊坡在人工降雨作用8 h后的工況，根據(jù)現(xiàn)場(chǎng)試驗(yàn)結(jié)果，雨水在模型邊坡土體內(nèi)形成了暫態(tài)飽和區(qū)，其滲透深度約為0.5 m，為方便建模及計(jì)算，飽和與非飽和土體的界限以直線代替曲線，土體滲流區(qū)設(shè)置見圖5（b）。考慮到雨水作用于坡體之后，徑流帶動(dòng)顆粒會(huì)在地面流動(dòng)，所以將流體網(wǎng)格沿坡體前方進(jìn)行了拓展，以符合實(shí)際情況，如圖5（b）所示。在非飽和區(qū)域，其顆粒的細(xì)觀參數(shù)按照第2.2.1節(jié)三軸試驗(yàn)標(biāo)定的參數(shù)進(jìn)行設(shè)置，由于飽和區(qū)的土體強(qiáng)度較非飽和土體強(qiáng)度低，故對(duì)飽和土體的黏結(jié)強(qiáng)度進(jìn)行了折減，其余參數(shù)保持一致。

圖5 降雨入滲下邊坡DEM-CFD耦合計(jì)算模型Fig.5 Fluid-solid coupling（DEM-CFD）calculation model of slope under rainfall infiltration

根據(jù)模型邊坡人工降雨的流體特征及降雨強(qiáng)度，對(duì)作用于邊坡飽和區(qū)域的流體參數(shù)按照表2所示的參數(shù)進(jìn)行設(shè)置。

表2 流體參數(shù)設(shè)置Tab. 2 Setting of fluid parameters in saturated soil

為使流體流動(dòng)方向與實(shí)際一致，Z軸的流速與X軸的流速之比設(shè)為1.0∶1.5，流體流動(dòng)的方向與室外降雨試驗(yàn)坡體的徑流一致。

2.3 計(jì)算結(jié)果分析

2.3.1 煤系土邊坡顆粒運(yùn)動(dòng)的宏觀分析

1）邊坡土顆粒運(yùn)動(dòng)軌跡計(jì)算結(jié)果分析

圖6為模擬試驗(yàn)結(jié)束后從計(jì)算模型側(cè)面觀察到的土顆粒移動(dòng)情況。從圖6（a）可以看出，顆粒的移動(dòng)主要發(fā)生在飽和區(qū)，非飽和區(qū)局部有顆粒移動(dòng)，邊坡整體會(huì)出現(xiàn)一個(gè)滑動(dòng)帶，這一現(xiàn)象可以由顆粒的移動(dòng)清晰地看出。通過(guò)獲取計(jì)算邊坡模型的中截面，運(yùn)用CAD將邊坡的滑動(dòng)帶進(jìn)行繪制，如圖6（b）所示。從圖6可以發(fā)現(xiàn)，上部土體發(fā)生的位移約為0.9 m，下部土體發(fā)生的位移約為1.7 m，同時(shí)形成一條近似直線狀的滑動(dòng)面。坡頂顆粒被沖刷，沖刷比較嚴(yán)重區(qū)域位于斜坡上部，沖刷深度在0.3～0.4 m之間，這與室外邊坡模型試驗(yàn)結(jié)果比較一致。

圖6 邊坡側(cè)面顆粒運(yùn)動(dòng)軌跡示意圖Fig.6 Schematic diagram of the overall particle movement of the slope

2）邊坡土體應(yīng)力變化計(jì)算結(jié)果分析

為了解邊坡土體在降雨入滲條件下內(nèi)部受力變化，以便更好地分析土體內(nèi)部的破壞機(jī)理，模擬計(jì)算時(shí)在土體內(nèi)部設(shè)置了17個(gè)測(cè)量圓。對(duì)邊坡施加流體作用后，通過(guò)模擬計(jì)算得到邊坡土體X、Y、Z3個(gè)方向應(yīng)力變化監(jiān)測(cè)數(shù)據(jù)，然后經(jīng)surfer軟件后處理生成云圖，計(jì)算結(jié)果如圖7～9所示。

通過(guò)比較各時(shí)刻邊坡應(yīng)力演變過(guò)程發(fā)現(xiàn)：當(dāng)模型計(jì)算到20 000時(shí)步后，如圖7所示，此時(shí)斜坡和坡頂上的應(yīng)力較坡體內(nèi)部應(yīng)力小，應(yīng)力集中主要發(fā)生在非飽和區(qū)，說(shuō)明飽和區(qū)顆粒受流體作用，整體應(yīng)力變小，穩(wěn)定性變差；當(dāng)模型計(jì)算到60 000時(shí)步后，坡頂和斜坡上部的應(yīng)力下降較多，坡腳部分應(yīng)力增加，說(shuō)明在這段時(shí)間坡體內(nèi)部相對(duì)比較穩(wěn)定，而飽和區(qū)土體遭到持續(xù)破壞，見圖8；當(dāng)模型計(jì)算到100 000時(shí)步后，此時(shí)邊坡各方向應(yīng)力變化主要發(fā)生在坡頂和坡腳，在坡腳處土體應(yīng)力增加，坡頂部分土體應(yīng)力減少幅度較大（圖9），說(shuō)明在雨水作用下，坡頂土體被大量沖刷至坡腳堆積，導(dǎo)致坡腳處顆粒大規(guī)模堆積。

圖7 計(jì)算至20 000時(shí)步邊坡各方向應(yīng)力變化（中截面）Fig.7 Cloud chart of stress change in each direction of slope soil when the model is calculated to 20 000

圖8 計(jì)算至60 000時(shí)步邊坡各方向應(yīng)力變化（中截面）Fig.8 Cloud chart of stress change in each direction of slope soil when the model is calculated to 60 000

圖9 計(jì)算至100 000時(shí)步邊坡各方向應(yīng)力變化（中截面）Fig.9 Cloud chart of stress change in each direction of slope soil when the model is calculated to 100 000

2.3.2 煤系土邊坡失穩(wěn)破壞的細(xì)觀機(jī)理分析

1）顆粒力鏈分析

邊坡模型在降雨過(guò)程中力鏈的演變過(guò)程如圖10所示。從圖10可以看出，隨著降雨的持續(xù)進(jìn)行，邊坡飽和土體區(qū)域的力鏈分布變得比較稀疏，這說(shuō)明此部分土體顆粒間的接觸力較小，表現(xiàn)出極不穩(wěn)定的狀態(tài)。在上部流體和顆粒的作用下，力鏈在非飽和土體下部區(qū)域表現(xiàn)得比較密集和粗大，說(shuō)明邊坡內(nèi)部顆粒間的接觸力較大，此區(qū)域土體表現(xiàn)得比較穩(wěn)定。以此同時(shí)，力鏈沿坡體向前發(fā)生較大的延伸，斜坡變緩，坡頂及斜坡上力鏈減少，變得稀疏，說(shuō)明飽和區(qū)顆粒發(fā)生了較大移動(dòng)。當(dāng)模型計(jì)算至100 000時(shí)步后，非飽和區(qū)土體顆粒的力鏈基本沒有大的變化，如圖10（d）所示。但處于飽和區(qū)的坡面上顆粒進(jìn)一步減少，變得稀疏，說(shuō)明坡面上的顆粒繼續(xù)向下移動(dòng)，造成坡腳處顆粒堆積增多，致使坡腳處和地面上的力鏈變粗以及變得密集。

圖10 降雨作用下邊坡土體顆粒力鏈演變過(guò)程Fig.10 Evolution process of soil particle force chain in slope under rainfall

2）顆粒配位數(shù)變化分析

降雨過(guò)程中邊坡土體顆粒的配位數(shù)變化云圖如圖11所示。從圖11中可以看出，隨著降雨的持續(xù)作用，邊坡土體顆粒的配位數(shù)逐漸變小，其中位于飽和區(qū)域坡頂和坡面上土體顆粒的配位數(shù)降幅較大，而坡體內(nèi)部非飽和區(qū)域土體顆粒的配位數(shù)降幅較小，說(shuō)明在雨水作用的區(qū)域，顆粒間接觸不良，結(jié)構(gòu)穩(wěn)定性變差，在非飽和區(qū)的土體內(nèi)部，顆粒間接觸良好，該區(qū)域較為穩(wěn)定。當(dāng)模型計(jì)算到60 000時(shí)步后（圖11（c）、（d）），飽和區(qū)的配位數(shù)均變得較小，尤其是斜坡頂部，而坡體內(nèi)部和斜坡下部配位數(shù)較大，說(shuō)明上部顆粒的接觸間較差，在流體作用下，顆粒逐漸向下移動(dòng)至坡腳處堆積，坡腳處變得密實(shí)，穩(wěn)定性變好。

圖11 降雨作用下邊坡土體顆粒配位數(shù)變化過(guò)程Fig.11 Change process of soil particle coordination number in slope under rainfall

3）顆粒孔隙率變化分析

降雨過(guò)程中邊坡模型孔隙率變化計(jì)算結(jié)果見圖12。從圖12可以看出，隨著降雨持續(xù)進(jìn)行，整個(gè)邊坡土體的孔隙率都逐漸變大，其中，邊坡飽和區(qū)部分（主要是坡頂位置）的孔隙率增加較大，非飽和區(qū)部分的孔隙率增加幅度較小。從局部上看，坡腳處土體的孔隙率在降雨過(guò)程中增加的幅度較小，當(dāng)計(jì)算時(shí)步為60 000步以后，其孔隙率基本不再增加（最后增加至0.44），如圖12（c）所示，說(shuō)明在雨水作用下，坡頂處顆?；瑒?dòng)，致使顆粒在跛腳處堆積，導(dǎo)致坡腳處土體孔隙率變化不大。而斜坡頂上部土體的孔隙率增加尤為明顯，當(dāng)計(jì)算時(shí)步為100 000步時(shí)，其孔隙率增加至0.8，如圖12（d）所示，這說(shuō)明在雨水的作用下坡頂土體顆粒發(fā)生移動(dòng)，導(dǎo)致該區(qū)域土體顆粒變得疏松。

圖12 降雨作用下邊坡土體孔隙率變化過(guò)程Fig.12 Change process of porosity of slope soil under rainfall

2.4 討論

煤系土邊坡由于其易風(fēng)化，遇水容易軟化和崩解，在計(jì)算模擬時(shí)將其視為非連續(xù)介質(zhì)，采用離散元理論方法進(jìn)行計(jì)算分析是比較合理的。本文通過(guò)采用離散元方法模擬的計(jì)算結(jié)果與室外降雨試驗(yàn)結(jié)果進(jìn)行比較發(fā)現(xiàn)，兩種研究方法得出的煤系土邊坡的破壞模式都是雨水引起的沖刷破壞，邊坡滑動(dòng)面是近似的直線段，而且模型邊坡雨水沖刷的范圍和采用DEM-CFD耦合方法計(jì)算預(yù)測(cè)的滑動(dòng)面位置非常接近，這說(shuō)明本文采用的DEM-CFD流固耦合方法模擬煤系土邊坡的降雨失穩(wěn)破壞是合理的，這也驗(yàn)證了采用DEM-CFD流固耦合方法模擬分析巖土工程中的離散介質(zhì)運(yùn)動(dòng)規(guī)律是可行的[31]。

力鏈?zhǔn)欠从愁w粒運(yùn)動(dòng)宏觀力學(xué)性能的重要指標(biāo)，力鏈的粗細(xì)、位置變化、密集程度是顆粒受力及運(yùn)動(dòng)變化直觀的顯現(xiàn)。在降雨過(guò)程中，力鏈沿坡體向前產(chǎn)生較大的延伸，斜坡變緩，坡頂及斜坡上力鏈減少，且變得稀疏，說(shuō)明飽和區(qū)顆粒極不穩(wěn)定，整體發(fā)生了較大移動(dòng)。而在上部流體和顆粒的作用下，力鏈在下部區(qū)域表現(xiàn)得比較密集和粗大，說(shuō)明邊坡內(nèi)部接觸力較大，此區(qū)域表現(xiàn)得比較穩(wěn)定。在此次模擬邊坡降雨試驗(yàn)中，通過(guò)力鏈變化的分析能夠很好地解釋雨水作用下邊坡的破壞演變過(guò)程，對(duì)煤系土邊坡的破壞機(jī)理有了更深入的理解。

3 結(jié) 論

1）室外降雨試驗(yàn)中，煤系土邊坡坡率為1.0∶1.5時(shí)，其破壞的主要模式是雨水沖刷，在坡面形成深淺不一的沖蝕溝，最大沖蝕溝達(dá)0.36 m，發(fā)生在靠近坡中偏上的位置。本文采用DEM-CFD流固耦合方法模擬的煤系土邊坡破壞形式與室外降雨試驗(yàn)結(jié)果比較一致，邊坡滑動(dòng)面預(yù)測(cè)結(jié)果為近似的直線段，與實(shí)際模型邊坡雨水沖刷的范圍非常接近，這說(shuō)明采用DEM-CFD流固耦合方法對(duì)煤系土邊坡的穩(wěn)定性進(jìn)行分析是可行的。

2）邊坡土體顆粒的細(xì)觀參數(shù)，如力鏈、配位數(shù)以及孔隙率，在降雨過(guò)程中都會(huì)隨之發(fā)生變化，這些細(xì)觀參數(shù)與邊坡土體的宏觀力學(xué)表現(xiàn)直接關(guān)聯(lián)。因此，通過(guò)顆粒的細(xì)觀參數(shù)變化分析，可以很好地解釋雨水作用下煤系土邊坡的破壞演變規(guī)律，從微細(xì)觀角度分析非連續(xù)介質(zhì)煤系土邊坡的破壞機(jī)理是十分必要的。

3）本文研究結(jié)果表明，采用離散單元法和計(jì)算流體動(dòng)力學(xué)方法相結(jié)合模擬分析離散介質(zhì)邊坡的穩(wěn)定性是可行的，研究成果不僅為該區(qū)域煤系土邊坡的防護(hù)設(shè)計(jì)與施工提供理論依據(jù)，而且為從微細(xì)觀角度更好地分析離散介質(zhì)巖土工程的宏觀力學(xué)規(guī)律提供了一種新的思路。