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

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

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

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

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

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

1.1 工程地質(zhì)

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

1.2 滑坡概況

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖4 最大不平衡力與計算步的關(guā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ù)開挖至2022年排產(chǎn)境界以及最終境界的三維塑性區(qū)分布圖(如圖6(a)、圖6(b))可以看出,潛在的塑性破壞區(qū)域主要分布在N14勘探線西幫附近和N26勘探線東幫附近,與滑坡現(xiàn)狀較為吻合,尤其是東幫N26線的塑性變形區(qū)域在2018年發(fā)生了較大的多臺階的順層滑坡,危及礦山的安全生產(chǎn)。對比東西幫塑性區(qū)分布(圖6(c)～圖6(f))表明:①東幫塑性區(qū)的分布范圍以及深度都遠(yuǎn)大于西幫,說明東幫的順傾邊坡穩(wěn)定性小于西幫;②隨著邊坡開挖深部的不斷更加,東幫的塑性區(qū)集中區(qū)域的分布范圍和深度增加,貫通分布在整個最終境界邊幫;西幫邊坡塑性區(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向(即東西向)的位移云圖可以直觀地反映在開挖過程中東西幫邊坡向臨空面方向的發(fā)展趨勢。圖7(a)～圖7(f)表明:①隨著開挖進(jìn)行,東西幫的邊坡向臨空向位移的位移量逐漸增加,最大X向位移絕對值分別由0.392 m和0.321 m增加至1.267 m和1.313 m;②東幫向臨空面的X向最大位移值大于西幫,變形區(qū)域遍布整個東幫,西幫邊坡的變形區(qū)域則主要集中在邊坡下部,這主要是由于東幫順傾和西幫反傾的巖體結(jié)構(gòu),導(dǎo)致東西幫變形具有各向異性特征,較大的變形亦不利于東幫邊坡的整體穩(wěn)定性。因此,為確保礦山的安全生產(chǎn),應(yīng)對東幫變形及破壞區(qū)采取適當(dāng)加固防治措施。

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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

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