楊意德 楊天鴻 劉 洋 鄧文學(xué) 李 華 田益琳 葉會(huì)師

(1.東北大學(xué)資源土木工程學(xué)院,遼寧 沈陽(yáng) 110819;2.河北鋼鐵集團(tuán)司家營(yíng)研山鐵礦有限公司,河北 唐山 063701)

露天礦邊坡穩(wěn)定性分析是邊坡安全設(shè)計(jì)的關(guān)鍵環(huán)節(jié)[1],就目前邊坡穩(wěn)定性分析來(lái)看,多數(shù)專(zhuān)家、學(xué)者分析二維模型較多,其優(yōu)點(diǎn)為可了解露天邊坡某個(gè)剖面在水位、爆破振動(dòng)及人工開(kāi)挖等影響下的穩(wěn)定性[2-4],分析多個(gè)剖面便可了解實(shí)際邊坡一定范圍內(nèi)的穩(wěn)定性狀況,給研究者提供不少便利。然而,邊坡巖層的走向和傾向影響邊坡空間應(yīng)力分布,尤其是在層狀邊坡中,邊坡應(yīng)力應(yīng)變分析簡(jiǎn)化為平面應(yīng)變問(wèn)題分析有明顯的不足之處,邊坡破壞產(chǎn)生的滑坡絕大多數(shù)是復(fù)雜的三維問(wèn)題。

在三維邊坡穩(wěn)定性分析方面,Duncan[5]曾從20篇文獻(xiàn)中歸納了三維邊坡穩(wěn)定性分析的特點(diǎn)和局限性,曹蘭柱[6]、宋子嶺[7]等對(duì)不同的三維邊坡分析方法作對(duì)比分析,前人研究表明三維有限元和三維極限平衡方法的應(yīng)用均有諸多前提條件[8-9],如此一來(lái),也就限制了三維穩(wěn)定性分析理論的應(yīng)用范圍,從而使三維邊坡穩(wěn)定性分析不能有效發(fā)展和應(yīng)用。相對(duì)極限平衡法,有限元強(qiáng)度折減法在分析邊坡穩(wěn)定性時(shí)具有的優(yōu)勢(shì)在其可自動(dòng)搜索滑面,同時(shí)滿(mǎn)足了力的平衡條件并且考慮了巖土體應(yīng)力—應(yīng)變關(guān)系、巖土體和加固結(jié)構(gòu)的耦合作用,更為嚴(yán)謹(jǐn)。有限元強(qiáng)度折減法在二維穩(wěn)定性分析中已經(jīng)取得大量研究成果并且已經(jīng)付諸于工程實(shí)踐當(dāng)中[10-14],其在三維邊坡穩(wěn)定性分析中也有很多應(yīng)用。在有限元分析中,網(wǎng)格劃分至關(guān)重要,網(wǎng)格質(zhì)量的好壞直接關(guān)系到分析是否能夠順利、快速地完成,也關(guān)系到是否能夠得到高精度的分析結(jié)果[15-17]。FLAC3D快速拉格朗日有限差分法方法是將差分方程近似表示微分方法的一種數(shù)值方法,由于其計(jì)算收斂性較好,在巖土邊坡工程中應(yīng)用廣泛,許多學(xué)者通過(guò)工程實(shí)踐進(jìn)行了研究[18-21],但是當(dāng)需要建立較復(fù)雜的三維模型時(shí),前處理往往為主要工作[22],但由于FLAC3D并無(wú)配套的三維前處理建模及網(wǎng)格劃分軟件,因此,需借助其他的工具建模并劃分網(wǎng)格并通過(guò)接口程序轉(zhuǎn)換為FLAC3D可使用的數(shù)據(jù)格式。

以司家營(yíng)研山鐵礦為研究對(duì)象,基于礦山地質(zhì)剖面建立三維地質(zhì)模型和地表DTM地表模型,并借助Rhino軟件及Griddle插件對(duì)復(fù)雜地表以及巖性界面進(jìn)行多軌掃略放樣、布簾簡(jiǎn)化、幾何布爾運(yùn)算以及有限元網(wǎng)格劃分等工作,將現(xiàn)狀邊坡分步開(kāi)挖至最終境界,進(jìn)行FLAC3D三維數(shù)值模擬計(jì)算,從宏觀角度分析邊坡的潛在塑性區(qū)分布及其演化規(guī)律后,再借助協(xié)同監(jiān)測(cè)系統(tǒng)對(duì)局部區(qū)域進(jìn)行了驗(yàn)證,并得到潛在滑區(qū)的進(jìn)一步認(rèn)識(shí),為邊坡優(yōu)化設(shè)計(jì)及礦山安全生產(chǎn)提供了依據(jù)。

