蘇振寧，邵龍?zhí)?/p>

(大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室， 遼寧 大連 116024)

尾礦庫具有體量大，危險性高，失穩(wěn)破壞性強等特點，是重大危險源之一。合理保障尾礦壩的穩(wěn)定性和耐久性，對于保障礦業(yè)可持續(xù)生產(chǎn)，保證庫區(qū)下游的生產(chǎn)生活具有重要的工程價值。對尾礦壩邊坡一般選取典型斷面以平面應(yīng)變問題形式利用極限平衡法或強度折減法進(jìn)行穩(wěn)定分析[1-5]。但尾礦庫一般為建于山谷中的典型三維結(jié)構(gòu)體，二維穩(wěn)定分析不能真實反映尾礦壩的穩(wěn)定性。

陳順良等[6]采用三維有限元分析了尾砂開采對壩體位移和應(yīng)力的影響。鄧紅衛(wèi)等[7]建立三維數(shù)值模型，采用滲流-應(yīng)力耦合機(jī)理，研究了干灘長度對尾礦庫壩體穩(wěn)定性的影響。李洪梁等[8]利用無人機(jī)測量采集數(shù)據(jù)，建立三維模型計算浸潤線分布，采用極限平衡法分析壩體穩(wěn)定性。以上研究均對尾礦壩建立了三維模型，但穩(wěn)定性僅采用定性的應(yīng)力分析或二維安全系數(shù)結(jié)果。采用三維流固耦合模型有限元強度折減法能計算尾礦壩三維穩(wěn)定安全系數(shù)從而對穩(wěn)定性進(jìn)行分析[9-10]，還能對尾礦壩三維失穩(wěn)形態(tài)進(jìn)行研究[11]。但有限元強度折減法存在邊坡穩(wěn)定判據(jù)影響安全系數(shù)的問題[12]，并且尾礦壩為非均質(zhì)復(fù)雜邊坡，能否對所有土層采用同一折減系數(shù)仍有待研究[2]。

有限元極限平衡法[13]根據(jù)“真實”應(yīng)力場結(jié)果和基于極限平衡狀態(tài)的安全系數(shù)定義，可以采用優(yōu)化算法搜索區(qū)域內(nèi)的臨界滑動面。于斯瀅等[2]采用二維有限元極限平衡法分析了尾礦壩穩(wěn)定性，取得了很好的效果。三維有限元極限平衡法[14]同樣具備計算效率高，無需折減迭代計算，滑動面和安全系數(shù)明確的特點，能給出局部安全系數(shù)分布，能為尾礦壩安全設(shè)計提供可靠依據(jù)和技術(shù)支持。本文即采用三維有限元極限平衡法對所研究的尾礦壩邊坡開展三維穩(wěn)定分析，并與有限元強度折減法的分析結(jié)果相比較，以期更好指導(dǎo)工程實踐。

1 邊坡穩(wěn)定分析方法

1.1 二維極限平衡法

根據(jù)《尾礦庫安全技術(shù)規(guī)程》[15](AQ 2006—2005)要求對典型斷面使用瑞典圓弧法和Bishop法判斷壩體穩(wěn)定性。本文選取尾礦壩最大斷面分別采用瑞典圓弧法、Bishop法和Morgenstern-Price(M-P)法計算安全系數(shù)。瑞典圓弧法和Bishop法要求滑動面為圓形；M-P法滿足整體和每個土條的力與力矩平衡可適用于任意形狀的滑動面。極限平衡法計算軟件采用GeoStudio。

1.2 三維有限元強度折減法

有限元強度折減法在有限元計算過程中不斷將土體的強度參數(shù)進(jìn)行折減，直至邊坡到達(dá)臨界狀態(tài)。當(dāng)邊坡達(dá)到臨界狀態(tài)時，強度參數(shù)的折減系數(shù)也就是安全系數(shù)。該安全系數(shù)表征著材料的強度儲備，基于摩爾-庫侖強度準(zhǔn)則的強度折減法表達(dá)式為：

(1)

式中：c為土體黏聚力；φ為土體摩擦角；K為強度參數(shù)折減系數(shù)，當(dāng)邊坡達(dá)到臨界狀態(tài)，K即為安全系數(shù)。

在存在滲流的尾礦壩穩(wěn)定分析中，有限元強度折減可采用流固耦合模型，進(jìn)行滲流計算。強度折減法不需要提前設(shè)定滑動面，可通過臨界狀態(tài)時位移增量云圖直接判定滑動面。

1.3 三維有限元極限平衡法

