王梓鑒，楊 俊

（九江市水利電力規(guī)劃設計院，江西 九江 332000）

0 引言

尾礦庫是由金屬或非金屬等礦山選礦的廢棄物所構成的，由于我國對礦產資源的需求與日俱增，尾礦庫的數(shù)量也在不斷增多。近年來，地震、洪水等原因導致我國尾礦庫潰壩事故頻繁發(fā)生，尾礦料的傾瀉可能導致下游居民房屋受到沖擊、農田遭受破壞，同時尾礦料中的一些重金屬離子將會造成水污染，嚴重影響當?shù)鼐用竦纳a生活。目前國內外學者針對尾礦壩的穩(wěn)定性問題開展了大量有益的研究工作，如龍明魁[1]等對元江尾礦壩進行靜力穩(wěn)定性分析，研究尾礦壩的位移、應力變化規(guī)律；陳燕[2]等開展了不同降雨強度下的尾礦壩穩(wěn)定性分析，尾礦壩穩(wěn)定性會隨著降雨強度的增加而下降；李良振[3]進行滲流-應力-邊坡耦合，得出尾礦邊坡穩(wěn)定性系數(shù)，為尾礦壩具體施工提供了技術支持。然而，尾礦材料是性質極其復雜的地質介質，在施工過程中人為或非人為因素作用下，其物質組成與內部結構會發(fā)生一定的變化。實際上尾礦壩一般處于飽和的疏松狀態(tài)，顆粒細小，并且比重較大，由于顆粒的組成、礦物的成分等因素的影響，它們看似穩(wěn)定，實際上對擾動特別敏感，在地震中很容易發(fā)生液化和破壞變形[4]，因此，研究我國尾礦壩的抗震穩(wěn)定性就顯得尤為緊迫和必要。

由于尾礦壩的動力失穩(wěn)機制復雜，目前多數(shù)的研究集中于將動力問題簡化為擬靜力問題來計算尾礦壩的安全系數(shù)，如李再光[5]使用擬靜力法分析了南京鋼鐵尾礦壩在地震作用下的穩(wěn)定性。當考慮簡化地震作用效應時，用擬靜力法簡單有效，但地震作用下的動力特性卻得不到真實地反映。動力有限元分析法將每一時刻的動應力施加到相應的靜應力上，并對尾礦壩進行靜力分布狀態(tài)分析判斷靜力穩(wěn)定性，作為壩體動力反應的基礎，通過輸入地震時稱曲線，比較適用邊界條件和巖土體結構的復雜動力問題，從而對尾礦壩進行地震過程中的變形分析和安全評價。另外使用結構可靠度分析方法可以使得穩(wěn)定分析結果更具說服力。

本文結合動力有限元模型和結構可靠度理論，使用GEOSTUDIO軟件，并將Seep/w、Quake/w、Slope/w三個模塊相耦合，對江西省德安尾礦壩工程進行滲流分析、地震響應分析以及尾礦壩的穩(wěn)定計算，最后采用蒙特卡羅法計算德安尾礦壩的失效概率，分析相關研究成果，對于德安尾礦壩的安全生產、生態(tài)環(huán)境的維護以及周邊居民的健康生命安全具有十分重要的意義。

1 尾礦壩動力穩(wěn)定性分析原理

在尾礦庫地震失穩(wěn)狀態(tài)下，尾礦土的應力應變關系表現(xiàn)為非線性關系，本文采用等效線性本構模型[6]。等效線性模型是黏彈性模型的一種類型，特點在于將不同應變幅值下的滯回特性用阻尼比隨應變變化，即將土體視為黏彈性體，反映土體滯后性和動剪應力－應變關系非線性的兩個基本特征[7]。等效線性模型計算見式（1）。

式中，AL為應力應變滯回環(huán)面積；λ為等效阻尼比；AT為滯回環(huán)拐點與x軸、坐標原點O所圍成的直角三角形面積。

如果震動頻率很低，滯回環(huán)不包括黏性阻尼，會導致阻尼比很小，則需要修正公式進行修正，見式（2）、式（3）。

