基于數(shù)值模擬的露天礦高陡邊坡可靠度分析

尹 歸,陶干強(qiáng),朱忠華*,孫海帝

(1.南華大學(xué) 資源與環(huán)境安全工程學(xué)院,湖南 衡陽(yáng) 421001;2.湖南柿竹園有色金屬有限公司,湖南 郴州 423037)

0 引 言

隨著我國(guó)露天開(kāi)采逐步向深部發(fā)展，高陡邊坡穩(wěn)定性問(wèn)題愈發(fā)嚴(yán)重，其穩(wěn)定性已成為露天礦安全生產(chǎn)的關(guān)鍵技術(shù)問(wèn)題之一[1]，也是邊坡處治措施及安全評(píng)價(jià)的研究熱點(diǎn)。數(shù)學(xué)分析方法是目前學(xué)界與工程界進(jìn)行邊坡穩(wěn)定性評(píng)價(jià)的常用方法，曹蘭柱等[2]、韓萬(wàn)東等[3]分別采用Morgenstern-Price極限平衡法以及Spencer條分法和Bishop條分法對(duì)邊坡破壞機(jī)制和穩(wěn)定性狀況進(jìn)行了計(jì)算分析。但由于巖土體是一種具有非連續(xù)性、非均質(zhì)性的復(fù)雜材料，巖土體力學(xué)特性離散性較為顯著，且受限于技術(shù)水平，所獲取的力學(xué)參數(shù)具有較大的不確定性。這使得采用確定力學(xué)參數(shù)的常規(guī)數(shù)學(xué)分析方法的計(jì)算結(jié)果的可靠性與準(zhǔn)確性受到影響，而概率分析方法的引入，則為該問(wèn)題的解決提供了新途徑。蔣水華等[4]應(yīng)用隨機(jī)響應(yīng)面法對(duì)邊坡系統(tǒng)展開(kāi)了可靠度分析；唐小松等[5]則利用Bootstrap抽樣技術(shù)對(duì)邊坡可靠度展開(kāi)研究；彭興等[6]基于蒙特卡洛模擬實(shí)驗(yàn)理論，為邊坡可靠性分析提供了方法和框架。隨著研究逐步進(jìn)展，邊坡穩(wěn)定性分析已經(jīng)有了大量研究成果，計(jì)算體系也逐漸完整，但目前對(duì)邊坡展開(kāi)研究時(shí)，常局限于采用單一的方法進(jìn)行研究，若將兩種或多種分析方法綜合起來(lái)對(duì)邊坡穩(wěn)定性進(jìn)行分析，會(huì)更加有效。

邊坡穩(wěn)定性分析是分析邊坡穩(wěn)定程度及合理確定邊坡參數(shù)的重要方法，為邊坡處治和安全評(píng)價(jià)提供依據(jù)[7]。目前，邊坡穩(wěn)定性分析常用的方法有極限平衡分析方法、數(shù)值分析方法和可靠度分析方法。這些分析方法由于需要特定的假設(shè)條件，都具有一定的局限性。數(shù)值模擬方法能夠得到邊坡最可能的滑動(dòng)面和相應(yīng)的安全系數(shù)，但沒(méi)有統(tǒng)一的邊坡失穩(wěn)判據(jù)，不能判斷邊坡是否失穩(wěn)，無(wú)法得到邊坡的破壞概率；可靠度分析方法將邊坡部分參數(shù)視為隨機(jī)變量，用概率論方法分析邊坡穩(wěn)定性，可以得到邊坡破壞概率，但是計(jì)算過(guò)程需多次假定滑動(dòng)面，才能找出其最可能的滑動(dòng)面，耗時(shí)較長(zhǎng)。因此，本文將可靠度理論與數(shù)值模擬方法相結(jié)合進(jìn)行邊坡穩(wěn)定性研究，并以某露天礦高陡邊坡為算例，驗(yàn)證該方法的可行性。首先，根據(jù)露天礦邊坡的工程地質(zhì)特征，確定邊坡模型的材料屬性及材料的物理屬性，并以工程實(shí)際尺寸，使用三維數(shù)字礦山軟件建立邊坡三維表面模型，再利用C++編程實(shí)現(xiàn)限定Delaunay四面體網(wǎng)格剖分。剖分結(jié)束后，將剖分的網(wǎng)格數(shù)據(jù)轉(zhuǎn)化為后綴為flac3d或f3grid的計(jì)算模型，通過(guò)命令impgrid載入FLAC3D[8]，采用強(qiáng)度折減法模擬計(jì)算得到最有可能的滑動(dòng)面和相應(yīng)的安全系數(shù)；最后用驗(yàn)算點(diǎn)法對(duì)其最有可能的滑動(dòng)面進(jìn)行可靠度計(jì)算，確定可靠性指標(biāo)和破壞概率。

