黃 敏,崔年生,危劍林,何建華,孔繁成,黃英華
(1.長沙礦山研究院有限責任公司,湖南 長沙 410012;2.金屬礦山安全技術國家重點實驗室,湖南 長沙 410012;3.福建省新華都工程有限責任公司,福建 廈門 361000;4.湖南有色新田嶺鎢業(yè)有限公司,湖南 郴州 423000;5.紹興銅都礦業(yè)有限公司,浙江 紹興 312000)
開展金屬礦山地下開采引起的地表巖層移動變形規(guī)律研究,對評判地表建(構)筑物安全程度、圈定地表影響范圍、緩解礦山與周邊居民日益緊張的關系、征用土地、搬遷房屋以及保證礦山安全生產(chǎn)具有十分重要的意義[1-4]。當前,眾多科研工作者都對其進行了探索和分析,取得了許多有價值的科研成果,如王艷輝等[5]采用正交試驗及數(shù)值模擬計算方法對不同采高、不同冒落高度、不同巖體力學參數(shù)下的巖層移動規(guī)律進行了深入研究,具有一定的借鑒作用;夏開宗等[6]通過對具有陡傾結構面的程潮鐵礦西區(qū)的地表變形監(jiān)測資料及宏觀破壞特征進行了分析,得出礦區(qū)巖層移動分為采空區(qū)頂板巖體破壞擴展至地表引起的塌陷階段及采空區(qū)周邊圍巖向采空區(qū)的傾倒破壞階段,并得出了巖體主要發(fā)生水平移動的變形規(guī)律,呈現(xiàn)出先緩慢變形,后快速變形的特征,兩者之間有明顯的轉(zhuǎn)折點;趙海軍等[7]在自重應力和構造應力兩種工況下,采用數(shù)值模擬方法對急傾斜礦體開采的巖移特征與變形機理進行了詳細分析,得出了構造應力條件下地表出現(xiàn)雙沉降中心的,自重應力條件下只存在單沉降中心的現(xiàn)象,同時,兩者在地表移動變形量、移動變形影響區(qū)規(guī)模及地表宏觀變形破壞特征上存在較大差異;李文秀等[8]采用FLAC3D軟件對金寶山鉑鈀礦地下采礦過程引起的地表變形和巖層移動變形規(guī)律進行了三維仿真模擬,得出了地下開采引起的地表巖層移動是一個連續(xù)變化的過程并預測了地表移動范圍;袁仁茂等[9]采用理論力學及數(shù)值計算的方法對金川二礦區(qū)急傾斜厚大礦體地下開采引起的地表巖層移動傳遞模式及地表巖層移動規(guī)律進行了深入研究,并結合地表監(jiān)測結果,得出礦體厚度對地表巖移的顯現(xiàn)特征具有明顯的影響;廉海等[10]運用FLAC軟件對石人溝鐵礦地下開采過程進行分步開挖模擬,得到了不同開挖階段的地表動態(tài)移動曲線;陸清運等[11]等先利用地表監(jiān)測數(shù)據(jù)與三維數(shù)值模擬計算反演出材料力學參數(shù),后采用三維數(shù)值方法分析地下開采冒通地表與未冒通地表兩種工況下的地表變形特征,得到了巖層移動角的變化趨勢。
礦山地下開采對地表建(構)筑物的影響主要取決于地質(zhì)和采礦等因素,如礦床的開采深度、采高、傾角、采礦工藝、地質(zhì)結構、圍巖特性等。為探究地下開采的影響范圍,評估地表建構物安全程度,指導礦山安全生產(chǎn),本文采用FLAC3D有限元差分法軟件分析礦山開采后地表巖層移動變形情況,評判礦山地下開采對地表建(構)筑物的影響范圍和程度。
哈圖金礦屬中型地下金礦山,采用地下開采方式,設計生產(chǎn)能力為1 200 t/d。齊Ⅰ金礦區(qū)地表和深部分布礦脈有近百條,+934 m中段以上主要有L27、L7、L27-8、L27-14、L27-15等35種礦脈,礦體分散但品位高。礦體走向近東西向,平均厚度2.6 m,局部厚度8~15 m,礦體北傾,傾角50°~70°,屬熱液充填交代型礦床。礦體上下盤圍巖主要為蝕變玄武巖,少量蝕變凝灰?guī)r及石英脈巖。礦山采用淺孔留礦法采礦,電耙出礦結構,階段高度40 m,沿礦體走向布置采場,采場長30~50 m,頂柱高為5~6 m,底柱高5 m,人行天井寬1.2 m,人行天井兩邊留2~2.5 m礦柱。
礦區(qū)地表有少量的工業(yè)廠房、柏油路等建(構)筑物。井下采空區(qū)主要分布在934 m、974 m、1 014 m及1 054 m中段。
以各中段采空區(qū)影響范圍為研究區(qū)域,具體建模尺寸需要考慮FLAC3D模型大小對數(shù)值模擬結果的影響,并兼顧礦山開采實際及計算機計算速度,根據(jù)彈塑性力學理論-圣維南原理合理確定。據(jù)此對礦區(qū)地表、礦體及各采空區(qū)數(shù)據(jù)進行提取和分析。最終建立的模型范圍X方向:3 850~5 550;Y方向:5 100~5 800;Z方向:+750 m至地表。
由于礦體形態(tài)復雜,采空區(qū)數(shù)量眾多,F(xiàn)LAC3D軟件在復雜模型的構建上存在一定缺陷,因此采用AutoCAD、Dimine及Midas GTS NX等軟件聯(lián)合構建地質(zhì)模型。Dimine軟件僅用于三維實體模型生成,實現(xiàn)模型可視化。為了將三維模型應用到數(shù)值模擬仿真計算,選用Midas GTS NX軟件作為媒介進行模型傳遞。Dimine模型經(jīng)AutoCAD處理后導入Midas GTS NX中后進行面縫合、實體擴展、布爾運算、網(wǎng)格劃分等操作。地質(zhì)模型見圖1。
數(shù)值分析必須將模型計算區(qū)域進行離散化處理,通過轉(zhuǎn)換軟件將模型的單元、節(jié)點和組信息導入FLAC3D軟件中。為了合理規(guī)劃模型單元數(shù)量、提高計算效率、適應FLAC3D軟件計算要求,對各中段采空區(qū)及礦體進行加密處理,其中,各采空區(qū)單元尺寸3 m,圍巖36 m、次級單元18 m,最終生成的數(shù)值分析模型如圖2所示,模型單元總數(shù)約131.5萬個,節(jié)點總數(shù)約22.5萬個。
采用數(shù)值模擬分析礦山開采對地表及建構筑物的影響程度和范圍。數(shù)值模擬計算的前提是確定材料本構模型、力學參數(shù)、邊界條件等。
1) 本構模型。選擇摩爾-庫侖本構模型,另外,開挖部分采用空單元(null)模型。
2) 材料力學參數(shù)。確定模型的材料屬性是FLAC3D軟件計算的關鍵。根據(jù)礦圍巖室內(nèi)試驗,考慮巖體節(jié)理裂隙等客觀因素,并通過RocLab軟件折減計算,得到各巖體力學參數(shù),見表1。
3) 邊界條件。數(shù)值分析中初始應力場的施加十分關鍵,施加有偏差的應力場會導致數(shù)值計算的失敗。本次按巖體自重應力施加荷載,在模型底部施加Z向約束,四周分別施加X向約束,Y向約束,地表為自由面。
4) 計算方案。本次數(shù)值模擬根據(jù)礦區(qū)934~1 054 m中段實際開采區(qū)域進行模擬計算。整個模擬計算分五步進行,首先初始平衡,之后對1 054 m中段、1 014 m中段、974 m中段、934 m中段分四個步驟進行開挖計算。
圖1 地質(zhì)模型Fig.1 Geological model
圖2 數(shù)值分析模型Fig.2 Numerical analysis model
表1 模型力學參數(shù)Table 1 Model mechanics parameters
巖性容重/(kN/m3)體積模量/GPa剪切模量/GPa抗拉強度/MPa泊松比內(nèi)聚力/MPa內(nèi)摩擦角/(°)礦體28.504.504.100.180.152.6148.35蝕變玄武巖28.106.096.670.220.274.2043.74
按照地表移動變形值大小及其對建(構)筑物及地表的影響程度,可將地表移動盆地劃分為最外邊界、移動邊界和裂縫邊界三個邊界(圖3)。對應描述地表移動的參數(shù)主要有邊界角、移動角及裂縫角。
1) 最外邊界:一般指下沉1 mm的點圈定的邊界,如圖3中的ACBD。
2) 移動邊界:是指臨界變形值確定的邊界,根據(jù)《建筑物、水體、鐵路及主要井巷煤柱留設與壓煤開采規(guī)程》(簡稱“三下規(guī)程”)中規(guī)定:使用破壞等級與地表變形的對應關系最低等級的三個參數(shù)臨界值來劃定地表影響范圍,即傾斜i=3.0 mm/m、曲率K=0.2 mm/m2、水平變形ε=2.0 mm/m。最終劃定的移動邊界如圖3中的A’C’B’D’。
3) 裂縫邊界:是指移動盆地的最外側的裂縫圈定的邊界,如圖3中的A”C”B”D”。
根據(jù)“三下規(guī)程”,地下開采地表建筑物的損壞等級由地表水平變形、曲率、傾斜三個指標判定。地表變形的大小和建筑物本身抵抗變形的能力決定了地表建(構)筑物的損壞程度。地表傾斜反映地表移動盆地沿某一方向的坡度;曲率反映彎曲程度的變化;水平變形表示單位長度線段的拉伸或壓縮。
4.2.1 地表位移分析
礦山地表整體位移等值線圖如圖4所示。為描述地表整體位移變化趨勢,選用0.01 m為趨勢邊界。隨著934~1 054 m中段之間多個采場回采后,巖層位移量逐漸增大,最大位移約為12 cm,主要集中在開挖區(qū)域的頂?shù)装逯醒胛恢?,主要為礦體的上盤位置。同時,隨著開挖的進行,位移影響范圍也逐漸增大,最終在地表形成位移約為1 cm的影響區(qū)域。該區(qū)域主要集中在L7脈板房至水塔之間,東西范圍約320 m,南北范圍約265 m。最終,開挖后地表全位移最大值為0.0108 m。
為進一步分析地下開采對地表的影響范圍,沿礦體走向選取剖面進行分析。如Y=5 450剖面,見圖5。由圖5可以看出,地下礦體開挖后的影響范圍將波及地表,主要集中在Y=4 500~4 900之間,位移量約為0.01 m。
圖3 地表移動盆地邊界示意圖Fig.3 The schematic diagram of the boundary of the surface subsidence basin
圖4 地表整體位移等值線圖Fig.4 Contour map of surface whole displacement
圖5 Y=5 450處整體位移等值線云圖Fig.5 Contour map of whole displacement in Y=5 450
4.2.2 地表傾斜分析
地表傾斜見圖6~7。在東西向傾斜方面,地下開采在地表形成了一正一負兩個傾斜區(qū)域,以X=4 600~4 700為界,分別位于東西兩側。正傾斜、負傾斜區(qū)域的傾斜絕對值最大分別為4e-5m/m、3.5e-5m/m,均未達到3 mm/m,小于最低等級要求。在南北向傾斜方面,地表形成了一正一負兩個傾斜區(qū)域,以Y=5 450~5 550為界,分別位于南北兩側。正傾斜、負傾斜區(qū)域的傾斜絕對值最大分別為4e-5m/m、4.5e-5m/m,均未達到3 mm/m。
圖6 地表東西向傾斜Fig.6 Surface east-west tilt
圖7 地表南北向傾斜Fig.7 Surface north-south tilt
4.2.3 地表曲率分析
地表曲率見圖8~9。在東西向曲率方面,地表形成了一正一負兩個曲率區(qū)域,正曲率、負曲率區(qū)域的曲率絕對值最大分別為4e-6/m、3.5e-6/m,均未達到2e-4/m。在南北向曲率方面,地表形成了一正一負兩個曲率區(qū)域,正曲率、負曲率區(qū)域的曲率絕對值最大分別為1.2e-5/m、1.6e-5/m,均未達到2e-4/m。
圖8 地表東西向曲率Fig.8 Surface east-west curvature
圖9 地表南北向曲率Fig.9 Surface north-south curvature
4.2.4 地表水平變形分析
地表水平變形見圖10~11。在東西向水平變形方面,地表形成了一正一負兩個水平變形區(qū)域,以正變形區(qū)域居中,負變形區(qū)域位于其兩側,正變形、負變形區(qū)域的變形絕對值最大分別為3e-3mm/m、2.1e-2mm/m,均未達到2 mm/m。在南北向水平變形方面,地表形成了一正一負兩個水平變形區(qū)域,正變形、負變形區(qū)域的變形絕對值最大為2.2e-2mm/m、3.6e-2mm/m,均未達到2 mm/m。
為了界定地下開采對地表的影響范圍,根據(jù)國家的相關標準及規(guī)范對巖層移動形成的地表位移、傾斜、曲率及水平變形進行分析,并通過地表變形值的比較,最終圈定地表移動影響范圍。
礦區(qū)934~1 054 m中段礦體開挖以后,采空區(qū)主要影響區(qū)域是礦體上盤,同時由于礦體沿東西走向較長,礦體在X=4 500~4 900之間比較集中,且有L27-8厚礦體(礦脈平均厚度達到4.68 m,最厚處可達16.3 m)存在,因此,在地表的開挖影響區(qū)域主要集中在X=4 500~4 900之間,這與地表整體位移及豎向位移分析結果一致,同時也與急傾斜礦體地下開采圍巖的變形破壞規(guī)律[7,9]相吻合。因此,根據(jù)地表整體位移及地表豎向位移的分析,按照地表移動邊界中最外邊界(一般指下沉10 mm的點圈定的邊界)的圈定原則,以下沉10 mm作為臨界值圈定地表影響范圍(圖12),可見圈定的地表影響范圍具有一定的可信度。
圖10 地表東西向水平變形Fig.10 Surface east-west horizontal deformation
圖11 地表南北向水平變形Fig.11 Surface north-south horizontal deformation
圖12 地表影響范圍Fig.12 The affected area of the surface
1) 借助AutoCAD、Dimine、Midas GTS NX等軟件聯(lián)合構建三維數(shù)值計算仿真模型,并通過轉(zhuǎn)換軟件生成FLAC3D文件,這一方法提高了建模效率,能夠滿足復雜地質(zhì)模型的構建及三維計算分析的要求。
2) 地表傾斜、曲率及水平變形均小于磚混結構建筑物破壞等級要求的臨界變形值,表明地表傾斜、曲率及水平變形均不會對地表建(構)筑物造成影響。
3) 地下采空區(qū)主要影響區(qū)域是礦體上盤,最終開挖影響區(qū)域主要集中在礦體分布集中區(qū)域?qū)纳媳P位置,這與地表位移分析結果一致,同時也與急傾斜礦體地下開采圍巖的變形破壞規(guī)律相吻合。
4) 以地表移動邊界中最外邊界(下沉10 mm為臨界值)圈定地表影響范圍具有一定的可信度。