式中，λv為黏性阻尼比；δi、δi+1為相鄰兩個周期的位移振幅。

動力分析的步驟與靜力方法大致相似，都要先把尾礦壩壩體離散為單元，然后對單元進行動力分析，動力有限元方程式為：

式中，[M ]為質量矩陣，可用集中質量法求得；[C]為單元阻尼矩陣；[K ]為單元剛度矩陣，通過常規(guī)有限元求得；a為加速度；u為速度；s為單元節(jié)點的位移；｛F（ t） ｝為單元節(jié)點上的動力荷載，可以表示為式（5）：

式中，｛Fb｝為單元結構上的體力；｛Fs｝為邊界上的面力；｛Fn｝是節(jié)點上集中力；｛Fg｝為地震荷載引起的地應力。采用Rayleigh阻尼假設單元阻尼矩陣，即阻尼矩陣是剛度矩陣和質量矩陣的線性組合[8]：

式中，α、β 為比例系數(shù)。另外 α=λγω0，β=λγ/ω0，λγ為由單元結構剪應變所決定的阻尼比，ω0為壩體與尾礦庫無阻尼系統(tǒng)的頻率。

目前，通常采用Wilson線性加速度法求解式（4），然后逐步積分，求出每一時刻節(jié)點上的位移｛st｝、速度｛ut｝和加速度｛at｝。

2 可靠度分析方法

可靠度是指規(guī)定的時間范圍內，工程在規(guī)定的條件下所能夠完成預期功能的概率，對評價工程整體具有一定的指導意義，而無法完成預期功能的概率就定義為失效概率[9]?？煽慷鹊姆治龇椒ㄓ泻芏啵鏙C法、一次矩陣法、蒙特卡羅法等，本文采用蒙特卡羅法（MCS）。MCS方法采用大量重復抽樣的基本思想，通過模擬實驗次數(shù)來獲取失效概率值，當設定的次數(shù)越高，所得到的失效概率值也就越精確。失效概率Pf=m/n，式中m為失效次數(shù)，n為實驗次數(shù)。使用MCS對尾礦壩土體輸入?yún)?shù)進行隨機抽樣后，可將輸入?yún)?shù)特性值分別賦給尾礦壩穩(wěn)定性有限元分析模型，再進行尾礦壩穩(wěn)定性有限元分析。然而大量的穩(wěn)定性分析計算繁瑣、效率低下，本文通過將MATLAB軟件編寫程序與GEOSTUDIO軟件實行對接，并在DOS環(huán)境下進行相關批處理計算，無需打開GEOSTUDIO軟件界面，有效地避免了重復操作有限元軟件，大大提高了計算效率。隨后通過批處理計算獲得相應的n個臨界安全系數(shù)FSmin和滑動面等信息來計算尾礦壩的失效概率。

3 工程案例

江西省九江市德安縣擁有鉛鋅銻礦礦床一座，為了解決選礦廢渣的排放堆積問題建立尾礦庫。德安尾礦庫采用上游法筑壩，壩頂標高64.30m，初期壩坡比為1：2，有 3 座堆積子壩，1#子壩上、下游坡比為 1：2，2#、3#子壩的上下游坡比為1：4，堆積壩目前標高75.00m。

表1 尾礦壩各土層物理力學參數(shù)

圖1 德安尾礦壩模型

3.1 有限元模型及計算參數(shù)

本文通過 GEOSTUDIO中 Seep/w、Quake/w與Slope/w模塊耦合建立工程地質模型來進行模擬分析。為了能更直觀展現(xiàn)，德安尾礦壩有限元模型一共劃分為2 314個節(jié)點和2 211個網(wǎng)格大小為1.00m的四邊形和三角形混合單元，并設置了A、B、C、D四個地震響應歷程點監(jiān)測地震過程中的位移和加速度，工程地質模型和特征點位置見圖1。邊界調節(jié)為：上游水位（左側）22.50m，下游水位（右側）10.00m，各土層所采用的物理力學參數(shù)見表1。首先在Seep/w模塊建立德安尾礦壩壩體模型，進行飽和滲流分析，隨后將滲流分析結果導入Quake/w模塊，輸入地震波，進行動力響應分析。