三維有限元極限平衡法基于“真實”的有限元應(yīng)力場，有限元計算采用彈塑性-流固耦合模型，在計算中無需折減土體強度參數(shù)。該方法通過構(gòu)造三角形網(wǎng)格滑動面，在有限元應(yīng)力場的基礎(chǔ)上根據(jù)整體安全系數(shù)定義計算滑動面安全系數(shù)，并通過優(yōu)化方法尋找區(qū)域內(nèi)最小安全系數(shù)以及其所對應(yīng)的滑動面。該方法優(yōu)點在于采用真實應(yīng)力場計算，無需迭代計算，沒有收斂性問題，物理意義明確，能給出準(zhǔn)確的最危險滑動面和對應(yīng)的安全系數(shù)以及滑動面局部安全系數(shù)分布，計算簡單易于推廣和工程應(yīng)用。本文采用MATLAB軟件編寫了三維有限元極限平衡法程序，用以計算尾礦壩邊坡安全系數(shù)。有限元極限平衡法由以下三個部分組成：

1.3.1 安全系數(shù)定義

三維有限元極限平衡法局部安全系數(shù)和整體安全系數(shù)分別被定義為：

(2)

(3)

式中：d為滑動面上滑動方向；τf為抗滑剪應(yīng)力矢量，方向與滑動方向d相反，其值根據(jù)摩爾庫侖強度準(zhǔn)則計算；τ為滑動面上剪應(yīng)力矢量；s為滑動面面積。局部安全系數(shù)是抗剪應(yīng)力與剪應(yīng)力在滑動方向上的投影的比值，其物理意義是一點在滑動方向上距離極限平衡狀態(tài)的衡量。通過邊坡整體極限平衡條件和積分中值定理[13]，可得整體安全系數(shù)，其表征滑動面上各點局部安全系數(shù)的中值。

對于存在孔隙水壓力的邊坡的穩(wěn)定分析，強度計算應(yīng)采用基于有效應(yīng)力的方法，摩爾-庫侖強度準(zhǔn)則為：

τf=σ′tanφ′+c′

(4)

式中：σ′為滑動面上的有效正應(yīng)力；φ′為有效摩擦角；c′為有效黏聚力。

有限元在流固耦合計算中應(yīng)力采用的是有效應(yīng)力與孔隙水壓力相結(jié)合的形式。所以從有限元中導(dǎo)出的應(yīng)力就是有效應(yīng)力，可以使用式(3)直接計算。只需保證強度參數(shù)為有效應(yīng)力強度參數(shù)。

1.3.2 滑動面構(gòu)造

根據(jù)定義式(3)，有限元極限平衡法需要對滑動面上的變量進(jìn)行積分，數(shù)值積分需要對積分域進(jìn)行離散，采用三角形網(wǎng)格形式對滑動面進(jìn)行劃分。

采用改進(jìn)橢球形表征滑動面。相比于使用圓心坐標(biāo)(x，y，z)和三軸半徑(a，b，c)作為參數(shù)的傳統(tǒng)橢球滑動面，改進(jìn)橢球形滑動面參數(shù)增加了橢球面繞z軸的旋轉(zhuǎn)角度α和中間柱體寬度l，如圖1所示。考慮到實際工程中受地形以及建模方向影響，橢球形滑動面的水平主軸方向未必是x和y坐標(biāo)軸方向，所以增加了旋轉(zhuǎn)角度α。在三維問題中，當(dāng)邊坡沿垂直于坡體斷面方向延伸，邊坡的安全系數(shù)越趨近于二維平面應(yīng)變結(jié)果，其滑動面也更趨近于管狀。所以增加了中間柱體寬度l。

圖1 改進(jìn)橢球形滑動面

1.3.3 最危險滑動面搜索

三維有限元極限平衡法的滑動面搜索方法選擇粒子群算法作為優(yōu)化方法。粒子群算法是一種智能優(yōu)化算法，具有良好的計算性能和出色的應(yīng)用效果，在邊坡最危險滑動面的搜索中有良好的效果[16]。

本文使用MATLAB的全局優(yōu)化工具箱進(jìn)行粒子群算法的最危險滑動面搜索，改進(jìn)橢球形滑動面參數(shù)作為搜索參數(shù)，整體安全系數(shù)作為目標(biāo)函數(shù)，具體流程見圖2。

圖2 粒子群算法流程圖

2 尾礦庫三維計算模型

2.1 工程概況

某待建尾礦庫采用上游筑壩法建設(shè)工藝，庫區(qū)內(nèi)匯水面積8.9 km2，平均坡度達(dá)到8.92%，溝長約3.68 km。初期壩壩址選擇在溝口處，初期壩采用庫內(nèi)開采的風(fēng)化巖石筑透水壩，頂寬6 m，內(nèi)外坡均為1∶2，壩頂標(biāo)高1 665 m，壩底標(biāo)高1 615 m，最大壩高50.0 m。后期子壩壩采用尾礦加高，坡比1∶5.0，終期標(biāo)高為1 813 m，總壩高198 m。當(dāng)壩高達(dá)到最終堆積標(biāo)高1 813 m標(biāo)高時，為二等尾礦庫。

