基于耦合模型的尾礦壩穩定性對比分析

嚴 瓊 吳順川 李 龍 王強志

( 北京科技大學土木與環境工程學院，北京100083)

礦山資源的開采給社會帶來豐富原材料的同時也給社會帶來了一定的隱患和危害，尾礦庫就是其一。近年來尾礦庫潰壩事故時有發生，破壞性巨大，造成不可挽回的生命和財產損失。因此，尾礦庫穩定性的研究愈發重要。現代科學技術的發展及計算機的應用已使人們在尾礦壩穩定性研究中有一些進展，然而，這些研究多采用有限元或有限差分的連續體模型算法［1］，而尾礦庫宏觀失穩破壞都是細觀顆粒累積變形發展的結果，故已往的算法已不能滿足分析的需要。

連續－離散耦合的算法及巖土體宏細觀受力和變形相結合的分析自20 世紀80 年代初由Park 和Felippa 等［2］提出以來，已得到長足發展。目前，連續－離散耦合算法的開發，同時結合尾礦壩各影響因素的細觀力學分析已成為工程界的熱點研究課題。

1 連續－離散耦合計算方法原理及應用

1.1 連續－離散耦合原理概述

連續－離散耦合是巖土體數值模擬中計算方法的耦合，本離散、連續域分別采用顆粒單元法和有限差分法進行計算。顆粒單元( 以PFC 程序實現) 作為耦合離散域，被嵌入的有限差分單元( 以FLAC 程序實現) 在建模時需要預留與離散域一致的空間，離散域以單層顆粒圍繞為邊界。

圖1 為耦合原理的簡單圖示，圖1( a) 為預留離散域空間的連續域網格，箭頭表示連續域節點受到來自離散域的作用力，圖1( b) 為耦合離散域范圍，周邊深灰色顆粒為單層顆粒組成的耦合邊界，黑色線段和黑點組成的沿深灰色顆粒走向的方框為節段序列，黑點所在的顆粒為控制顆粒，箭頭表示顆粒的瞬時速度［5］。

圖1 耦合連續域與離散域示意Fig.1 Schematic diagram of coupling continuum and discrete domains

1.2 連續－離散耦合作用理論

節段端點坐標和長度分別為x0，x1和l，顆粒作用力為P，顆粒質心為xp。作用于節段端點的力為F0和F1:

其中，t、n 分別為節點坐標系切向和法向的分量;m0，m1，m2和m3分別為與單位向量對應的變量。

節段端點坐標和長度分別由x0，x1和l 表示，顆粒坐標由xp表示，端點的速度由v0和v1表示。顆粒速度滿足公式

在此，ξ 是局部坐標，為

圖2 為映射顆粒與節段節點反作用力與相應速率示意圖。

圖2 映射顆粒與節段節點反作用力和相應速率示意Fig.2 Schematic diagram of the mapping particles，reaction force of segmented point and its associated velocity

1.3 耦合算法及細觀力學分析在尾礦壩中的應用

張鐸等［6］采用離散－連續耦合算法模擬了某尾礦壩邊坡在尾礦充填前后潛在滑移帶附近的宏細觀力學特征，驗證了尾礦壩在荷載增加中的漸進破壞的細觀力學解釋方法。張延軍［6］采用了不同含水量對抗剪強度影響的模型，利用滲流耦合變形的彈塑性有限元法，證實由于滲流作用影響，邊坡變形加大，最后形成潛在滑動面。周冬磊等［6］采用巖石力學伺服試驗系統，對煤礦底板的巖樣進行全應力－應變滲透試驗，得到巖樣裂隙在閉合、萌生、擴展過程中的滲透特性。繆海波等［6］應用滲流分析模塊，模擬在初期壩不同透水情況、同長度條件下尾礦壩的滲流場變化規律，獲得相應工況下的危險滑動面與危險系數。尹光志等［6］利用尾礦細微觀力學與形變觀測試驗裝置，對定量研究尾礦壩體細觀結構特征和力學特性及探索水引起的尾礦壩災變機制具有重要意義。