1 邊坡數(shù)值模擬

1.1 四面體網(wǎng)格剖分

網(wǎng)格剖分是有限元分析前處理的關(guān)鍵技術(shù)，四面體網(wǎng)格相比于六面體網(wǎng)格而言，具有更好的幾何適應(yīng)性，而且靈活性好、結(jié)構(gòu)簡(jiǎn)單、自動(dòng)化程度高。對(duì)模型進(jìn)行剖分時(shí)，只要生成的四面體數(shù)量足夠多，就能無(wú)限接近于模型邊界，能夠很好地模擬復(fù)雜的地質(zhì)體，不容易引起應(yīng)力集中，且計(jì)算模擬的精度有所提高。因此四面體網(wǎng)格在三維空間得到廣泛應(yīng)用，是當(dāng)今網(wǎng)格生成領(lǐng)域的研究熱點(diǎn)。Delaunay四面體網(wǎng)格剖分是目前最為流行的剖分方法之一。

限定Delaunay四面體剖分就是把一個(gè)帶邊界約束三維區(qū)域空間劃分成一個(gè)四面體網(wǎng)格，對(duì)三維模型進(jìn)行限定Delaunay四面體剖分，先將模型劃分為若干個(gè)限定線(xiàn)段和限定平面，而且在最后生成的四面體網(wǎng)格中這些限定線(xiàn)段和限定平面必須存在；并將限定的線(xiàn)段和平面分別分成小線(xiàn)段和小三角形，再對(duì)這些小線(xiàn)段和小三角形的頂點(diǎn)組成的點(diǎn)集進(jìn)行Delaunay四面體剖分，限定Delaunay四面體則就是由這個(gè)點(diǎn)集的頂點(diǎn)構(gòu)成[9]。

一般情況下，頂點(diǎn)集中的任意兩點(diǎn)在三維區(qū)域都可見(jiàn)，若兩點(diǎn)之間存在約束面，且兩點(diǎn)連線(xiàn)與該面的夾角不為0，則認(rèn)為兩頂點(diǎn)不可見(jiàn)，可見(jiàn)性受到約束面片的阻擋。如圖1所示，點(diǎn)p和點(diǎn)q是區(qū)域中的兩個(gè)頂點(diǎn)，位于面片f的兩側(cè)，他們的連線(xiàn)與f相交，p、q不可見(jiàn)。

圖1 限定Delaunay四面體規(guī)則示意圖

一個(gè)限定Delaunay四面體的四個(gè)頂點(diǎn)都是相互可見(jiàn)的，其外接球不包含除四面體頂點(diǎn)外的其他頂點(diǎn)。如圖1所示，點(diǎn)a、b、c、p構(gòu)成四面體是限定Delaunay四面體。

