許志發(fā)，龍 翠，唐永駿，高亞偉，趙懷剛，王光進

（1.昆明理工大學 國土資源工程學院，云南 昆明 650093；2.云南鋁業(yè)股份有限公司安全環(huán)保部，云南 昆明 650502；3.昆明理工大學建筑工程學院，云南 昆明 650504）

0 引 言

礦山三大控制性工程之一的尾礦庫是一種特殊的工業(yè)建筑物，是礦山企業(yè)在生產(chǎn)建設過程中的重要設施。中國是鎢資源大國，儲量占世界總量的68%，排名世界第一[1]。但隨著經(jīng)濟的發(fā)展，對鎢的需求量不斷增長，鎢尾礦庫的數(shù)量和等級逐年增加[2]。因此，尾礦庫的運行狀況不僅與礦山能否順利進行生產(chǎn)建設，還關(guān)系到尾礦庫下游的周邊環(huán)境安全和群眾的生命財產(chǎn)安危[3-4]。自1910年以來，已經(jīng)有數(shù)百起尾礦庫事故在世界各地發(fā)生[5]。其中三分之一的尾礦庫屬“病險庫”，因尾礦庫是具有高勢能的人造危險源，一旦失事，易造成非常嚴重的事故。例如，2008年9月8日，位于山西省的新塔礦業(yè)有限公司發(fā)生尾礦庫事故造成281人死亡，直接經(jīng)濟損失9 619.2萬元的特別重大事故；2010年9月21日，位于廣東的信宜紫金礦業(yè)有限公司所屬的銀巖錫礦尾礦庫發(fā)生尾礦庫潰壩事故造成壩體下游4人死亡，直接經(jīng)濟損失約1 900萬元，這些事故產(chǎn)生了極其惡劣的社會影響[6-8]。這些都是人為責任事故，如果能夠規(guī)避尾礦庫建設過程中的違規(guī)行為，并對筑壩過程中的應力變形進行實時有效的監(jiān)控，異常變形和壩體開裂等異?，F(xiàn)象能夠被及早發(fā)現(xiàn)，筑壩作業(yè)中的不規(guī)范得到調(diào)整，礦山的生產(chǎn)建設和下游群眾的生命財產(chǎn)安全將得到有效的保障。

尾礦庫是由各種不同粒徑尾礦顆粒的堆砌形成的人造土體，近年來國內(nèi)外學者對其應力變形特性進行了大量的試驗和探討。Floresc[9]分別對墨西哥兩種尾粉砂進行共振柱試驗，應變范圍為 10－5～10－2，用Davidenkov模型對試驗結(jié)果進行擬合，得出了兩個參數(shù)的試驗數(shù)據(jù)；袁澤川[10]運用剪應變歸一化和等效線性模型的方法研究圍壓對無黏性土、黏性土和一般性土的最大剪切模量的影響，得到尾礦在不同圍壓下的阻尼比和動模量衰減曲線；曹進海等[11]針對余震對尾礦壩穩(wěn)定性的影響問題，以某上游式尾礦壩為例，基于完全非線性有限差分法，從不同余震烈度、網(wǎng)格變形、位移變化及超孔壓比等方面對尾礦壩進行數(shù)值模擬分析；柳厚祥等[12]研究了尾礦顆粒在靜力、動力和共振情況下的物理力學性質(zhì)，模擬了在不排水有效應力法和排水有效應力法下，對尾礦庫在地震作用下的應力變形進行分析；周健[13]將尾礦庫看作三維兩相體，研究了振動孔隙水壓力增長規(guī)律，進而運用三維兩相有效應力動力分析方法對尾礦庫的應力形變進行分析。

研究選取四川省某鎢尾礦庫，利用Geo-Studio對尾礦庫的應力變形運用有限元法進行模擬，分析尾礦庫在原設計標高和在二期設計標高堆筑下尾礦庫整體的應力應變等狀態(tài)，為尾礦庫的設計和施工提供指導性依據(jù)。

1 應力變形數(shù)學模型

