基于強(qiáng)度折減法的尾礦庫(kù)滲流穩(wěn)定性分析

張 宇，張樹(shù)茂，崔 旋，周漢民

(北京礦冶科技集團(tuán)有限公司，北京100160)

尾礦庫(kù)作為重大危險(xiǎn)源之一，在國(guó)際災(zāi)害事故排名中位列第18位[1-2]，地庫(kù)水位超高引起的失穩(wěn)破壞是引起尾礦庫(kù)事故的主要原因之一[3-4]。上游法筑壩方法往往浸潤(rùn)線較高，大部分壩體均處于飽和狀態(tài)，非常不利于庫(kù)區(qū)的穩(wěn)定性和安全性[5-6]。采取合理的方法計(jì)算上游法尾礦庫(kù)浸潤(rùn)線位置并評(píng)估庫(kù)區(qū)穩(wěn)定性，對(duì)保證庫(kù)區(qū)安全具有重要的理論意義和實(shí)際工程價(jià)值。

極限平衡法和強(qiáng)度折減法是目前使用最為廣泛的兩種尾礦庫(kù)穩(wěn)定性分析方法。極限平衡分析方法具有原理簡(jiǎn)單、物理意義明確等優(yōu)點(diǎn)，是工程中應(yīng)用最為廣泛的穩(wěn)定性分析方法。但是采用極限平衡法進(jìn)行尾礦庫(kù)穩(wěn)定性分析時(shí)需預(yù)先假設(shè)滑裂面形狀和初始滑裂面位置，且無(wú)法同時(shí)考慮滲流和變形等因素對(duì)尾礦庫(kù)穩(wěn)定性的影響，具有一定的局限性。

隨著計(jì)算機(jī)技術(shù)的發(fā)展，以強(qiáng)度折減法為代表的有限元穩(wěn)定性分析方法得到迅速發(fā)展。采用強(qiáng)度折減法進(jìn)行土體穩(wěn)定性分析時(shí)，不需要預(yù)先假設(shè)滑裂面的形狀和位置，可以同時(shí)獲得整個(gè)土體的應(yīng)力場(chǎng)和變形場(chǎng)，彌補(bǔ)了極限平衡法的不足[7]。該方法的在尾礦庫(kù)穩(wěn)定性分析中應(yīng)用的關(guān)鍵，是如何合理地選擇失穩(wěn)破壞判別準(zhǔn)則。目前土體穩(wěn)定性分析中使用最為廣泛的失穩(wěn)破壞準(zhǔn)則有“數(shù)值收斂準(zhǔn)則”、“位移突變準(zhǔn)則”和“滑裂面貫通準(zhǔn)則”，各準(zhǔn)則在尾礦庫(kù)失穩(wěn)判別中的適用性尚有待考證。

基于此，本文以國(guó)內(nèi)某上游法堆筑尾礦庫(kù)為研究對(duì)象，模擬了干灘長(zhǎng)度為100 m的尾礦庫(kù)滲流場(chǎng)及變形場(chǎng)的分布特性。在此基礎(chǔ)上，采用強(qiáng)度折減法分析了尾礦庫(kù)的穩(wěn)定性，評(píng)價(jià)了不同失穩(wěn)判別準(zhǔn)則的有效性和局限性，并提出了一種適用于上游法尾礦庫(kù)的失穩(wěn)判別準(zhǔn)則，以期為對(duì)滲流作用下尾礦庫(kù)的穩(wěn)定性分析提供技術(shù)支撐。

1 計(jì)算軟件與計(jì)算理論

1.1 ABAQUS軟件簡(jiǎn)介

本文采用通用有限元計(jì)算軟件ABAQUS進(jìn)行了尾礦庫(kù)滲流穩(wěn)定性分析。該軟件提供了非飽和土滲流應(yīng)力耦合的分析功能[8]，并包含十分豐富的材料模型、單元模式和荷載邊界條件，易于二次開(kāi)發(fā)，對(duì)巖土工程有較好地適用性[9]。

1.2 飽和-非飽和滲流理論

由于庫(kù)內(nèi)回水及降雨的影響，非干堆型尾礦庫(kù)可以看做不可壓縮的固相尾礦顆粒和液相水，以及一定數(shù)量的可壓縮氣體的三相體系。浸潤(rùn)面以下，尾砂被認(rèn)為是液—固兩相體系，遵循以比奧固結(jié)理論為基礎(chǔ)的飽和土滲流本構(gòu)方程；浸潤(rùn)面以上，非飽和尾砂為固—?dú)狻喝囿w系，其滲透本構(gòu)由滲透系數(shù)函數(shù)和飽和度(含水量)函數(shù)進(jìn)行表征[10]。