本文采用C++語(yǔ)言編寫(xiě)四面體剖分程序，剖分算法數(shù)據(jù)結(jié)構(gòu)主要包括點(diǎn)、面和四面體數(shù)據(jù)結(jié)構(gòu)。點(diǎn)數(shù)據(jù)結(jié)構(gòu)主要記錄剖分后生成的數(shù)據(jù)點(diǎn)的數(shù)量和各個(gè)數(shù)據(jù)點(diǎn)的三維坐標(biāo)；面數(shù)據(jù)結(jié)構(gòu)主要記錄輸入模型的邊界面?zhèn)€數(shù)和邊界面的頂點(diǎn)信息，以及四面體剖分后生成三角面?zhèn)€數(shù)和每個(gè)三角面頂點(diǎn)信息；四面體結(jié)構(gòu)記錄剖分后的四面體信息，包括生成的四面體個(gè)數(shù)和每個(gè)四面體的4個(gè)頂點(diǎn)的信息。算法首先生成包含整個(gè)輸入模型的大四面體，然后對(duì)模型限定邊界進(jìn)行離散，將每條限定邊分成若干個(gè)線(xiàn)段的集合，在限定邊離散的基礎(chǔ)上對(duì)限定面進(jìn)行平面域上的三角化，生成所有的限定點(diǎn)集(包括限定點(diǎn)，限定線(xiàn)段的端點(diǎn)和三角網(wǎng)格的頂點(diǎn))的Delaunay四面體網(wǎng)格，使其滿(mǎn)足平面的限定條件。然后對(duì)限定邊和限定面進(jìn)行恢復(fù)，如果某條線(xiàn)段在四面體集合中沒(méi)有出現(xiàn)，則加入該線(xiàn)段中點(diǎn)細(xì)分，直到受限邊出現(xiàn)在剖分結(jié)果中，若在限定面上生成的三角形被干涉，則將該三角形的外心點(diǎn)加入到四面體網(wǎng)格中，從而達(dá)到恢復(fù)受限邊和面的目的。最后判斷生成的網(wǎng)格所有四面體單元是否滿(mǎn)足Delaunay外接球準(zhǔn)則，若不滿(mǎn)足，則對(duì)網(wǎng)格進(jìn)行優(yōu)化調(diào)整，使其滿(mǎn)足Delaunay外接球準(zhǔn)則。

1.2 強(qiáng)度折減法基本原理

強(qiáng)度折減法將邊坡的安全系數(shù)定義為使邊坡剛好達(dá)到臨界破壞狀態(tài)時(shí)，對(duì)其強(qiáng)度參數(shù)進(jìn)行折減的程度[10]。是將確定的邊坡原始黏結(jié)力C和內(nèi)摩擦角φ的正切值同時(shí)除以折減系數(shù)K[11-14]：

(1)

通過(guò)不斷試算折減系數(shù)K，最后得到使邊坡產(chǎn)生破壞時(shí)的折減系數(shù)K就是邊坡的安全系數(shù)。強(qiáng)度折減法采用二分法迭代計(jì)算折減系數(shù)，使用二分法迭代折減系數(shù)就是不斷將折減系數(shù)下限值Kmin和上限值Kmax之間的范圍縮小到前一次迭代時(shí)的一半，當(dāng)折減系數(shù)上限值Kmax與下限值Kmin小于給定計(jì)算精度η時(shí)，停止計(jì)算。根據(jù)上述折減系數(shù)二分迭代算法，折減系數(shù)迭代次數(shù)Imin可用式(2)表示：

(2)

式中：Kmax、Kmin分別為折減系數(shù)上限值和下限值；η為給定誤差。

2 可靠度計(jì)算方法

邊坡可靠性分析方法是以定值安全系數(shù)法為基礎(chǔ)，根據(jù)極限平衡原理建立狀態(tài)方程，用以計(jì)算邊坡可靠度概率的分析方法[15]?？煽恐笜?biāo)法包括中心點(diǎn)法和驗(yàn)算點(diǎn)法，但中心點(diǎn)法存在難以對(duì)隨機(jī)變量實(shí)際分布加以考慮等缺點(diǎn)[16-17]，故采用驗(yàn)算點(diǎn)法進(jìn)行計(jì)算。

影響邊坡穩(wěn)定因素有很多，如降雨、地震或爆破作用、邊坡巖體性質(zhì)、是否存在結(jié)構(gòu)面以及結(jié)構(gòu)面大小等，這些不確定的變量就是隨機(jī)變量，在邊坡可靠性計(jì)算中用于構(gòu)建描述邊坡?tīng)顟B(tài)的函數(shù)模型：

Z=g(X)=g(X1,X2,…,Xn)

(3)

當(dāng)功能函數(shù)式(3)比較復(fù)雜時(shí)，為近似求解功能函數(shù)的均值與方差，求解結(jié)構(gòu)可靠度指標(biāo)，需引入標(biāo)準(zhǔn)化變量，若Xi為某一變量的特定值，對(duì)應(yīng)的標(biāo)準(zhǔn)化變量xi為：