有限元法、離散元法、拉格朗日元法及邊界元法等在尾礦庫、邊坡等巖土工程都有較為廣泛的應用。目前，尾礦庫與邊坡已廣泛應用基于有限元法的數(shù)值模擬軟件對應力應變和穩(wěn)定性進行分析計算。有限元法是將復雜的物體或者系統(tǒng)劃分為有限不重疊的單元，這些單元相互間獨立又互相聯(lián)結(jié)。所有作用的力在每個單元上都等效到節(jié)點上，最后將各變量或其導數(shù)的節(jié)點值和插值函數(shù)組成的線性表達式代替微分方程中的變量。有限元分析法可以用數(shù)值模擬軟件來分析難以人工計算的巨量數(shù)和方程組。

1.1 有限元法基本原理

在進行分析之前需要對物體或系統(tǒng)進行網(wǎng)格劃分，在二維有限元分析中可以選用三角形、矩形、四邊形等多邊形進行系統(tǒng)分析。在這幾種網(wǎng)格劃分中運用三角形進行網(wǎng)格劃分，具有網(wǎng)格密度大、求解精度高等優(yōu)點。常應變狀和剛體位移態(tài)是單元的位移函數(shù)必須滿足的條件，位移必須連續(xù)[14]。因此，函數(shù)u（x，y），v（x，y）可寫為式（1）所示。

將位移值和坐標值代入式（1），得到6個求解方程，將其聯(lián)立獲得各個系數(shù)的求解值。將其再代入式（1）得到式（2）。

式中：｛｝δ =[u1，v1，u2，v2，u3，v3]T為 節(jié) 點 位 移 矩 陣 ；稱為形狀函數(shù)矩陣，其，i＝1，2，3稱為形狀（插值）函數(shù)。

確定位移函數(shù)之后，求得單元的應變和應力。單元應變?nèi)缡剑?）所示。中

式中：[B]為單元應變矩陣。

將式（3）代入，求得單元應力如式（4）所示。

式中：[D]為彈性矩陣。

設壩體發(fā)生虛位移，網(wǎng)格單元虛位移為{δ*}，與之的網(wǎng)格單元虛應變?yōu)閧ε*}，則有：

式中：A 為單元面積，m2；t為單元厚度，m；S 為單位虛位移，m；T為轉(zhuǎn)置矩陣。

經(jīng)變換得：

式中：{U}為總體位移列陣；[K]為總體剛度矩陣；{P}為總體荷載列陣。

總體剛度方程通過邊界約束條件進行修正后，總體位移列陣{U}可以求解得到，然后計算各單元的應力和應變分量，進而對物體或系統(tǒng)的應力和應變進行模擬分析。

1.2 初始模量

不同初始應力下，土的初始模量Ei是變化的。根據(jù)三軸試驗，若把Ei和相應的σ3用大氣壓力Pa歸一成無量綱數(shù)，可得到關(guān)系式（7）：

式中：k、n是取決于土的性質(zhì)參數(shù)，分別稱為模量系數(shù)和模量指數(shù)，可根據(jù)三軸試驗測出。

1.3 體積模量和泊松比

材料變形中，反應體積變化的參數(shù)一般為體積模量B或泊松比μ，兩者之間的關(guān)系見式（8）所示。

體積模量和泊松比也可以表示為與應力狀態(tài)有關(guān)，即：

式中：kb、m、G、F、D 皆為試驗常數(shù)；c為內(nèi)聚力，kPa；φ 為內(nèi)摩擦角，（°）。

2 尾礦庫的應力變形模擬分析

2.1 Geo-Studio概況

加拿大GEO-SLOPE公司在20世紀70年代開發(fā)了Geo-Studio，該軟件是可以應用于巖土、地質(zhì)、采礦等各行各業(yè)數(shù)值分析的軟件之一。它具備操作方便、模塊化計算、建模簡便、計算分析快等優(yōu)點，礦業(yè)、巖體工程等許多領(lǐng)域已廣泛應用。SIGMA/W是Geo-Studio系列軟件的重要組成部分，是應力變形分析中專業(yè)的有限元分析軟件，即可實現(xiàn)線彈性分析，也可以完成非線性彈塑性分析等高度復雜的各種分析。Geo-Studio軟件采用增量形式進行求解，可以選擇線彈性、各向異性彈性、雙曲線彈性、彈塑性等多種材料模型，能夠很好地進行邊坡、壩體等問題的模擬計算分析。

2.2 應力變形模擬

2.2.1 尾礦庫概況

