基于耦合模型的尾礦壩穩(wěn)定性對比分析

嚴(yán) 瓊 吳順川 李 龍 王強(qiáng)志

( 北京科技大學(xué)土木與環(huán)境工程學(xué)院，北京100083)

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

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

1 連續(xù)－離散耦合計(jì)算方法原理及應(yīng)用

1.1 連續(xù)－離散耦合原理概述

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

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

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

1.2 連續(xù)－離散耦合作用理論

節(jié)段端點(diǎn)坐標(biāo)和長度分別為x0，x1和l，顆粒作用力為P，顆粒質(zhì)心為xp。作用于節(jié)段端點(diǎn)的力為F0和F1:

其中，t、n 分別為節(jié)點(diǎn)坐標(biāo)系切向和法向的分量;m0，m1，m2和m3分別為與單位向量對應(yīng)的變量。

節(jié)段端點(diǎn)坐標(biāo)和長度分別由x0，x1和l 表示，顆粒坐標(biāo)由xp表示，端點(diǎn)的速度由v0和v1表示。顆粒速度滿足公式

在此，ξ 是局部坐標(biāo)，為

圖2 為映射顆粒與節(jié)段節(jié)點(diǎn)反作用力與相應(yīng)速率示意圖。

圖2 映射顆粒與節(jié)段節(jié)點(diǎn)反作用力和相應(yīng)速率示意Fig.2 Schematic diagram of the mapping particles，reaction force of segmented point and its associated velocity

1.3 耦合算法及細(xì)觀力學(xué)分析在尾礦壩中的應(yīng)用

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

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

2 失穩(wěn)尾礦壩工程案例簡介

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

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

3 尾礦壩耦合模型構(gòu)建與方案比選

3.1 應(yīng)力－滲流耦合模型

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

圖3 襄汾尾礦壩失穩(wěn)飽和狀態(tài)幾何剖面Fig.3 Geometric profile of instability ＆saturation state of Xiangfen tailings dam

圖4 節(jié)點(diǎn)孔隙水壓力監(jiān)測結(jié)果Fig.4 Monitoring result of the grid-point pore pressure

圖5 襄汾尾礦壩應(yīng)力－滲流計(jì)算結(jié)果Fig.5 The strees-seepage computed results of Xiangfen tailings dam

表1 計(jì)算模型參數(shù)Table 1 Parameters of the model

3.2 連續(xù)－離散耦合模型

3.2.1 細(xì)觀參數(shù)選取

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

表2 PFC 模型主要計(jì)算參數(shù)Table 2 Main calculation parameters of PFC model

圖6 雙軸試驗(yàn)應(yīng)力－應(yīng)變曲線Fig.6 Stress-strain curve of biaxial tests

3.2.2 耦合方案比選

耦合模型的主旨是通過耦合域的離散模型，從細(xì)觀力學(xué)的角度，分析宏觀變形特征。不同的耦合方案，選取不同的耦合區(qū)域，耦合介質(zhì)表現(xiàn)出的接觸力和位移狀況不同，但不會(huì)從整體上影響邊坡的穩(wěn)定性分析結(jié)果。

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

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

3.2.3 耦合模型建立

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

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

連續(xù)－離散耦合方法中連續(xù)域作為模型主體，在耦合初始化計(jì)算時(shí)，耦合域的離散元顆粒試樣也需要賦以一定的應(yīng)力狀態(tài)，以契合和反饋連續(xù)域節(jié)點(diǎn)經(jīng)虛擬做功傳遞的能量。離散元的初始應(yīng)力場在試樣生成后，采用雙軸應(yīng)力加載方式獲得。

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

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

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

應(yīng)力－滲流耦合方式作為較成熟的尾礦壩數(shù)值分析方法，可從應(yīng)力和水力結(jié)合的方式，展示水力坡降的動(dòng)態(tài)效果影響，但無法從細(xì)觀力學(xué)角度解釋巖土顆粒位移形變的成因，而連續(xù)－離散的耦合方法正彌補(bǔ)了這一缺陷。

表3 尾礦壩計(jì)算狀態(tài)和時(shí)步對照表Table 3 The computational state and time-step contrast of the tailings dam

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

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

4.2 連續(xù)－離散耦合模型的失穩(wěn)表征協(xié)同分析

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

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

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

5 結(jié) 論

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

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

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

總之，計(jì)算結(jié)果分析表明，采用的基于耦合模型的尾礦壩穩(wěn)定性分析，可有效揭示山西襄汾尾礦壩潰壩的宏細(xì)觀變形規(guī)律，具有一定的科學(xué)應(yīng)用價(jià)值。

［1］ 馬萃林，朱 明，王建華. 有限差分法FLAC 在邊坡穩(wěn)定性分析中的應(yīng)用［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，F(xiàn)elippa 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］ 張 鐸，劉 洋，吳順川，等. 基于離散-連續(xù)耦合的尾礦壩邊坡破壞機(jī)理分析［J］.巖土工程學(xué)報(bào)，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］ 張延軍.邊坡滲流耦合變形分析方法的研究及其應(yīng)用［J］.吉林大學(xué)學(xué)報(bào):地球科學(xué)版，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］ 周冬磊，王連國，黃繼輝，等. 裂隙巖體應(yīng)力滲流耦合規(guī)律及對底板隔水性能研究［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］ 繆海波，殷坤龍，郭付三，等. 金屬礦山尾礦壩滲流場模擬及穩(wěn)定性數(shù)值分析［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］ 尹光志，張千貴，魏作安，等.尾礦細(xì)微觀力學(xué)與形變觀測試驗(yàn)裝置的研制與應(yīng)用［J］. 巖石力學(xué)與工程學(xué)報(bào)，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］ 董成會(huì).尾礦堆積壩潰壩分析及數(shù)值模擬研究［D］.南昌:南昌大學(xué)，2012.

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