目前，離散－連續耦合的分析已經從方法上證實了其有效性，也在邊坡穩定性研究中有了一定的嘗試，但在考慮尾礦壩體水力作用引起的失穩破壞研究中，還未有深入的探索。本研究以應力－滲流耦合分析為基礎，結合連續－離散耦合建模理論，對失穩破壞的尾礦壩建立數值分析耦合模型，對比動水作用和靜水作用對尾礦壩的影響，對尾礦壩失穩的成因和機理從細觀力學角度進行了合理解釋。該方法彌補了現有尾礦壩長期采用連續元進行穩定性分析和連續－離散耦合計算中水力影響因素欠缺的不足，為邊坡及尾礦壩穩定性分析方法進行了又一補充。

2 失穩尾礦壩工程案例簡介

發生事故的山西襄汾尾礦庫建于20 世紀80 年代，1992 年停止使用。尾礦庫被封閉后，曾采取碎石填平、黃土覆蓋壩頂、植樹綠化下游壩坡、庫區上方建設排洪明渠等閉庫處理措施。新塔礦業公司通過拍賣購買了鐵礦產權，本應重新修建新的尾礦庫，但礦方卻在舊庫上挖庫排放尾礦，造成尾礦庫大面積液化，壩體失穩，并引發了這起特別重大的潰壩事故。

原廢棄尾礦庫建有二級壩，初期壩壩底高程約900 m，壩頂高程934 m，壩高約34 m，初期壩長約100 m，儲礦量26.8 萬m3;次級壩壩底高程934 m，壩頂高程950 m，壩高16 m，壩長約120 m，儲礦量約10.3 萬m3，尾礦庫總庫容約36.8 萬m3。尾礦庫潰壩后，尾礦庫次級壩后緣120 m 范圍內、深10 ～15 m的尾礦在平面上呈扇形依一定的坡度隨水流一起沖出庫區泄向下游，其泄出總量達19 萬m3，將尾礦庫下游的大量建筑物掩埋與破壞。

3 尾礦壩耦合模型構建與方案比選

3.1 應力－滲流耦合模型

浸潤線作為尾礦壩的生命線，其高度與壩體穩定性休戚相關。從圖3 尾礦壩失穩飽和狀態剖面圖可以看到，浸潤線過高，干灘長度過短，壩體已具高勢能飽和狀態，是導致潰壩的直接原因。根據該尾礦壩潰壩資料，結合出水點水位高度，采用Flac 軟件構建應力－滲流計算模型。根據該壩體運營模式，模型采取底層壩體先期平衡再進行上層堆載加壓的方式。土體本構模型選取摩爾－庫侖模型，邊界條件包括應力和水力邊界，應力場采用巖土體自重應力場，計算時步通過監測浸潤線附近某點的孔隙水壓力獲得( 見圖4，總時步約44 000 步) 。根據模型的幾何剖面特征( 如圖3) 和文獻［10］的土層參數，模型計算參數見表1，應力－滲流模型和滲流、浸潤線分布圖如圖5。

圖3 襄汾尾礦壩失穩飽和狀態幾何剖面Fig.3 Geometric profile of instability ＆saturation state of Xiangfen tailings dam

圖4 節點孔隙水壓力監測結果Fig.4 Monitoring result of the grid-point pore pressure

圖5 襄汾尾礦壩應力－滲流計算結果Fig.5 The strees-seepage computed results of Xiangfen tailings dam

表1 計算模型參數Table 1 Parameters of the model

3.2 連續－離散耦合模型

3.2.1 細觀參數選取

耦合模型的計算參數包括宏觀參數和細觀參數，宏觀參數取值見表1 的沉積尾砂數據。而細觀參數的選取是耦合計算中較為關鍵的步驟，不但要保證與連續域變形一致，也要確保能有效表征細觀受力特征。根據宏觀參數和材料特征，用實驗室模型進行宏細觀參數的反復調試，本研究分別采用離散元和有限差分模擬雙軸壓縮試驗，當二者數值試驗的結果接近時，可近似認為宏觀與細觀參數相互對應。離散元模型采用平行黏結模型，細觀參數取值見表2，雙軸試驗的結果見圖6。

表2 PFC 模型主要計算參數Table 2 Main calculation parameters of PFC model

圖6 雙軸試驗應力－應變曲線Fig.6 Stress-strain curve of biaxial tests

3.2.2 耦合方案比選

耦合模型的主旨是通過耦合域的離散模型，從細觀力學的角度，分析宏觀變形特征。不同的耦合方案，選取不同的耦合區域，耦合介質表現出的接觸力和位移狀況不同，但不會從整體上影響邊坡的穩定性分析結果。