尾礦庫在2012年完成設計和施工，由于原設計不能滿足礦山擴大生產(chǎn)的要求，2017年進行尾礦庫增高擴容工程的設計。地貌單元為山地中山，溝谷相間地形，兩側(cè)山坡地勢較緩，隨地形變化25°～35°，表面大部為第四系殘坡積層覆蓋，局部可見基巖出露，谷底為第四系沖洪積層覆蓋。尾礦庫堆積壩一期設計最終標高2 035 m（總壩高92 m），使用段溝底縱坡8.7%。因一期設計不能滿足未來尾礦的排放需求，因此對尾礦庫進行增高擴容，增高擴容后的設計壩高為147 m（即加高55 m），加高后的壩頂平臺標高為2 090 m。根據(jù)《尾礦設施設計規(guī)范》（GB50863—2013）[15]對尾礦庫的分類，原尾礦庫為三等庫，加高壩高后，按壩高考慮尾礦庫為二等，故尾礦壩等尾礦庫主要構(gòu)筑物的級別為二級。

2.2.2 尾礦庫穩(wěn)定性計算

庫容增加后尾礦庫的穩(wěn)定性應該滿足規(guī)范要求，因此針對庫容增高后的穩(wěn)定性進行研究。由于尾礦庫所在區(qū)域的地震基本烈度為Ⅷ度，所以需考慮地震作用下壩體的穩(wěn)定性。本文選擇三種運行工況開展尾礦庫穩(wěn)定性分析，一是尾礦庫處于正常運行工況；二是在洪水運行工況，其洪水位按二等尾礦壩的最小安全超高（1.0 m）考慮；三是在特殊運行工況，即考慮地震基本烈度Ⅷ度時的穩(wěn)定性。前兩種工況運用Slide軟件進行尾礦庫穩(wěn)定性分析。采用畢肖普法對尾礦庫進行穩(wěn)定性分析，計算結(jié)果如圖1、表1所示。結(jié)果顯示在壩高抬高至2 090 m時，尾礦庫穩(wěn)定性仍滿足規(guī)范要求。

圖1 2 090 m時兩種工況畢肖普法計算結(jié)果Fig.1 Bishop method results of two conditions at 2 090 m

表1 高程2 090 m時的抗滑穩(wěn)定性計算結(jié)果表Tab.1 Results of anti-sliding stability at an elevation of 2 090 m

地震工況下的尾礦庫動力穩(wěn)定性分析采用Geo-Studio系列軟件的QUAKE/W模塊進行模擬。其可以用來解決因地震而產(chǎn)生的動態(tài)載荷以及碰撞作用下的沖擊載荷等作用下的動力穩(wěn)定性問題。應用該模塊，可以同時考慮因地震產(chǎn)生而作用在土體中的孔隙水壓力增長和動應力，能夠研究地震對土工建筑物強度穩(wěn)定性的影響。在動力穩(wěn)定性分析過程中，選用汶川地震波和人工合成地震波兩組典型地震波對尾礦庫開展地震工況的邊坡穩(wěn)定性分析。

采用有限元極限平衡法對地震作用各時段進行穩(wěn)定性計算，得到汶川地震波和人工合成地震波兩組典型地震波對壩體作用的最小安全系數(shù)。由圖2可知，兩種地震波作用過程中，壩體安全系數(shù)均大于1，汶川地震波作用過程中最小值為1.320，人工合成地震波作用過程中最小值為1.212，均符合規(guī)范要求。

2.2.3 幾何模型與計算參數(shù)

壩體應力變形有限元分析計算剖面選擇最不安全剖面（即沿谷底計算剖面）在堆高至2 043 m和2 090 m時開展其應力變形分析，而材料分區(qū)仍然采用二維穩(wěn)定性的概化分區(qū)（即距子壩30 m的位置作為尾細砂與尾粉砂的分界線；距子壩180 m的位置作為尾粉砂與尾粉土的分界線；距子壩350 m的位置作為尾粉土與尾粉質(zhì)黏土的分界線）。同時，計算的2 043 m和2 090 m高程時的浸潤線埋深采用的是剛好滿足規(guī)范規(guī)定最小安全系數(shù)時的臨界浸潤線埋深，即尾礦庫在堆高至2 043m時其壩體的臨界浸潤線埋深值為14.0 m，而堆高至2 090 m時壩體的臨界浸潤線埋深值為25.0 m。依據(jù)室內(nèi)土工試驗數(shù)據(jù)分析，同時與礦山提供的勘察報告進行參照和國內(nèi)外尾礦庫設計中相關(guān)參數(shù)，進行全面的對比分析，模擬分析中所需的各材料的參數(shù)如表2所示。

