司家營研山鐵礦采場三維邊坡穩定性數值模擬及協同監測分析

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

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

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

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

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

1 工程地質概況及邊坡滑坡現狀

1.1 工程地質

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

1.2 滑坡概況

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

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

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

2 三維地質模型建立及網格劃分

2.1 三維地質模型建立

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

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

2.2 數值網格劃分

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

圖3 有限元網格劃分Fig.3 Finite element mesh generation process

3 基于FLAC3D的三維邊坡穩定性計算

3.1 力學參數確定及邊界條件

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

表2 巖體物理力學參數指標Table 2 Physical and mechanical parameters of rock mass

3.2 平衡應力計算

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

圖4 最大不平衡力與計算步的關系曲線Fig.4 Relation curve between maximun unbalance force and calculating time step

圖5 現狀Z方向應力云圖剖面Fig.5 Stress in direction Z of mining status

3.3 塑性區分析

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

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

3.4 位移云圖分析

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

4 協同監測分析

4.1 協同監測介紹

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

4.2 數據分析

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

圖9 協同監測數據分析Fig.9 Analysis of collaborative monitoring data

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

5 結 論

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

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

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

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

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