單就耦合域面積來說，越大越可能更多地顯示顆粒體特征，但由于顆粒離散元試樣從生成到耦合初始化，需要經過各種參數配比和調試，耦合域大小并不代表能獲得合理的宏細觀數據。此外，耦合區域所在位置的選取更是耦合方案的關鍵。對比圖7 的各個耦合域方案，結合圖8 壩體非耦合方式加載的位移變化，方案a 既能表現一級壩平衡后的特征，又能突出二級壩加載后的變形影響，可作為合理耦合域的選定范圍。

圖7 尾礦壩耦合模型不同耦合域方案Fig.7 Respective discrete domain of tailings dam coupling model

3.2.3 耦合模型建立

連續－離散耦合模型的建立包含2 部分:連續域建模形式與應力－滲流模型的應力計算部分基本一致，采用FLAC 計算;離散域作為連續域中的一部分，采用PFC 計算，既通過邊界顆粒與連續域進行數據交換，與連續域保持著受力和形變的一致性，也在該區域內顯示出巖土體獨有的細觀力學特征。

圖8 尾礦壩模型加載過程及位移矢量Fig.8 Loading procedure and displacement distribution graph of tailings dam model

連續－離散耦合方法中連續域作為模型主體，在耦合初始化計算時，耦合域的離散元顆粒試樣也需要賦以一定的應力狀態，以契合和反饋連續域節點經虛擬做功傳遞的能量。離散元的初始應力場在試樣生成后，采用雙軸應力加載方式獲得。

本案例作為含水的高勢能飽和土體壩，水力因素必須考慮。但由于目前連續－耦合計算方法中，水力作用的研究尚不太成熟，因此，在應力－滲流計算結果的基礎上，通過初步確定壩體浸潤線，其上下分別采用不含水及含水飽和土體的計算參數，其中飽和土體應用有效應力原理，考慮靜水浮力作用進行耦合計算。總體設定的計算時步為40 500 步，分設4 個監測時步對應壩體各計算狀態( 見表3) ，計算參數參照表1 和表2，耦合計算模型見圖9。

4 計算結果與分析

4.1 不同耦合模型的宏觀變形對比分析

應力－滲流耦合方式作為較成熟的尾礦壩數值分析方法，可從應力和水力結合的方式，展示水力坡降的動態效果影響，但無法從細觀力學角度解釋巖土顆粒位移形變的成因，而連續－離散的耦合方法正彌補了這一缺陷。

表3 尾礦壩計算狀態和時步對照表Table 3 The computational state and time-step contrast of the tailings dam

圖9 尾礦壩耦合計算中的模型Fig.9 The models in coupling calculation of the tailings dam

在本案例分析中，連續－離散耦合將動水作用簡化為靜水作用，從2 種耦合模型結果對比、位移發展和塑性變形效果方面，可對本案例的尾礦壩潰壩成因給予詳盡的分析。由圖10、圖11，2 種分析方法位移矢量及塑性變形圖表明:應力－滲流耦合模型計算平衡的塑性變形區域與連續－離散耦合模型基本一致，壩體在2 種計算模型下的宏觀變形結果基本相同，且與該尾礦壩潰壩實例一致( 如圖12) 。

4.2 連續－離散耦合模型的失穩表征協同分析

本研究中連續－離散耦合的數值計算方法采用有限差分FLAC 和離散元PFC 協同計算，可獲得對已失穩襄汾尾礦壩分析對象較為契合的結果。連續－離散耦合模型在飽和土體中根據有效應力原理，采用浮重度的計算參數，計算結果的位移已達失穩態勢。通過對比2 個模型的位移和塑性變形圖( 圖10、圖11) ，可以看出: 基于連續－耦合模型的靜水作用計算方法，可以有效表征含水尾礦壩體失穩的情況，雖位移矢量大小與應力－滲流模型有差別，但壩體失穩的發展趨勢較為一致，用于本例的宏細觀分析，可達到較好的分析效果。

圖10 尾礦壩位移矢量Fig.10 The displacement vector diagram of the tailings dam

圖11 尾礦壩塑性變形Fig.11 The plastic deformation diagram of the tailings dam

5 結 論