2.2.4 邊界條件與幾何劃分

考慮到尾礦庫的實際情況，在分析中邊界條件的確定如圖3所示（一期：2 043 m，二期：2 090 m），由于在實際中庫尾一般為山體、硬巖等自然地形，因此將左側(cè)邊界確定為水平方向的約束，即U=0，V≠0；尾礦庫底部在堆壩前都需要提前將淤泥、腐殖土、雜填土等清除干凈，使底部形成一個相對穩(wěn)定的地基，因此將底部邊界的水平和垂直約束都設定為約束，即U=0，V=0。尾礦庫原堆積壩高2 043 m和最終堆積壩高2 090 m時應力變形計算邊界條件和網(wǎng)格劃分結(jié)果如圖3所示。

表2 尾礦壩應力變形計算參數(shù)Tab.2 Calculation parameters for tailings dam stress deformation

圖3 沿谷底計算剖面不同堆積壩高應力變形計算圖Fig.3 Computational map of high stress deformation of different accumulation dams calculated along the valley floor

2.2.5 尾礦庫應力變形數(shù)值計算結(jié)果

在計算中考慮到地下水對壩體應力變形的作用，分別對正常工況下壩高堆至2043m（壩高100 m）和最終壩高2 090 m時（壩高147 m）的水平、豎直方向位移，水平、豎直方向總應力和剪應力總應力進行應力和變形數(shù)值模擬，計算結(jié)果見圖4、圖5。

圖4 壩高堆至2 043 m時正常工況下位移等值云圖、總應力云圖及剪應力總應力云圖Fig.4 Displacement equivalent cloud image,total capacity cloud image and shear stress total stress cloud image under normal conditions of up to 2 043 m dams

3 結(jié) 論

本文采用數(shù)值模擬軟件對某鎢尾礦庫不同堆高下的穩(wěn)定性和應力變形進行模擬，分析其穩(wěn)定性、沉降和應力應變，研究了其變化規(guī)律。

（1）分析庫容增加下正常運行工況、洪水運行工況和地震運行工況下三種不同工況尾礦庫的穩(wěn)定性，分析結(jié)果表明庫容和壩高增加的情況下尾礦庫的穩(wěn)定性符合設計規(guī)范要求。

（2）壩高堆至2 043 m（總壩高100 m）和最終壩高2 090 m（總壩高147 m）時，壩體變形規(guī)律基本一致。壩體水平向變形主要向著壩體下游，最大水平向位移出現(xiàn)在壩坡中上部，豎向變形自上往下逐漸降低，最大豎向變形出現(xiàn)在壩坡頂部。

（3）不同堆高情況下在初期壩內(nèi)坡腳都出現(xiàn)一定的應力集中現(xiàn)象，但壩體內(nèi)大部分區(qū)域剪應力不大。

圖5 最終壩高2 090 m時正常工況下位移等值云圖、總應力云圖、剪應力總應力云圖Fig.5 Displacement equivalent cloud image,total capacity cloud image and shear stress total stress cloud image under normal conditions of up to 2 090 m dams

（4）堆積時向下固結(jié)沉降和尾礦沉降壓縮導致側(cè)向膨脹，引起水平方向位移。隨壩體增高，壩體變形逐漸增大，且下層堆積壩會在上層堆積后出現(xiàn)局部“隆起”。同時，干灘面長度減小也會使壩體水平位移和垂直位移有所增大。

[1] 鄔海濱，李繼福，袁勤智，等.回收某硫化礦尾礦中白鎢的選礦試驗研究[J].中國鎢業(yè)，2017，32（3）：36-41.WU Haibin，LI Jifu，YUAN Qinzhi，et al.Beneficiation of scheelite from a tailings of a sulphide ore[J].China Tungsten Industry，2017，32（3）：36-41.