(4)

其中μXi和σXi分別為Xi的均值和標(biāo)準(zhǔn)差。

標(biāo)準(zhǔn)化變量的平均值點(diǎn)是標(biāo)準(zhǔn)空間中的坐標(biāo)原點(diǎn)。因此邊坡的狀態(tài)方程可表示為：

Z=g(x)=g(x1,x2,…,xn)

(5)

(6)

Z的均值和標(biāo)準(zhǔn)差為：

(7)

(8)

則對(duì)于所有i值，用標(biāo)準(zhǔn)化變量表示解為：

x*=-αiβ

(9)

從原點(diǎn)到x*的距離即是可靠指標(biāo)。即：

(10)

由于驗(yàn)算點(diǎn)和β皆為未知，需用迭代法計(jì)算，采用Fiessler迭代計(jì)算步驟如下：

1)建立狀態(tài)方程Z=g(x)；

2)導(dǎo)出g(x)的表達(dá)式；

4)令x=0，β=0；

6)計(jì)算g(x)；

7)由下式算出Z的標(biāo)準(zhǔn)差：

(11)

8)用下式計(jì)算新的x值：

(12)

9)按式(10)計(jì)算β值；

重復(fù)步驟5)到9)直到數(shù)值收斂到允許誤差內(nèi)為止。

迭代計(jì)算收斂后，得到可靠性指標(biāo)β，進(jìn)而求得邊坡破壞概率。

3 工程應(yīng)用

3.1 地質(zhì)特征

本文研究的露天邊坡位于云南省易門(mén)銅廠，該露天礦屬于高原山區(qū)地形，中間被銅廠溝切割，海拔2 180～2 375 m，地形總體是中間低，兩邊高，地形坡度一般為15°～35°之間。礦體走向長(zhǎng)1 600 m，延伸大于300 m，礦體傾角在42°～45°之間，屬于傾斜礦體。境內(nèi)山巒重疊，地質(zhì)環(huán)境薄弱，邊坡穩(wěn)定性較差。邊坡位于該露天礦東南側(cè)，邊坡設(shè)計(jì)臺(tái)階高度為15 m，臺(tái)階坡面角為65°，總體邊坡高度為330 m，邊坡角度為45°。經(jīng)鉆孔揭露，構(gòu)成邊坡體的主要巖石有泥質(zhì)白云巖、灰白色白云巖、青灰色白云巖、火成膠結(jié)角礫巖等。邊坡所處位置地質(zhì)大致分為4個(gè)巖層，第一層至第四層分別為火成膠結(jié)角礫巖、青灰色白云巖、灰白色白云巖和泥質(zhì)白云巖，巖石節(jié)理裂隙不發(fā)育，水文地質(zhì)條件簡(jiǎn)單，地表水系不發(fā)育，地下水埋藏較深，因此本文未考慮地下水作用。

3.2 數(shù)值模擬

考慮到計(jì)算機(jī)計(jì)算能力和內(nèi)存大小的限制，建立塊段模型時(shí)，由于泥質(zhì)白云巖、灰白色白云巖、青灰色白云巖三種巖性力學(xué)參數(shù)十分相似，可合并成一種巖性。當(dāng)邊坡臺(tái)階數(shù)目較多時(shí)，采用FLAC3D對(duì)邊坡進(jìn)行穩(wěn)定性計(jì)算時(shí)，會(huì)出現(xiàn)單個(gè)臺(tái)階破壞或臺(tái)階聯(lián)合破壞，對(duì)邊坡整體破壞產(chǎn)生一定的影響，在研究邊坡的穩(wěn)定性時(shí)，本論文側(cè)重于研究邊坡的整體破壞，為避免臺(tái)階破壞對(duì)計(jì)算造成一定影響，將邊坡面轉(zhuǎn)化為平面，不考慮臺(tái)階的影響。計(jì)算中使用的力學(xué)參數(shù)如表1所示。另外，為配合數(shù)值計(jì)算的需要，需考慮邊坡巖體穩(wěn)定影響范圍和數(shù)值模擬軟件計(jì)算能力，根據(jù)地層確定的數(shù)值計(jì)算范圍見(jiàn)表2。依據(jù)數(shù)值計(jì)算范圍構(gòu)建邊坡三維模型，進(jìn)行限定Delaunay四面體剖分，將剖分結(jié)果導(dǎo)入到FLAC3D中，得到的數(shù)值模擬分析模型如圖2所示，模型朝采場(chǎng)方向的為白云巖組巖性模型，其余部分為火成膠結(jié)角礫巖組。另外，為便于計(jì)算結(jié)果對(duì)比，使用FLAC3D建立六面體網(wǎng)格如圖3所示。