大型有限元軟件ABAQUS基于固定網(wǎng)格，按照非飽和土力學(xué)理論進(jìn)行求解，浸潤(rùn)面取為孔隙水壓力為零處，在進(jìn)行飽和-非飽和滲流時(shí)具有較好的計(jì)算精度和效率[11]。利用ABAQUS進(jìn)行不飽和滲流的材料設(shè)定時(shí)，滲透系數(shù)的關(guān)系可由一個(gè)與飽和度相關(guān)的折減系數(shù)進(jìn)行修正，如下式所示。

kunsat=ks·ksat

ks=(Sr)3Sr<1.0

ks=1.0Sr≥1.0

式中：kunsat為非飽和滲透系數(shù)，ks為折減系數(shù)，ksat為飽和滲透系數(shù)，Sr為土體飽和度。

在非飽和土體中，基質(zhì)吸力直接受材料飽和度的影響，二者關(guān)系可由吸濕(脫水)曲線表征，通常以基質(zhì)吸力(負(fù)孔壓)—飽和度的關(guān)系曲線給出。

1.3 強(qiáng)度折減理論

強(qiáng)度折減法是目前巖土行業(yè)公認(rèn)的可靠的、較為先進(jìn)的土坡穩(wěn)定分析計(jì)算方法。該方法采用抗剪強(qiáng)度折減系數(shù)Fr代替?zhèn)鹘y(tǒng)意義上的穩(wěn)定安全系數(shù)Fs，又稱為強(qiáng)度儲(chǔ)備安全系數(shù)。計(jì)算中假定不同的強(qiáng)度折減系數(shù)Fr，根據(jù)折減之后的強(qiáng)度參數(shù)進(jìn)行有限元分析，觀察計(jì)算是否收斂。在整個(gè)計(jì)算過(guò)程中不斷增加Fr，達(dá)到臨界破壞時(shí)的強(qiáng)度折減系數(shù)Fr就是邊坡穩(wěn)定安全系數(shù)Fs，折減后的抗剪強(qiáng)度參數(shù)可表達(dá)為：

cm=c/Fr

φm=arctan(tanφ/Fr)

式中：c和φ為土體抗剪強(qiáng)度，cm和φm為維持平衡所需的或土體實(shí)際發(fā)揮的抗剪強(qiáng)度；Fr為強(qiáng)度折減系數(shù)。

2 計(jì)算模型及材料參數(shù)

某尾礦庫(kù)采用初期壩+上游法堆壩，總壩高210 m。初期壩為碾壓堆石壩，壩高60 m，壩頂寬為4 m，上游邊坡坡比1∶1.6，下游邊坡坡比1∶1.8；堆積壩平均堆積邊坡為1∶5；兩側(cè)巖體為中風(fēng)化白云質(zhì)灰?guī)r，尾礦庫(kù)分區(qū)如圖 1所示。尾礦壩二維有限元模型被離散為3 153個(gè)節(jié)點(diǎn)及3 020個(gè)單元，采用具有孔壓自由度的四節(jié)點(diǎn)平面應(yīng)變單元和少量退化的三角形單元進(jìn)行模擬。

圖1 尾礦庫(kù)分區(qū)及有限元網(wǎng)格示意圖Fig.1 Zoning and finite element mesh of tailing reservoir

計(jì)算中有限元的網(wǎng)格固定于土骨架，氣體或液體在流過(guò)網(wǎng)格的過(guò)程中滿足流體連續(xù)方程。土體采用基于有效應(yīng)力指標(biāo)的摩爾庫(kù)倫彈塑性模型模擬，液體滲流采用Forchheimer滲透定律模擬，參數(shù)見(jiàn)表1及圖2，所有材料泊松比取為0.3。需要注意的是，一般巖土工程勘察報(bào)告以及《尾礦設(shè)施設(shè)計(jì)規(guī)范》主要提供土體的壓縮模量Es，而壓縮模量難以直接應(yīng)用于有限元數(shù)值計(jì)算中。賈堤[12]等人探討了利用場(chǎng)地的巖土工程勘察報(bào)告數(shù)據(jù)獲得數(shù)值分析中土體彈性模量的三種方法，建議采用公式E=αEs來(lái)估算彈性模量E，其中α建議取為8.2，Es為土體壓縮模量，本文采用該方法計(jì)算尾砂材料的彈性模量。

表1 計(jì)算參數(shù)表Table 1 Calculation parameters of various materials

圖2 土水特征曲線Fig.2 Characteristic curves of soil water

3 非飽和滲流條件下尾礦庫(kù)流固耦合分析