為確保壩體的滲透穩(wěn)定，降低浸潤線，初期壩頂以上每升高15 m，增設(shè)水平排滲盲溝，尾礦堆積壩上共需設(shè)9條水平排滲盲溝。排滲盲溝與堆積壩壩面水平距離為175 m。排滲盲溝由水平排滲盲溝和垂直盲溝的導(dǎo)水管組成，盲溝在尾礦沉積灘上開挖，平行壩軸線，每50 m盲溝設(shè)一條導(dǎo)水管，盲溝由反濾料外包土工布組成。

2.2 尾礦庫有限元模型

尾礦料的物理力學(xué)性質(zhì)和尾礦庫堆積材料的概化分層根據(jù)相同礦區(qū)另一尾礦庫勘察資料和工程經(jīng)驗綜合確定。因為礦物開采自統(tǒng)一礦區(qū)且采用相同的選礦工藝，所以此勘察可以作為新建尾礦庫中尾礦料參數(shù)的取值依據(jù)。尾礦材料分為尾粉砂，尾粉土和尾粉質(zhì)黏土(見表1)。尾礦庫區(qū)地形特征和尾礦材料概化分層，以及排滲盲溝設(shè)計建立尾礦庫幾何模型如圖3所示，其中A-A是所考察的典型斷面位置。

表1 尾礦材料計算參數(shù)表

圖3 尾礦庫幾何模型

二維模型采用GeoStudio軟件分析，穩(wěn)態(tài)滲流計算采用SEEP/W模塊，極限平衡法安全系數(shù)計算采用SLOPE/W模塊。三維模型采用ABAQUS軟件進(jìn)行流固耦合計算。尾礦庫內(nèi)上游根據(jù)不同干灘長度設(shè)置水壓邊界條件，水壓力設(shè)為0 kPa，排滲盲溝和初期壩設(shè)置自由滲出邊界。初期壩和尾粉砂設(shè)置為非飽和狀態(tài)下的參數(shù)，尾礦庫底部為約束位移邊界條件，采用摩爾庫侖強度準(zhǔn)則。計算分別采用三種工況：(1) 不考慮滲流的工況；(2) 洪水工況設(shè)計干灘長度是100 m；(3) 正常生產(chǎn)工況設(shè)計干灘長度是200 m。

3 計算結(jié)果

3.1 滲流計算對比

將三維與二維滲流計算結(jié)果在斷面處對比，如圖4所示。根據(jù)計算，在洪水工況下，三維和二維計算得到的浸潤線未觸及靠近坡頂?shù)膬蓷l排滲盲溝。在正常工況下，三維浸潤線觸及靠近初期壩的四條排滲盲溝，二維浸潤線均未觸及排滲盲溝。同時，在兩種情況下，靠近初期壩位置和初期壩內(nèi)，二維計算浸潤線埋深大于三維結(jié)果。造成以上情況的原因是隨著靠近下游方向，谷口過水面積減小，二維分析無法反映這一特性，造成浸潤線埋深較深。滲流計算結(jié)果顯示，排滲盲溝是確保浸潤線埋深的關(guān)鍵構(gòu)造措施，此尾礦庫管理應(yīng)加強排滲盲溝的監(jiān)控，確保排滲通暢。

圖4 二維與三維尾礦庫滲流場計算結(jié)果

3.2 穩(wěn)定性計算結(jié)果對比

3.2.1 二維極限平衡法

二維極限平衡法計算結(jié)果列于表2中，滑動面如圖5所示。三種計算方法中，M-P法的安全系數(shù)最小，這是因為M-P法不要求滑動面為球形滑動面。與M-P法相比Bishop法誤差在8.3%～11.8%之間，瑞典圓弧法誤差在0.9%～5.1%之間。二維計算結(jié)果均滿足規(guī)范要求。從滑動面可見，M-P法計算得到的滑動面位于下層尾粉質(zhì)黏土區(qū)域面積大，因為尾粉質(zhì)黏土與其他尾礦料相比強度參數(shù)更低。

圖5 正常工況尾礦庫斷面位置的滑動面比較

表2 不同方法安全系數(shù)比較

3.2.2 三維有限元強度折減法

不考慮滲流以及洪水工況和正常工況下，邊坡的破壞形式與滑動區(qū)域相似，所以取正常工況為例進(jìn)行分析。圖6(a)為計算斷面處以計算不收斂的前一個增量步的位移增量云圖表示的滑動面俯視圖。

圖6 正常工況滑動區(qū)域俯視圖