3.2 地震作用下尾礦壩的動力響應分析

尾礦壩動力穩(wěn)定計算分析需要在計算程序中輸入地震動時程曲線（地震加速度時程曲線），且不同的地震動時程曲線對壩體動力計算結果影響較大。根據(jù)《中國地震動峰值加速度區(qū)劃圖》（GB 18306-2015），該尾礦所在區(qū)地震設防烈度7度，設計基本地震加速度值為0.10g，需要對加載地震波進行人工合成。根據(jù)德安尾礦壩的場地信息，以場地反應譜為目標譜生成人工加速度時程曲線，卓越周期為0.2s，時長為10s，時間步長為0.02s，如圖2所示。

圖2 地震加速度時程曲線

本次地震歷程響應點選取A、B、C、D四點，分別代表尾礦壩初期壩的壩底、壩頂和堆積壩的壩中、壩頂四個部分。根據(jù)該尾礦壩所在的位置加入地震參數(shù)，進行位移響應分析，位移時間圖見圖3。

從圖中可以得到4個響應點中最先產生X方向位移的是堆積壩壩頂?shù)腄點，然后是堆積壩壩中的C點，最后是初期壩的A點和B點。位于尾礦壩頂點的位移峰值大于其它各點的峰值，四個地震響應點的位移變化趨勢基本相一致?？梢灶A判，在同一地震波的影響下，尾礦壩壩頂受地震的影響最大。隨著地震事件的增加，A、B、C、D四點處的峰值位移都在不斷的增加。其中1.7s和6.7s時分別產生了向左和向右的最大位移，為了更加直觀的看到尾礦壩的情況，分析了這兩個時刻的尾礦壩相對位移圖（見圖4、圖5）。

圖3 X-位移時間圖

分析A、B、C、D四點X方向速度見圖6，可以看出A點的峰值速度出現(xiàn)在1.98s，其值為0.076m/s；B點峰值速度出現(xiàn)在4.50s，其值為0.100m/s；C點峰值速度出現(xiàn)在4.05s，其值為0.067m/s；D點峰值速度出現(xiàn)在4.33s，其值為0.120m/s，四個點的峰值速度有所差別，但是B、C、D點峰值速度出現(xiàn)的時間點大致相同，堆積壩壩頂?shù)腄點在X方向速度波動最為劇烈，各響應點最大位移從壩頂方向到壩底逐漸減小，四點總體的震動趨勢與地震波行較為相似，而尾礦壩不同位置在同一地震波影響下反映出的劇烈程度不一樣。圖7為地震過程中尾礦壩的總應力隨時間變化的分布圖，震后總應力大小基本沒有發(fā)生變化，堆積壩壩頂D的總應力最大，初期壩壩底A的總應力最低。

圖4 1.2s左側最大位移

圖5 8.8s右側最大位移

圖6 地震歷程中速度時間變化圖

圖7 地震歷程中總應力時間變化圖

圖8 最危險滑面平均加速度時程圖

隨后對尾礦壩進行穩(wěn)定性分析，將Seep/w、Quake/w和Slope/w模塊進行耦合，在地震動力分析的基礎上，使用進入進出法分析、計算尾礦壩的安全系數(shù)，不同時刻尾礦壩最危險滑面平均加速度和壩體安全系數(shù)時程曲線如圖8和圖9所示。在初始狀態(tài)下該尾礦壩的安全系數(shù)為1.36，一開始安全系數(shù)趨于穩(wěn)定，之后便有較大的起伏，在4.7s時，安全系數(shù)最低為0.92，根據(jù)首次破壞準則，此時尾礦壩可能發(fā)生潰壩危險。且安全系數(shù)時程圖跟地震波型圖亦十分相似，安全系數(shù)的峰值與地震波的峰值存在一定的滯后性。

3.3 可靠度的計算

對尾礦壩進行可靠度分析，結合德安尾礦壩工程資料，并參考鄭欣[10]、張揚等[11]，確定土體黏聚力及摩擦角的均值及標準差，結果如表2所示。