浸潤(rùn)線的位置是決定尾礦庫(kù)滲流情況及穩(wěn)定性的關(guān)鍵因素[13-15]。本文考慮了飽和度及基質(zhì)吸力對(duì)尾砂滲透系數(shù)的影響，圖3給出了干灘長(zhǎng)度分別為100 m時(shí)尾礦庫(kù)浸潤(rùn)面位置分布情況，即總水頭為0的水壓力等值線分布情況。可以看出，尾礦庫(kù)浸潤(rùn)線基本沿下游壩坡面分布，在壩中部位置略有抬高，最小埋深5 m。由于初期壩為透水堆石壩，因此浸潤(rùn)面均經(jīng)過(guò)初期壩上游坡面流向下游壩外。

圖3 干灘長(zhǎng)度為100 m時(shí)尾礦庫(kù)浸潤(rùn)線分布Fig.3 Distribution of phreatic line in tailings reservoir with the beach width of 100 m

圖4和圖5分別給出滲流條件下尾礦庫(kù)變形場(chǎng)分布及變形矢量圖?？梢钥闯觯畲笏较蜃冃沃饕l(fā)生在庫(kù)區(qū)1/3高度的下游坡面處，最大豎向變形主要發(fā)生在庫(kù)頂灘面處。尾砂的變形矢量主要指向庫(kù)底，即庫(kù)區(qū)內(nèi)尾砂主要發(fā)生向下的豎向變形；上游坡面處尾砂略微出現(xiàn)指向上游的水平變形。

圖4 尾礦庫(kù)變形分布圖Fig.4 Distribution of deformations in tailings reservoir

圖5 尾礦庫(kù)變形矢量圖Fig.5 Deformation vector diagram of tailings reservoir

4 尾礦庫(kù)滲流穩(wěn)定性分析

在采用強(qiáng)度折減法對(duì)尾礦庫(kù)進(jìn)行穩(wěn)定性分析時(shí)，強(qiáng)度折減過(guò)程尾礦庫(kù)失穩(wěn)判別準(zhǔn)則的選擇至關(guān)重要，目前常用的判別準(zhǔn)則有：1)“數(shù)值收斂準(zhǔn)則”：以數(shù)值計(jì)算收斂與否作為評(píng)價(jià)標(biāo)準(zhǔn)，其與有限元的算法有關(guān)；2)“位移突變準(zhǔn)則”：以特征部位的位移拐點(diǎn)作為評(píng)價(jià)標(biāo)準(zhǔn)；3)“滑裂面貫通準(zhǔn)則”：以是否形成連續(xù)的貫通區(qū)作為評(píng)價(jià)標(biāo)準(zhǔn)。其中，影響數(shù)值計(jì)算收斂與否的因素較多，如數(shù)值計(jì)算方法的穩(wěn)定性、計(jì)算模型網(wǎng)格質(zhì)量等，采用“數(shù)值收斂準(zhǔn)則”作為穩(wěn)定性判據(jù)可能帶來(lái)誤差[16]。因此本文主要采用“位移突變準(zhǔn)則”和“滑裂面貫通準(zhǔn)則”作為失穩(wěn)判別準(zhǔn)則，對(duì)干灘長(zhǎng)度為100 m的尾礦庫(kù)進(jìn)行了穩(wěn)定性評(píng)價(jià)，計(jì)算結(jié)果如圖6～7所示。

圖6給出尾礦庫(kù)外坡面坡頂和坡中特征點(diǎn)處位移與強(qiáng)度折減系數(shù)的關(guān)系。可以看出，坡頂特征點(diǎn)1的強(qiáng)度折減系數(shù)Fr由1.41增大到1.43的過(guò)程中，位移迅速增加而發(fā)生突變，該點(diǎn)對(duì)應(yīng)的安全系數(shù)取Fs為1.41。同理，坡中特征點(diǎn)2的安全系數(shù)Fs應(yīng)取為1.43?？梢钥闯?，采用“位移突變準(zhǔn)則”不同的特征點(diǎn)對(duì)應(yīng)的安全系數(shù)不同，難以選取某個(gè)單一的值作為整個(gè)尾礦庫(kù)的安全系數(shù)。另外，由于位移是一種標(biāo)量指標(biāo)，難以反映位移方向?qū)吰路€(wěn)定性的影響。

圖6 特征點(diǎn)位移與強(qiáng)度折減系數(shù)的關(guān)系Fig.6 Relationship of displacement and stability factor of characteristic points

圖7 尾礦庫(kù)滑裂面塑性應(yīng)變分布Fig.7 Plastical strain distribution of slipping surface of tailings reservoir