1 工程地質(zhì)概況及邊坡滑坡現(xiàn)狀

1.1 工程地質(zhì)

司家營(yíng)研山鐵礦位于河北唐山灤州市,是國(guó)內(nèi)大型露天鐵礦之一。礦區(qū)被第四系地層大面積覆蓋,基巖露頭除在礦區(qū)東部和尚山～鐵石山一帶有較連續(xù)的分布外,其他均為零星出露。礦區(qū)采場(chǎng)內(nèi)中元古界長(zhǎng)城系大紅峪組(Pt2chd)石英砂巖作為沉積蓋層呈角度不整合覆蓋在新太古界單塔子巖群白廟子組黑云變粒巖(Ar3gnt)、混合巖化黑云變粒巖(Ar3mgnt)、鉀長(zhǎng)石化云母片巖(Ar3sch)等變質(zhì)巖系及太古代沉積變質(zhì)鐵礦(BIF)之上。除水平覆蓋于巖層上部的第四系之外,其他地層均為傾斜巖層,總體傾向西,其中大紅峪組石英砂巖傾角15°±5°,變粒巖系傾角40°±5°,構(gòu)造線(xiàn)近南北,形成了東幫順傾,西幫反傾的巖層結(jié)構(gòu)。根據(jù)目前開(kāi)采現(xiàn)狀,東邊坡由變粒巖系和第四系沖積層組成,西邊坡由石英砂巖和第四系沖積層組成,而南北兩邊坡則由其東半部的變粒巖系和西半部的石英砂巖構(gòu)成。

1.2 滑坡概況

近幾年司家營(yíng)研山鐵礦發(fā)生過(guò)多次規(guī)模不一的邊坡滑坡,滑坡的空間分布及部分滑坡現(xiàn)狀照片如圖1,其中規(guī)模較大的滑坡統(tǒng)計(jì)如表1所示。在經(jīng)過(guò)現(xiàn)有資料統(tǒng)計(jì)和現(xiàn)場(chǎng)勘查后發(fā)現(xiàn),司家營(yíng)研山鐵礦東幫主要由第四系土層、黑云變粒巖及白云母片巖組成,黑云變粒巖按照其風(fēng)化程度的不同從上到下又可以劃分成強(qiáng)風(fēng)化黑云變粒巖,中風(fēng)化黑云變粒巖及微風(fēng)化黑云變粒巖,且東幫臨近新河,白云母片巖遇水容易泥化,不同風(fēng)化程度的黑云變粒巖分界面傾向與邊坡傾向一致,極易發(fā)生順層滑坡破壞。西幫邊坡穩(wěn)定性?xún)?yōu)于東幫,主要滑坡分布在第四系表土層。

圖1 司家營(yíng)研山鐵礦歷史滑坡分布Fig.1 Historical landslide distribution map of Sijiaying Yanshan Iron Mine

表1 司家營(yíng)研山鐵礦露天邊坡歷史滑坡區(qū)域統(tǒng)計(jì)Table 1 Statistics on the historical landslide area of the open slope of Sijiaying Yanshan Iron Mine

2 三維地質(zhì)模型建立及網(wǎng)格劃分

2.1 三維地質(zhì)模型建立

基于礦山開(kāi)采現(xiàn)狀平面圖、勘探線(xiàn)剖面圖建立露天采場(chǎng)三維地質(zhì)模型。詳細(xì)的方法流程如圖2(a),具體步驟為:①基于3Dmine中的坐標(biāo)轉(zhuǎn)換將所有剖面分布于各對(duì)應(yīng)勘探線(xiàn),形成各剖面在三維空間上的分布狀態(tài),提取各巖層分界線(xiàn);②根據(jù)各巖層分界線(xiàn),利用Rhino建模軟件中NURBS曲面建模功能,通過(guò)多軌掃略放樣生成各不同巖層間分界面;③根據(jù)現(xiàn)狀、采剝計(jì)劃和最終境界平面圖,生成三角網(wǎng)DTM面,基于Rhino布簾工具建立簡(jiǎn)化的現(xiàn)狀面、2022年采剝計(jì)劃面以及最終境界設(shè)計(jì)平面(如圖2(b)、圖2(c)所示);④建立相應(yīng)尺度的立方體,基于Rhino布爾運(yùn)算的幾何切割功能進(jìn)行巖層及現(xiàn)狀境界等分界面的切割,最后形成司家營(yíng)鐵研山鐵礦三維地質(zhì)模型,如圖2(d)所示。