[2] 李茂林，王 旭.甘肅某白鎢浮選尾礦再回收白鎢試驗研究[J].中國鎢業(yè)，2015，30（2）：17-20.LI Maolin，WANG Xu.Study on the recovery of scheelite from a scheelite flotation tailing in Gansu province[J].China Tungsten Industry，2015，30（2）：17-20.

[3] 魏作安，尹光志，沈樓燕，等.探討尾礦庫設計領(lǐng)域中存在的問題[J].有色金屬（礦山部分），2002，54（4）：44-45，48.WEIZuoan，YINGuangzhi，SHENLouyan，etal.Discusstheproblems existing in the design of tailings ponds [J].Nonferrous Metals（Mining Section），2002，54（4）：44-45，48.

[4] MILTON A，COOKE J A，JOHNSON M S.A comparison of cadmium in ecosystems on metalliferous mine tailings in Wales and Ireland[J].Water，Air and Soil Pollution，2004，153（1）：157-172.

[5] OZCANNT，ULUSAY RESAT，ISIKNS.A study on geotechnical characterization and stability of downstream slope of a tailings dam to improve its storage capacity （Turkey）[J].Environmental Earth Sciences，2013，69（6）：1871-1890.

[6] 謝旭陽，田文旗，王云海，等.我國尾礦庫安全現(xiàn)狀分析及管理對策研究[J].中國安全生產(chǎn)科學技術(shù)，2009，5（2）：5-9.XIE Xuyang，TIAN Wenqi，WANG Yunhai，et al.Analysis of safety status and management countermeasures of tailings in China[J].

[7] 張興凱，王啟明，相桂生.金屬非金屬尾礦庫安全現(xiàn)狀及分析[J].中國安全生產(chǎn)科學技術(shù)，2006（2）：60-62.ZHANG Xingkai，WANG Qiming，XIANG Guisheng.Safety status and analysis of metal nonmetal tailings storage[J].Journal of Safety Science and Technology，2006（2）：60-62.

[8] 李青石，李庶林，陳際經(jīng).試論尾礦庫安全監(jiān)測的現(xiàn)狀及前景[J].中國地質(zhì)災害與防治學報，2011，22（1）：99-106.LI Qingshi，LI Shulin，CHEN Jijing.Current status and prospect of tailings dump safety monitoring[J].The Chinese Journal of Geological Hazard and Control，2011，22（1）：99-106.

[9] FLORES O C，ROMO M P O.Dynamic behavior of tailings[C]//Shamsher Prakash.Geotechnical Engineering Commons.California：University of Missouri—Rolla，2001：1-6.

[10] 袁澤川.尾礦的動力應力-應變特性 [J].工業(yè)建筑，2017，47（10）：139-145.YUAN Zechuan.Dynamic stress-strain characteristics of tailings[J].Industrial Construction，2017，47（10）：139-145.

[11] 曹進海，胡 軍.基于FLAC3D余震對尾礦壩穩(wěn)定性影響的研究[J].中國鎢業(yè)，2017，32（4）：68-73.CAO Jinhai，HU Jun.Study on the influence of FLAC3Daftershocks on the stability of tailings dam [J].China Tungsten Industry，2017，32（4）：68-73.

[12] 柳厚祥，廖 雪，李 寧，等.高尾礦壩的有效應力地震反應分析[J].振動與沖擊，2008，27（1）：65-70，92.LIU Houxiang，LIAO Xue，LI Ning，et al.Effective stress seismic response analysis of high tailings dam[J].Journal of Vibration and Shock，2008，27（1）：65-70，92.

[13] 周 健.尾礦壩在地震作用下的三維兩相有效應力動力分析[J].工程抗震，1995（1）：39-43.ZHOU Jian.Three-dimensional two-phase effective stress dynamic analysisoftailingsdamunderseismicaction[J].EarthquakeResistant Engineering，1995（1）：39-43.

[14] 李慧謙.尾礦庫堆壩試驗研究及穩(wěn)定性分析[D].重慶：重慶大學，2008.LI Huiqian.Experimental study and stability analysis of tailings dam[D].Chongqing：Chongqing University，2008.

[15]中華人民共和國國家標準編寫組.GB50863—2013尾礦設施設計規(guī)范[S].北京：中國計劃出版社，2013：1-50.