(1) 應力－滲流耦合方法作為成熟的尾礦壩穩定性研究方法，相比連續－離散耦合方法，在宏觀上能較合理地表征邊坡失穩態勢，也能為后者提供浸潤線和滲流方向等研究基礎數據。結果表明，二者在壩體位移的失穩發展和壩體塑性變形方面呈現較為一致的趨勢。

圖12 襄汾尾礦壩潰壩實例Fig.12 The break case of XiangFen tailings dam

(2) 根據有效應力原理和應力－滲流耦合的浸潤線計算結果，對壩體進行土層飽和與非飽和參數的劃分，引入浮重度參數，可對尾礦壩進行連續－離散耦合的計算分析，但所需的平衡時步比應力－滲流耦合的計算時步少，位移矢量大小也明顯偏小，但其位移矢量值已達失穩狀態。

總之，計算結果分析表明，采用的基于耦合模型的尾礦壩穩定性分析，可有效揭示山西襄汾尾礦壩潰壩的宏細觀變形規律，具有一定的科學應用價值。

［1］ 馬萃林，朱 明，王建華. 有限差分法FLAC 在邊坡穩定性分析中的應用［J］.中國礦山工程，2008，37(5) :19-22.

Ma Cuilin，Zhu Ming，Wang Jianhua.Application of finite difference method FLAC in analyzing slope stability［J］.China Mine Engineering，2008，37(5) :19-22.

［2］ Park K C，Felippa C A.Partitioned analysis of coupled systems［C］∥Belytschko T，Hughes T J R，et al. Computational Methods for Transient Analysis.Amsterdam:［s.n.］，1983:157-219.

［3］ Felippa C A，Park K C. Staggered transient analysis procedures for coupled-field mechanical systems: formulation［J］. Computer Methods in Applied Mechanics and Engineering，1980，24(1) :61-111.

［4］ Onate E，Rojek J.Combination of discrete element and finite element methods for dynamic analysis of geomechanics problems［J］. Computer Methods in Applied Mechanics and Engineering，2004，193:3087-3128.

［5］ Itasca Consulting Group Inc..PFC2D version 4.0 Verification Problems and Example Applications［M］. Minnesota: Itasca Consulting Group Inc.，2008.

［6］ 張 鐸，劉 洋，吳順川，等. 基于離散-連續耦合的尾礦壩邊坡破壞機理分析［J］.巖土工程學報，2014，36(8) :1473-1482.

Zhang Duo，Liu Yang，Wu Shunchuan，et al. Failure mechanism analysis of tailing dams based on coupled discrete and continuous method［J］. Chinese Journal of Geotechnical Engineering，2014，36(8) :1473-1482.

［7］ 張延軍.邊坡滲流耦合變形分析方法的研究及其應用［J］.吉林大學學報:地球科學版，2006，36(1) :103-107.

Zhang Yanjun.The study of the analysis method of seepage coupled deformation in slope and its application［J］.Journal of Jilin University:Earth Science Edition，2006，36(1) :103-107.

［8］ 周冬磊，王連國，黃繼輝，等. 裂隙巖體應力滲流耦合規律及對底板隔水性能研究［J］.金屬礦山，2011(11) :53-57.

Zhou Donglei，Wang Lianguo，Huang Jihui，et al.Coupled behavior of stress and permeability in fractured rock masses and its study of water isolation［J］.Metal Mine，2011(11) :53-57.

［9］ 繆海波，殷坤龍，郭付三，等. 金屬礦山尾礦壩滲流場模擬及穩定性數值分析［J］.金屬礦山，2010(3) :134-138.

Miao Haibo，Yin Kunlong，Guo Fusan，et al.Seepage field simulation and numerical analysis of stability of metal mine tailings Dam［J］.Metal Mine，2010(3) :134-138.

［10］ 尹光志，張千貴，魏作安，等.尾礦細微觀力學與形變觀測試驗裝置的研制與應用［J］. 巖石力學與工程學報，2011，30( 5) :926-934.

Yin Guangzhi，Zhang Qiangui，Wei Zuoan，et al. Development and application of observation testing apparatus for micromechanics and deformation of tailings［J］.Chinese Journal of Rock Mechanics and Engineering，2011，30(5) :926-934.

［11］ 董成會.尾礦堆積壩潰壩分析及數值模擬研究［D］.南昌:南昌大學，2012.

Dong Chenghui. Analysis of Dam-break and Numerical Simulation Study on Tailings Stacking Dam［D］.Nanchang:Nanchang University，2012.