圖2 司家營(yíng)研山鐵礦地質(zhì)模型Fig.2 Geological model of Sijiaying Yanshan Iron Mine

2.2 數(shù)值網(wǎng)格劃分

復(fù)雜有限元網(wǎng)格劃分需保證界面之間網(wǎng)格節(jié)點(diǎn)的重合,即分界面兩側(cè)網(wǎng)格在該面上完全重合,此為三維復(fù)雜幾何網(wǎng)格建立的關(guān)鍵。由于Rhino的幾何分割存在一定的精度誤差,相鄰幾何分界面實(shí)際上無(wú)法完全重合,本研究首先劃分所有界面的表面網(wǎng)格(如圖3(a)),在確保表面網(wǎng)格連續(xù)、分界面上網(wǎng)格重合后(基于Rhino二次開(kāi)發(fā)的Griddle插件進(jìn)行網(wǎng)格劃分,可自動(dòng)對(duì)一定范圍內(nèi)的兩NURBS曲面上的網(wǎng)格進(jìn)行合并),進(jìn)一步劃分三維網(wǎng)格。司家營(yíng)研山鐵礦模型的整體尺寸為長(zhǎng)4 693m,寬為3 065m,高為1 066 m,采用四面體網(wǎng)格,共劃分單元數(shù)為16 065 000個(gè),節(jié)點(diǎn)數(shù)為2 830 020個(gè),導(dǎo)入FLAC3D中如圖3(b)所示,共計(jì)完成19個(gè)不同的復(fù)雜幾何網(wǎng)格劃分。

圖3 有限元網(wǎng)格劃分Fig.3 Finite element mesh generation process

3 基于FLAC3D的三維邊坡穩(wěn)定性計(jì)算

3.1 力學(xué)參數(shù)確定及邊界條件

模擬所需巖體力學(xué)參數(shù)主要來(lái)源于工程勘察報(bào)告以及相關(guān)的文獻(xiàn)中[23],本文不詳細(xì)贅述。主要巖體的物理力學(xué)參數(shù)指標(biāo)取值見(jiàn)表2所示。計(jì)算采用邊界條件為:對(duì)底部邊界采用位移約束,對(duì)四周采用法向位移約束。共分三步進(jìn)行平衡應(yīng)力計(jì)算,第一步采用彈性本構(gòu)模型計(jì)算現(xiàn)狀平衡應(yīng)力場(chǎng),第二步采用彈塑性本構(gòu)模型,在上一步基礎(chǔ)上開(kāi)挖至2022年的排產(chǎn)境界,計(jì)算其平衡應(yīng)力場(chǎng)并判斷塑性區(qū),第三步采用彈塑性本構(gòu)模型,并在第二步基礎(chǔ)上開(kāi)挖至最終的設(shè)計(jì)境界,進(jìn)行進(jìn)一步的計(jì)算分析。

表2 巖體物理力學(xué)參數(shù)指標(biāo)Table 2 Physical and mechanical parameters of rock mass

3.2 平衡應(yīng)力計(jì)算

邊坡穩(wěn)定性計(jì)算主要取決于最大不平衡力收斂以及不平衡比率情況,FLAC3D默認(rèn)當(dāng)不平衡比率達(dá)到10-5計(jì)算達(dá)到平衡并停止計(jì)算。圖4即為對(duì)現(xiàn)狀邊坡模型計(jì)算過(guò)程中的最大不平衡力與計(jì)算步的關(guān)系曲線(xiàn),該曲線(xiàn)變化趨于穩(wěn)定表明模型收斂較好,計(jì)算到一定的步數(shù)最大不平衡力便趨于零。圖5(a)和圖5(b)為現(xiàn)狀模型典型剖面的Z方向應(yīng)力云圖,邊坡垂直應(yīng)力在高程上呈現(xiàn)漸變趨勢(shì),底部應(yīng)力比頂部應(yīng)力大,呈現(xiàn)等值線(xiàn)分布,單斜構(gòu)造的巖層傾斜分布導(dǎo)致兩側(cè)應(yīng)力分布的不對(duì)稱(chēng),較符合現(xiàn)場(chǎng)實(shí)際狀況。