圖7給出了強(qiáng)度折減系數(shù)Fr分別為1.41和1.43時(shí)尾礦庫(kù)內(nèi)塑性應(yīng)變的分布情況。這兩種工況下庫(kù)區(qū)內(nèi)均出現(xiàn)了大范圍滑裂帶，滑裂面位置大致相同，貼近基巖面的深層滑動(dòng)。這是由于尾礦庫(kù)的分區(qū)導(dǎo)致各層尾砂的抗剪強(qiáng)度逐層減弱，強(qiáng)度最低的尾粉質(zhì)黏土主要沉積在庫(kù)底，形成抗剪薄弱面。需要注意的是，由于構(gòu)成初期壩的堆石體具有透水性好、強(qiáng)度高的特點(diǎn)，因此塑性應(yīng)變區(qū)無(wú)法通過(guò)初期壩而形成貫通面。即由于初期壩的存在，在采用強(qiáng)度折減法進(jìn)行上游法尾礦庫(kù)穩(wěn)定性分析的時(shí)候，采用“滑裂面貫通準(zhǔn)則”難以對(duì)尾礦庫(kù)是否失穩(wěn)進(jìn)行準(zhǔn)確判別。

本文提出一種新的尾礦庫(kù)失穩(wěn)評(píng)價(jià)準(zhǔn)則，即“變形矢量圖準(zhǔn)則”，該方法以初期壩附近尾砂出現(xiàn)指向壩外、及沿著初期壩上游坡向上的變形矢量為整個(gè)尾礦庫(kù)失穩(wěn)的判別依據(jù)。圖8給出Fr=1.41、Fr=1.43、Fr=1.53工況下尾礦庫(kù)整體變形矢量圖。Fr=1.41時(shí)，尾礦靠近外坡附近堆積尾砂主要產(chǎn)生指向庫(kù)底及庫(kù)外的變形，但變形矢量方向大部分未指出庫(kù)外，僅在初期壩附近呈現(xiàn)平行于堆積外坡面的變形；Fr=1.43時(shí)，庫(kù)區(qū)中下部尾砂整體出現(xiàn)指向庫(kù)外的變形，由于初期壩的阻礙作用，靠近初期壩的尾砂產(chǎn)生沿著初期壩坡向上的變形，尾礦庫(kù)發(fā)生失穩(wěn)破壞。隨著Fr增大為1.53，失穩(wěn)破壞繼續(xù)發(fā)展。因此，該尾礦庫(kù)的安全系數(shù)Fs取為1.41。通過(guò)上述分析可知，“變形矢量圖準(zhǔn)則”既可以直觀地呈現(xiàn)整個(gè)尾礦庫(kù)的變形情況，又可以給出滑動(dòng)體的分布情況，從而得到整體尾礦庫(kù)的穩(wěn)定性綜合安全系數(shù)，克服了單獨(dú)采用“位移突變準(zhǔn)則”和“滑裂面貫通準(zhǔn)則”的不足。

圖8 不同強(qiáng)度折減系數(shù)Fr下尾礦庫(kù)變形矢量圖Fig.8 Deformation vectors of tailings reservoir with various Fr

5 結(jié)論

本文采用強(qiáng)度折減法對(duì)某上游法尾礦庫(kù)飽和—非飽和滲流條件下穩(wěn)定性進(jìn)行了系統(tǒng)的分析，對(duì)比了不同失穩(wěn)判別準(zhǔn)則在尾礦庫(kù)穩(wěn)定性分析中的適用性和局限性，在此基礎(chǔ)上提出了一種新的判別準(zhǔn)則。主要結(jié)論如下：

1)本文基于飽和—不飽和滲流理論，對(duì)尾礦庫(kù)滲流條件下的浸潤(rùn)線分布和變形特性進(jìn)行了分析。研究結(jié)果表明，滲流作用下尾礦庫(kù)豎直方向主要發(fā)生豎向沉降，隨著高度的增加豎向沉降增大；水平向主要發(fā)生指向庫(kù)內(nèi)的變形，且在外坡面1/3處達(dá)到最大。

2)尾砂在庫(kù)區(qū)內(nèi)呈現(xiàn)強(qiáng)度較高的粗顆粒先沉積，強(qiáng)度較低的細(xì)顆粒后沉積的分區(qū)特性，因此越靠近庫(kù)底及庫(kù)內(nèi)岸坡尾砂的抗剪強(qiáng)度指標(biāo)越低，尾礦庫(kù)一旦失穩(wěn)，將發(fā)生沿著庫(kù)底岸坡的深層失穩(wěn)破壞。

3)本文建議采用“變形矢量圖準(zhǔn)則”作為上游法尾礦庫(kù)的失穩(wěn)判別準(zhǔn)則。該方法以初期壩附近尾砂出現(xiàn)指向壩外、及沿著初期壩上游坡向上的變形矢量為標(biāo)志，可綜合地判斷尾礦庫(kù)的整體穩(wěn)定安全系數(shù)。