圖2 邊坡數(shù)值模擬分析四面體模型

圖3 邊坡數(shù)值模擬分析六面體模型

表1 巖石物理力學(xué)參數(shù)

表2 模型尺寸參數(shù)表

為保持所建立的邊坡模型受力平衡，采用位移約束的方式將模型邊界條件設(shè)置為：前、后、左和右側(cè)邊界均為截離邊界，模型前后側(cè)邊界固定Y方向位移，模型左右側(cè)邊界固定X方向位移，模型底部固定Z方向位移，坡面為自由面。

采用強(qiáng)度折減法對(duì)邊坡模型展開(kāi)分析，邊坡剪切應(yīng)變?cè)隽吭茍D如圖4所示，圖中還包括邊坡的安全系數(shù)和速度矢量圖，速度矢量圖顯示出了邊坡產(chǎn)生滑動(dòng)時(shí)的滑動(dòng)趨勢(shì)。

從圖4中可以清楚地看到塑性貫通區(qū)域，通過(guò)速度矢量圖進(jìn)一步驗(yàn)證了這一判斷：處于滑動(dòng)面外側(cè)各網(wǎng)格點(diǎn)的速度均大于其他區(qū)域[18]，表明該邊坡已出現(xiàn)滑動(dòng)，邊坡穩(wěn)定結(jié)構(gòu)發(fā)生破壞，同時(shí)從圖4中可以得到此時(shí)相應(yīng)安全系數(shù)為1.80。圖5為六面體剖分的邊坡模型模擬結(jié)果，其安全系數(shù)為1.70，且在坡頂上方出現(xiàn)了應(yīng)力集中現(xiàn)象。與FLAC3D自帶的六面體網(wǎng)格剖分程序相比，四面體網(wǎng)格在進(jìn)行剖分時(shí)，網(wǎng)格單元數(shù)量多，雖然模擬過(guò)程中耗時(shí)比六面體網(wǎng)格長(zhǎng)，但計(jì)算精度相比于六面體而言計(jì)算精度會(huì)有所相應(yīng)的提高。而且四面體網(wǎng)格具有較為均勻的疏密程度，能夠很好地模擬復(fù)雜的地質(zhì)體，不容易引起應(yīng)力集中。

圖4 邊坡剪切應(yīng)變?cè)隽吭茍D及速度矢量圖

圖5 六面體網(wǎng)格剖分下邊坡剪切應(yīng)變?cè)隽吭茍D

圖6為露天礦邊坡強(qiáng)度折減后邊坡位移等值云圖和位移矢量圖。位移等值線(xiàn)圖形態(tài)表現(xiàn)為從邊坡頂部斜切進(jìn)入，并與坡面相交，等值線(xiàn)拐點(diǎn)處的切線(xiàn)近似與坡面平行，位移最大部分集中在邊坡中上部靠近坡面處，而且可以明顯看出邊坡上部靠近坡面位移矢量近乎與坡面平行，表現(xiàn)為“剪切”，中部和下部位移矢量在漸近坡趾處表現(xiàn)為“剪出”。這些現(xiàn)象都表明，邊坡潛在破壞以圓弧形剪切破壞為主。

圖6 邊坡位移等值云圖和位移矢量圖

3.3 可靠性分析