圖4 最大不平衡力與計(jì)算步的關(guān)系曲線(xiàn)Fig.4 Relation curve between maximun unbalance force and calculating time step

圖5 現(xiàn)狀Z方向應(yīng)力云圖剖面Fig.5 Stress in direction Z of mining status

3.3 塑性區(qū)分析

塑性區(qū)分布是判斷邊坡破壞的主要判據(jù)。根據(jù)開(kāi)挖至2022年排產(chǎn)境界以及最終境界的三維塑性區(qū)分布圖(如圖6(a)、圖6(b))可以看出,潛在的塑性破壞區(qū)域主要分布在N14勘探線(xiàn)西幫附近和N26勘探線(xiàn)東幫附近,與滑坡現(xiàn)狀較為吻合,尤其是東幫N26線(xiàn)的塑性變形區(qū)域在2018年發(fā)生了較大的多臺(tái)階的順層滑坡,危及礦山的安全生產(chǎn)。對(duì)比東西幫塑性區(qū)分布(圖6(c)～圖6(f))表明:①東幫塑性區(qū)的分布范圍以及深度都遠(yuǎn)大于西幫,說(shuō)明東幫的順傾邊坡穩(wěn)定性小于西幫;②隨著邊坡開(kāi)挖深部的不斷更加,東幫的塑性區(qū)集中區(qū)域的分布范圍和深度增加,貫通分布在整個(gè)最終境界邊幫;西幫邊坡塑性區(qū)域分布演化則由淺部第四系土層轉(zhuǎn)移至深部的基巖上,但分布區(qū)域范圍有限,因此西幫的整體邊坡安全系數(shù)高于東幫。

圖6 2022年排產(chǎn)境界及最終境界塑性區(qū)分布Fig.6 Distribution of plastic zone of 2022 production schedule pit limit and the ultimate pit limit

3.4 位移云圖分析

X向(即東西向)的位移云圖可以直觀地反映在開(kāi)挖過(guò)程中東西幫邊坡向臨空面方向的發(fā)展趨勢(shì)。圖7(a)～圖7(f)表明:①隨著開(kāi)挖進(jìn)行,東西幫的邊坡向臨空向位移的位移量逐漸增加,最大X向位移絕對(duì)值分別由0.392 m和0.321 m增加至1.267 m和1.313 m;②東幫向臨空面的X向最大位移值大于西幫,變形區(qū)域遍布整個(gè)東幫,西幫邊坡的變形區(qū)域則主要集中在邊坡下部,這主要是由于東幫順傾和西幫反傾的巖體結(jié)構(gòu),導(dǎo)致東西幫變形具有各向異性特征,較大的變形亦不利于東幫邊坡的整體穩(wěn)定性。因此,為確保礦山的安全生產(chǎn),應(yīng)對(duì)東幫變形及破壞區(qū)采取適當(dāng)加固防治措施。

4 協(xié)同監(jiān)測(cè)分析

4.1 協(xié)同監(jiān)測(cè)介紹

通過(guò)數(shù)值模擬結(jié)果及現(xiàn)場(chǎng)情況來(lái)看,研山鐵礦東幫有較大的滑坡風(fēng)險(xiǎn),因此對(duì)研山鐵礦東幫N26線(xiàn)進(jìn)行應(yīng)力計(jì)、微震協(xié)同監(jiān)測(cè),以進(jìn)一步評(píng)估東幫邊坡的穩(wěn)定性,監(jiān)測(cè)布置如圖8所示。微震監(jiān)測(cè)系統(tǒng)為6通道,在-157、-187 m各布置3個(gè),錨桿應(yīng)力計(jì)在-157、-187 m平臺(tái)各布置1個(gè)孔,每個(gè)孔分別布置4個(gè)應(yīng)力計(jì),每2個(gè)應(yīng)力計(jì)之間間隔為4 m。

4.2 數(shù)據(jù)分析

