基于MIDAS的尾礦邊坡三維流固耦合穩(wěn)定性分析

李良振，張延軍，秦勝伍，鄧浩，倪金

1.吉林大學 建設工程學院，長春 130026；2.中國地質(zhì)調(diào)查局 沈陽地質(zhì)調(diào)查中心， 沈陽 110034

0 引言

礦產(chǎn)資源作為現(xiàn)代工業(yè)的一大支柱，在國民經(jīng)濟的發(fā)展中具有重要作用。一般礦石開采之后，提取出目前對人類有用的成分，將暫時無法提取的礦物質(zhì)以及廢渣選取合適位置進行集中排放，從而形成尾礦壩[1]。尾礦壩是一種特殊的工業(yè)建筑物，其安全運行對礦山的效益、周圍人民的生命財產(chǎn)安全以及環(huán)境具有重要的影響。許多尾礦庫處于非常險要的地理位置，如江河湖泊、重大交通設施、密集居民區(qū)和耕地的上游。一旦壩體發(fā)生破壞就會造成重大的財產(chǎn)損失和人員傷亡[2]。因此，研究尾礦壩邊坡穩(wěn)定性對于生命安全和企業(yè)效益具有重要意義。

在滲流-應力耦合及邊坡穩(wěn)定性方面，國內(nèi)外大量學者形成了一套完備的計算理論。在滲流方面，Darcy通過砂桶實驗，提出了一定流速條件下的達西定律，從而奠定了層流滲流的理論基礎[3-4]。Perttson提出邊坡穩(wěn)定性計算方法-圓弧法，在 Taylor和 Fellenius的改進下形成瑞典圓弧法[5-7]。此后，相關學者根據(jù)破壞面形狀以及劃分條塊后條間力處理方式的不同，相繼建立了 Bishop 法，Janbu 法、morgenstern-price 法等[8]。在流固耦合方面，Terzaghi 提出了有效應力原理，得到了總應力、有效應力以及孔隙水壓力之間的關系，從而提出了一維固結模型[9-10]。Zienkiewicz 和 shiomi 在線性固結理論的基礎上考慮了非線性問題，改進了 Biot 公式，使得公式可以解決非線性問題[11-13]。李錫肈等考慮到飽和土固結效應的土體和其結構相互作用及相互影響的問題[14-15]，其研究使得土體固結在飽和土和非飽和土方面都有了相應的本構模型。為后續(xù)的直接計算和數(shù)值模擬提供了相對完備的理論基礎和數(shù)學模型。

但是，國內(nèi)外對于尾礦壩的滲透穩(wěn)定性分析相對較少，且一般只分析滲透情況，對于壩體邊坡穩(wěn)定性只給出定性判定，缺乏定量分析。因此，筆者在滲流分析的基礎上，結合吉林集安石墨尾礦壩實際情況和庫區(qū)工程地質(zhì)條件，采用有限元方法對該尾礦壩的流固耦合穩(wěn)定性進行綜合分析。

1 尾礦壩三維數(shù)值模型建立

本次模型建立基于谷歌大比例尺高清地形云圖，通過區(qū)域等高線獲取22 500個地形控制點，使用曲面擬合的方法重建區(qū)域內(nèi)地面。模型高152 m，長3 250 m,寬1 410 m(圖1)。初期壩采用透水堆石壩，壩底標高為505 m，壩頂標高為536 m，壩高31 m。后期堆積壩最終標高648 m，壩高120 m,坡度31.2°?？値烊葸_1 991.3萬m3。

圖1 尾礦壩三維網(wǎng)格模型Fig.1 3D mesh model of tailings dam

本次模型網(wǎng)格共生成的網(wǎng)格節(jié)點121 859個，生成的網(wǎng)格單元181 467個，保證了計算精度要求。由于該地區(qū)有較多石墨尾礦壩，且尾礦顆粒成分及組成大致相同。因此通過測量其他石墨尾礦壩的尾礦顆粒及初期壩物理力學參數(shù)，得到本模型的物理力學參數(shù)，如表1所示。

表1 尾礦壩體各部分物理力學參數(shù)

2 尾礦壩滲流分析

2.1 基本假設

(1)假設在滲流過程中，固體顆粒不發(fā)生運移，滲流僅僅是水的位移。

