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

(1.南昌工程學院 江西省水利土木工程基礎設施安全重點實驗室，江西 南昌 330099；2.南昌工程學院 土木與建筑工程學院，江西 南昌 330099)

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

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

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

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

1 DEM-CFD流固耦合計算基本原理

1.1 流體-顆粒相互作用

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

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

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

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

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

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

式（3）中的經驗系數(shù) χ為：

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

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

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

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

1.2 流體-流體相互作用

1）連續(xù)方程

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

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

2）動量方程

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

式中，p為流體壓力。

1.3 顆粒運動控制方程

顆粒的運動方程是在牛頓第二定律的基礎上推導出來的，為了考慮流體對顆粒的作用，在控制方程中引入顆粒與流體間的相互作用力ffluid，得到：

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

1.4 CFD-DEM耦合流程與程序實現(xiàn)

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

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

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

2.1 煤系土邊坡室外模型試驗

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

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

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

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

2.2.1 模型參數(shù)標定

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

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

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

圖4 三軸試驗細觀參數(shù)標定Fig.4 Meso parameters calibrated by simulated triaxial test

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

2.2.2 邊坡計算模型建立

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

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

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

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

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

為使流體流動方向與實際一致，Z軸的流速與X軸的流速之比設為1.0∶1.5，流體流動的方向與室外降雨試驗坡體的徑流一致。

2.3 計算結果分析

2.3.1 煤系土邊坡顆粒運動的宏觀分析

1）邊坡土顆粒運動軌跡計算結果分析

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

圖6 邊坡側面顆粒運動軌跡示意圖Fig.6 Schematic diagram of the overall particle movement of the slope

2）邊坡土體應力變化計算結果分析

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

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

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

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

圖9 計算至100 000時步邊坡各方向應力變化（中截面）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)破壞的細觀機理分析

1）顆粒力鏈分析

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

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

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

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

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

3）顆?？紫堵首兓治?/p>

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

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

2.4 討論

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

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

3 結 論

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

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

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