采用強(qiáng)度折減法使邊坡達(dá)到臨界狀態(tài)時(shí)，其潛在滑動(dòng)面附近等值線(xiàn)最為密集[19]。根據(jù)邊坡剖面位移等值云圖，確定等值線(xiàn)最為密集區(qū)域的位移值的范圍，利用FLAC3D程序內(nèi)置的編程語(yǔ)言FISH編制程序，搜索并保存等值線(xiàn)最為密集區(qū)域的數(shù)據(jù)點(diǎn)坐標(biāo)，利用Origin處理數(shù)據(jù)進(jìn)行擬合，得到的結(jié)果如圖7所示，其中各個(gè)離散點(diǎn)為等值線(xiàn)密集區(qū)域的點(diǎn)，曲線(xiàn)為該邊坡簡(jiǎn)化的潛在滑動(dòng)面。在此基礎(chǔ)上對(duì)該邊坡展開(kāi)可靠度計(jì)算。采用Fiessler迭代計(jì)算步驟編制邊坡可靠度計(jì)算程序，程序用C++語(yǔ)言編程實(shí)現(xiàn)。該露天礦邊坡剖面滑動(dòng)面示意圖見(jiàn)圖8，極限狀態(tài)函數(shù)為式(13)。

圖7 等值線(xiàn)密集區(qū)域離散點(diǎn)及其擬合線(xiàn)圖

圖8 邊坡計(jì)算剖面示意圖

Z=g(φ1,φ2)=W1×cosα×tanφ1+

W2×cosα×tanφ2+C2L2+C1L1-

(W1+W2)×sinα=0

(13)

式中：φ1、φ2—內(nèi)摩擦角，°；W1、W2—潛在滑動(dòng)面質(zhì)量，N/m3；C1、C2—黏聚力，N/m2；L1、L2—潛在滑動(dòng)面長(zhǎng)度，m；α—滑動(dòng)面傾角，°。

在諸多不確定因素中，物理力學(xué)參數(shù)對(duì)邊坡的安全穩(wěn)定性有著至關(guān)重要的影響，而且為了減少計(jì)算的復(fù)雜性，故選取對(duì)邊坡穩(wěn)定性影響較大的因素作為隨機(jī)變量，取φ1、φ2、α為隨機(jī)變量，該邊坡巖體相關(guān)力學(xué)參數(shù)如表3所示。

表3 邊坡巖體相關(guān)參數(shù)

采用可靠指標(biāo)法對(duì)該邊坡的破壞概率進(jìn)行迭代計(jì)算，迭代7次后，可靠度指標(biāo)β值達(dá)到收斂。求得可靠度指標(biāo)β=2.381 8，其中基本隨機(jī)變量X1為白云巖內(nèi)摩擦角，X2為火成巖內(nèi)摩擦角，X3為滑動(dòng)面傾角。根據(jù)Pf=φ(-β)=1-φ(β)，求得破壞概率Pf=0.865 6%。現(xiàn)有研究成果表明，露天礦總體邊坡可接受的破壞概率一般在0.3%～1.0%之間[20]。表明該邊坡在安全系數(shù)達(dá)到1.80的情況下，破壞概率Pf為0.865 6%，屬于穩(wěn)定邊坡，破壞概率較低。

4 結(jié) 論

1)本研究基于數(shù)值模擬平臺(tái)，建立了數(shù)值模擬和可靠度理論相結(jié)合的綜合分析方法，對(duì)邊坡穩(wěn)定性展開(kāi)分析。彌補(bǔ)了單一分析方法不足以有效評(píng)定邊坡穩(wěn)定性的缺陷。

2)對(duì)邊坡模型采用限定Delaunay四面體網(wǎng)格剖分，生成了高質(zhì)量的四面體網(wǎng)格模型。將四面體剖分與六面體剖分的數(shù)值模擬結(jié)果進(jìn)行對(duì)比分析，表明四面體剖分的優(yōu)越性。

3)采用FLAC3D內(nèi)置的FISH語(yǔ)言編程，在強(qiáng)度折減法的基礎(chǔ)上，自動(dòng)搜索保存滑動(dòng)面附近的數(shù)據(jù)點(diǎn)，具有簡(jiǎn)便實(shí)用和速度快的特點(diǎn)。

4)將該方法應(yīng)用于某露天礦邊坡穩(wěn)定性分析，得到該邊坡在安全系數(shù)達(dá)到1.80時(shí)，仍有0.865 6%的破壞風(fēng)險(xiǎn)，屬于穩(wěn)定邊坡。表明了該方法的可行性。