(2)假設滲流過程中，滲流速度在一定區(qū)間內(nèi)，即雷諾數(shù)(Re)上限在1～10之間，水流為層流，滲流滿足達西定律。

(3)假設滲透系數(shù)為常數(shù)，忽略溫度的變化

(1)

式(1)中，ρ為液體的密度，g為地球重力加速度，μ為水的動力粘滯系數(shù)，k為滲透率。

2.2 計算原理

根據(jù)Darcy定律,飽和砂體滲流方程：

(2)

式(2)中，K為滲透系數(shù)，Q為總流量，A為過水面積。H1,H2為水頭，L為滲透路徑長度。

在三維介質(zhì)中，飽和滲流微分方程為：

(3)

式(3)中，Kx,Ky,Kz為三個方向的滲透系數(shù)，為水頭。

在三維空間坐標系中，則有速度分量

(4)

Gd=i×γw

(5)

式(5)中，i為水力坡度，γw為水的容重，Gd為動水壓力。

通過水力比降，可以判定滲透力分布特點以及滲透力的大小。

2.3 計算結果

尾礦壩運行過程中，壩頂水位隨降水量的變化而變化，主要有干旱水位、正常水位和洪水位3種水位。水位變化影響壩頂干灘長度(壩坡頂?shù)剿痪€的距離)。根據(jù)測量多個臨近石墨礦尾礦壩干灘長度隨水位的變化情況，得到干灘長度隨水位變化情況如表2所示。

表2 干灘長度隨洪水位變化

考慮到洪水位時為最危險水位，如果壩體在洪水位時是安全的，則干旱水位和正常水位時一定安全。所以，只需要計算洪水時的滲流情況，并設置安全措施以保證壩體邊坡穩(wěn)定性。

通過計算，得到孔隙水壓力分布如圖2所示。孔隙水壓力的零勢面即為浸潤線位置。浸潤線切過尾礦壩邊坡，說明有水溢出。溢出高度41 m，溢出區(qū)域長度115±3 m。邊坡溢出水會沖刷尾礦壩邊坡坡面，對于壩坡穩(wěn)定性極為不利，如果不設置排滲措施，尾礦壩邊坡一定是危險的。

圖2 尾礦壩浸潤線分布Fig.2 Distribution of infiltration line of tailings dam

圖3 尾礦壩應力坡降分布Fig.3 Distribution of stress gradient of tailings dam

應力坡降范圍是0～1.247 4 m2·s-2。水力坡降等值區(qū)呈現(xiàn)帶狀分布，應力坡降最大值集中在水頭邊界處和堆積壩坡腳處，坡腳處為最大值位置(圖3)。根據(jù)滲透力，可以判定，可能發(fā)生破壞的部位在水頭邊界處和堆積壩坡腳處，其中堆積壩坡腳處破壞將最嚴重。由于上游水頭邊界處破壞為尾細砂破壞，而且破壞位置距離壩頂達到125 m，為安全地帶。而坡腳處由于滲透力產(chǎn)生破壞則會對尾礦壩邊坡產(chǎn)生非常不利影響。

為保證壩體安全，設置排滲體用以改變滲流路徑和滲流量，從而提高壩體穩(wěn)定性。排滲體為砂墊層，滲透性比尾細砂大。每個砂墊層呈水平狀，降低了施工的難度(圖4)。砂墊層厚1 m，寬約100 m，根據(jù)規(guī)范要求，結合模型基巖形狀和尾礦壩形態(tài)，一共設置11塊排滲體。排滲體物理力學參數(shù)如表3所示。

圖4 壩體排滲體布置Fig.4 Arrangement of seepage body of dam body

名稱天然含水率/%天然密度/g·cm-3孔隙比比重滲透系數(shù)/cm·s-1內(nèi)摩擦角/°內(nèi)聚力/kPa承載力特征值/kPa排滲體131.830.5462.491.0×10-22520110

設置排滲體之后，浸潤線埋深明顯加大，橫向埋深增加至101±1.5 m，縱向埋深增加至35±0.5 m，尾礦壩安全性極大提高(圖5)。此時，水由初期壩溢出，溢出高度5±0.5 m，符合規(guī)范要求。初步判斷尾礦壩邊坡安全。