圖9 最危險滑面安全系數(shù)時程圖

表2 黏聚力及摩擦角參數(shù)統(tǒng)計表

尾礦壩破壞可靠度分析方法使用了MCS法，在Seep/w模塊中以參數(shù)均值建立尾礦壩有限元分析模型、劃分有限元網(wǎng)格、設置邊界條件，在Quake/w模塊中建立動力模型，在Slope/w模塊建立穩(wěn)定性分析模型，在Seep/w、Quake/w和Slope/w模塊耦合計算的基礎上，通過設置滑移面進入進出口的位置來進行邊坡穩(wěn)定分析，求出安全系數(shù)。將尾礦壩穩(wěn)定性有限元分析模型另存為名為“安全系數(shù).xml”的計算源文件。當尾礦壩土體輸入?yún)?shù)表示為標準正態(tài)隨機變量之后，可將輸入?yún)?shù)特性值分別賦給尾礦壩穩(wěn)定性有限元分析模型，生成新的尾礦壩穩(wěn)定性分析“安全系數(shù).xml”計算文件。最后使用MATLAB軟件在DOS環(huán)境下直接調用GEOSTUDIO軟件內部程序對新生成的“安全系數(shù).xml”文件進行一體化邊坡穩(wěn)定性有限元批處理分析，采用摩根斯坦-普萊斯法計算最危險滑動面，無需打開GEOSTUDIO軟件界面，大大地提高了計算效率。最后通過可靠度理論來計算失穩(wěn)風險概率，實驗次數(shù)為10 000次，倘若人為的對有限元軟件進行重復性的操作來計算10 000次的話，可能需要幾個星期，這里使用批處理方法只需普通電腦自動計算幾天，無需人為干涉。經計算，德安尾礦壩在7度地震作用下的失穩(wěn)風險概率為0.37%。

4 結論

（1）本文結合動力有限元模型和結構可靠度理論，分析了德安尾礦壩在7度地震作用下的動力穩(wěn)定性，探究了德安尾礦壩地震過程中的位移、速度、應力、安全系數(shù)等變化情況，比傳統(tǒng)的擬靜力方法更能反映尾礦壩在地震作用下的動力響應。在地震作用中，不同時刻壩體穩(wěn)定性隨地震加速度的變化而變化，隨著臨界滑面平均加速度的增加，尾礦壩安全系數(shù)會逐漸降低。

（2）運用蒙特卡羅隨機抽樣方法，并編寫MATLAB程序與GEOSTUDIO軟件相結合，將商業(yè)有限元軟件GEOSTUDIO當做黑匣子直接調用，進行一體化尾礦壩邊坡穩(wěn)定性有限元批處理計算，可以快速的進行了10 000次尾礦壩的穩(wěn)定性分析，大大提高了尾礦壩的計算效率。

（3）經可靠度分析，德安尾礦壩的失穩(wěn)概率為0.37%。盡管德安尾礦壩在7度地震工況下的失穩(wěn)概率較低，與尾礦壩設計規(guī)范給出的允許安全系數(shù)和破壞概率相比，該尾礦壩發(fā)生的失穩(wěn)破壞可能性較小，但是仍需要設計一些尾礦壩安全防護抗震措施。

（4）本文通過理論分析與數(shù)值模擬相結合的方法對地震作用下尾礦壩的穩(wěn)定性可靠度進行了計算分析，研究結果也具有一定的實際意義，但是分析過程中仍做了一些相應的簡化，因此，作者認為針對該方面的研究還有待進一步深入：1）影響尾礦壩穩(wěn)定性的不確定因素較多，本文僅考慮了尾礦材料參數(shù)粘聚力和內摩擦角的不確定性，未考慮到地震以及降雨作用的不確定性，因此，在今后的研究中可進一步考慮更多參數(shù)的不確定性，以使計算結果更加貼合實際；2）本文采用的GEOSTUDIO軟件為2D軟件，不能反映尾礦壩最真實的受力變形情況。因此，在今后的研究中可借助其他3D有限元軟件進行計算分析，結果將會更加接近實際。