對(duì)研山鐵礦東幫順層邊坡進(jìn)行了2021-04-29—2121-08-27的協(xié)同監(jiān)測(cè),微震監(jiān)測(cè)主要針對(duì)的是開(kāi)挖工作對(duì)圍巖的擾動(dòng)情況,每次的爆破開(kāi)挖都會(huì)對(duì)邊坡圍巖造成一定程度的損傷,而每次擾動(dòng)都會(huì)有對(duì)應(yīng)的微震信號(hào),如圖9(a)所示,不同灰度代表能量的大小,在-127 m以下有大量的微震監(jiān)測(cè)信號(hào),而淺部的微震信號(hào)多于深部,說(shuō)明淺部的損傷大于深部,且在坡體內(nèi)6～7m深度處出現(xiàn)較為集中的微震信號(hào)貫通,說(shuō)明爆破振動(dòng)在邊坡表層形成了一定深度的損傷區(qū)。

圖9 協(xié)同監(jiān)測(cè)數(shù)據(jù)分析Fig.9 Analysis of collaborative monitoring data

如圖9(c)所示,1～4號(hào)(-187 m)應(yīng)力計(jì)中整體應(yīng)力變化分別為 3.8、2.5、0.1、2.8 kN,降雨(7月13日)時(shí)最大應(yīng)力變化分別為 4.7、1.9、1.6、1.3 kN;如圖9(d)所示,5～8號(hào)(-157 m)應(yīng)力計(jì)中整體應(yīng)力變化分別為 5.6、6.9、1.1、0.8 kN,降雨時(shí)最大應(yīng)力變化 1.3、1.9、0.6、0.6 kN。 東幫布置的8個(gè)應(yīng)力計(jì)均有應(yīng)力變化,尤其是在降雨期間傳感器應(yīng)力有明顯變化,在雨后隨著坡內(nèi)靜水壓力下降后,應(yīng)力又逐漸恢復(fù)原水平,且上下2組應(yīng)力計(jì)均體現(xiàn)邊坡淺部應(yīng)力大于深部。綜上所述,微震與應(yīng)力計(jì)監(jiān)測(cè)到的結(jié)果相對(duì)應(yīng),淺部的損傷與應(yīng)力變化大于深部,爆破振動(dòng)對(duì)邊坡表層形成損傷區(qū),且降雨對(duì)于邊坡內(nèi)部應(yīng)力變化作用明顯。

5 結(jié) 論

在收集司家營(yíng)研山鐵礦的大量地質(zhì)剖面資料基礎(chǔ)上,結(jié)合礦山采剝計(jì)劃和最終境界設(shè)計(jì)方案,基于3Dmine和Rhino軟件,完成了礦山三維地質(zhì)模型的建立,并劃分出FLAC3D可直接使用的有限元網(wǎng)格以作東西幫整體采場(chǎng)的三維穩(wěn)定性分析,分兩步開(kāi)挖進(jìn)行數(shù)值模擬計(jì)算,并對(duì)東幫局部進(jìn)行協(xié)同監(jiān)測(cè)驗(yàn)證及進(jìn)一步分析,得到以下結(jié)論:

(1)利用 Rhino的 NURBS曲面建模功能及其Griddle插件的表面網(wǎng)格重合處理功能,可較好地實(shí)現(xiàn)露天礦復(fù)雜地質(zhì)體的三維建模及其有限元網(wǎng)格劃分,網(wǎng)格質(zhì)量可保證計(jì)算收斂性。

(2)潛在的塑性及變形區(qū)域,在東西幫均有分布,這主要是由單斜構(gòu)造影響下形成了邊坡順層與反傾2種截然不同的結(jié)構(gòu),應(yīng)力、變形及潛在破壞范圍均具有各向異性特征,主要分布于東幫邊坡,而西幫邊坡的整體穩(wěn)定性高于東幫。

(3)通過(guò)微震、應(yīng)力計(jì)協(xié)同監(jiān)測(cè)分析可知,研山東幫N26線(xiàn)有較明顯的損傷與應(yīng)力變化,與上述數(shù)值模擬結(jié)果吻合,且降雨對(duì)于邊坡內(nèi)部應(yīng)力變化作用明顯,淺部的應(yīng)力與損傷變化大于深部,應(yīng)采取相應(yīng)治理措施,尤其是淺部,以確保東幫邊坡的長(zhǎng)期穩(wěn)定性。

(4)本研究存在不足之處,未充分考慮層狀巖體原生各向異性的變形及強(qiáng)度差異,當(dāng)前同一類(lèi)型巖體力學(xué)參數(shù)按同一參數(shù)取值,因此,對(duì)采場(chǎng)整體三維邊坡穩(wěn)定性分析方法仍具有進(jìn)一步改進(jìn)的空間。