應力坡降范圍在0～1.100。水力坡降等值區(qū)呈現(xiàn)出帶狀分布(圖6)。應力坡降最大值集中在上游水頭邊界、中間砂墊層以及堆積壩坡腳處，其中，上游水頭邊界坡降最大，為1.100 m2·s-2。但是由于距離壩頂位置129 m，距離較遠，即使發(fā)生滲透破壞，對于尾礦壩邊坡的穩(wěn)定性影響不大。

圖5 尾礦壩浸潤線分布Fig.5 Distribution of infiltration line of tailings dam

圖6 尾礦壩應力坡降分布Fig.6 Distribution of stress gradient of tailings dam

3 滲流-應力-邊坡耦合

3.1 本構方程

忽略溫度場和化學場對應力場的影響，只考慮滲流和應力之間的耦合，方程為：

(6)

式(6)中，p為孔隙水壓力，t為時間，s為飽和度，v為流體容量，α為比奧系數(shù)，ε為體積應變。

3.2 強度折減法

考慮自重條件下的邊坡穩(wěn)定性，使用強度折減法進行計算。即將壩體材料的強度參數(shù)，按照一定比例系數(shù)，進行不斷折減，直至壩體達到極限破壞狀態(tài)。與此同時，壩體破壞部位會形成等效塑性應變區(qū)，塑性應變貫通的位置即為滑動面(帶)位置。

(7)

式(7)中，c′為折減后的內(nèi)聚力，φ′為折減后的內(nèi)摩擦角，c為材料初始的內(nèi)聚力，φ為材料初始的內(nèi)摩擦角。R為強度折減系數(shù)，即最終的安全系數(shù)是臨界破壞時的強度折減系數(shù)。

由于排滲體強度參數(shù)與尾細砂強度參數(shù)基本一致。為了計算簡便，將滲流計算結果直接附到原壩體進行計算。

3.3 計算結果

圖7 壩體最小主應力Fig.7 Minimum principal stress of dam body

最小主應力和最大主應力均為負值，說明均為壓應力，沒有拉應力(圖7，圖8)。最大剪應力在深度最大處達到最大值，往上呈現(xiàn)層狀逐級遞減，判斷剪應力主要受自重影響比較大(圖9)。

圖8 壩體最大主應力Fig.8 Maximum principal stress of dam body

圖9 壩體最大剪應力Fig.9 Maximum shear stress of dam body

根據(jù)塑性應變區(qū)分布，可以判定，壩頂處是最大值處，且壩坡上緣與坡腳處塑性應變區(qū)連通，形成閉合。三維邊坡滑動面與二維邊坡滑動面不同，為槽形滑動面。堆積區(qū)后部的塑性變形與周圍山體邊界大致呈平行條帶狀分布(圖10)。判斷后緣的變形是由于在相同的重力條件下，山體的強度高變

圖10 壩體塑性應變區(qū)分布Fig.10 Distribution of plastic strain zone in dam

形量小，堆積砂體較為軟弱變形量大，所以引起變形不協(xié)調(diào)，從而產(chǎn)生與山體接觸邊緣大致平行分布的塑性變形帶。綜上所述，塑性變形帶由初期壩前緣延伸至堆積區(qū)后緣，一旦發(fā)生破壞，滑體規(guī)模是十分巨大的。由于初期壩強度很高，所以塑性變形小，塑性變形帶被初期壩截斷。現(xiàn)設計的尾礦壩安全系數(shù)為1.536 4，遠遠大于規(guī)范要求，不會發(fā)生破壞。

4 結論

(1)尾礦壩邊坡中設置的排滲體主要對壩體浸潤線有直接影響，設置排滲體之后，浸潤線埋深顯著增加，達到穩(wěn)定狀態(tài)。

(2)排滲體可以改變滲流途徑，從而降低坡面以及坡腳的水力坡降，降低了滲透破壞的可能性。

(3)設置排滲體之后，滲流情況發(fā)生了變化。計算得到滑動帶是與山谷近乎平行的槽狀，由初期壩前緣延伸至堆積區(qū)后緣。尾礦壩邊坡穩(wěn)定性系數(shù)達到1.536 4，符合規(guī)范安全要求。