根據(jù)以往的研究[12]，采用特征點位移突變時的強度折減系數(shù)作為安全系數(shù)，是一種合理的有限元強度折減法失穩(wěn)判據(jù)。因為尾礦壩破壞模式導(dǎo)致特征點位移曲線相對平滑，突變點不明顯。將計算不收斂時x方向位移最大的點作為判斷位移突變的特征點，正常工況下特征點的水平位移與折減系數(shù)關(guān)系如圖7所示，水平位移與折減系數(shù)具有不同單位，不能直接判斷突變點，將兩者關(guān)系繪制于方形區(qū)域，并計算曲線段中的斜率變化角，將斜率變化角的極小值點作為水平位移的拐點。這種方法雖然獲得的拐點位置，但仍添加了人為因素(選擇方形區(qū)域繪制曲線)。根據(jù)這種方法計算得到的安全系數(shù)列于表2中，表2中同樣列出了以計算不收斂作為失穩(wěn)判據(jù)的安全系數(shù)，后者相對于前者的誤差在11.2%～11.5%。

圖7 根據(jù)特征點位移突變判據(jù)計算安全系數(shù)

3.2.3 三維有限元極限平衡法

三維有限元極限平衡法搜索得到的滑動區(qū)域如圖6(b)所示，安全系數(shù)列于表2中。有限元極限平衡法滑動面是改進(jìn)橢球形，當(dāng)橢球形底部超過模型區(qū)域，滑動面會自動調(diào)整回模型底部，所以在圖4中滑動面斷面并非橢圓，而是與二維M-P法優(yōu)化得到的滑動面形狀相似。三維有限元極限平衡法能計算局部安全系數(shù)，圖8為正常工況下的局部安全系數(shù)分布云圖，圖中可見不同材料強度參數(shù)對局部安全系數(shù)的影響，尾礦庫區(qū)模型底部的幾何也會對局部安全數(shù)產(chǎn)生影響?；瑒用娴撞康木植堪踩禂?shù)最小，靠近初期壩位置有最大的局部安全系數(shù)，說明此處具有較大的安全余量。

圖8 正常工況滑動面局部安全系數(shù)

3.2.4 對比分析

滲流會導(dǎo)致安全系數(shù)減小，二維極限平衡法中考慮滲流的影響要大于三維有限元方法。

三維有限元方法的安全系數(shù)要大于二維極限平衡法安全系數(shù)，有限元強度折減法安全系數(shù)要大于有限元極限平衡法安全系數(shù)。有限元強度折減法的安全系數(shù)相對偏于危險，尤其是采用不收斂作為判據(jù)時得到的安全系數(shù)。

在滑動面方面，圖4和圖6顯示了三維有限元強度折減法和三維有限元極限平衡法滑動區(qū)域以及滑動深度相似，在考察斷面，三維方法與M-P法也基本一致，僅在滑動面的滑入位置存在略微差異。

有限元強度折減法在確定精確的安全系數(shù)和滑動面都存在問題，并且尾礦庫內(nèi)材料復(fù)雜非均質(zhì)，對不同土層采用同樣的折減系數(shù)存在爭議。有限元極限平衡法能準(zhǔn)確的給出滑動面和對應(yīng)的安全系數(shù)，并且還能對局部安全系數(shù)進(jìn)行考察。兩種方法能相互印證，互為補充，可以更好的分析尾礦壩的三維穩(wěn)定性問題。

4 結(jié) 論

對某尾礦庫進(jìn)行三維數(shù)值建模，并在此基礎(chǔ)上采用有限元極限平衡法和有限元強度折減對尾礦壩進(jìn)行三維穩(wěn)定性分析，并與二維極限平衡法結(jié)果進(jìn)行了對比，具體結(jié)論如下：

(1) 三維穩(wěn)定分析方法相比于二維方法，能反映三維地形和三維滲流場對穩(wěn)定性的影響，避免了典型斷面選取對計算的影響。三維有限元極限平衡法和有限元強度折減法計算結(jié)果顯示尾礦壩潰壩呈現(xiàn)深層滑動整體破壞形式，得到尾礦壩破壞范圍近似橢球形。

(2) 三維有限元極限平衡法根據(jù)其安全系數(shù)定義式能給出明確的安全系數(shù)和對應(yīng)的臨界滑動面，解決了強度折減法中存在的對尾礦壩邊坡不同土層采用同一折減系數(shù)是否合理以及安全系數(shù)難以準(zhǔn)確判定等問題。三維有限元極限平衡法還能顯示出局部安全系數(shù)，更有利于指導(dǎo)工程實踐。

在實際工程中建議使用三維有限元極限平衡和有限元強度折減法相結(jié)合，綜合判斷尾礦壩的穩(wěn)定性。