基于FLAC3D 對開采引起地表沉降規律研究

王志凱 ,宋文龍

(1.五礦有色金屬股份有限公司,北京 100043;2.鑫聯環?？萍脊煞萦邢薰?北京 100027)

0 引言

地下開采造成了巖土體內部原有應力平衡破壞,往往導致井下采場圍巖冒落、變形、移動甚至斷裂[1-2]。圍巖的破壞和移動有時會波及地表,致使地表發生不同程度的變形和破壞,進而嚴重影響井下人員的作業安全[3]。因此,如何利用一定的理論和數值方法對礦山開采引起地表沉降規律進行深入分析研究是保障礦區安全和人員作業安全具有重要意義。

目前,在我國主要采用經驗公式法、工程類比法、現場監測法等方法對淺部礦產開采引起地表沉降影響范圍和影響程度的研究[4]。但是隨著開采深度的不斷加深,深部開采涉及“三高一擾動”問題,同時井下環境惡劣,上述方法常帶來各種弊端。隨著新技術、新理論、新方法的發展,尤其計算技術的應用,對井下開采造成地表移動與變形影響研究取得了許多進展。劉寶琛等[5]采用隨機介質理論分析了開采中巖石移動問題。李增琪等[6]通過層間遞推公式將多層體礦壓和巖層移動的計算得到了巖層移動的三維層體模型。唐春安等[7]研究開發RFPA2D能實現開采過程中覆巖變形、離層、斷裂及地表沉陷規律模擬。

為此,本文基于FLAC3D運用有限差分法從地表垂直位移、水平位移、斜率與曲率三個維度對某金礦開采引起地表沉降問題進行數值分析,探究其引起地表沉降的規律,以期為該礦安全開采提供一定保障。

1 工程背景

某金礦已開采至-730 m 中段,開拓深度已至-1 030 m 水平,深部礦體傾斜深至-1 500 m 以下。Ⅰ、Ⅴ號礦體是該金礦主礦體,兩者占礦山資源總量的96.42%。Ⅰ號礦體上部走向長約240 m,下部300 m,平均厚度30 m,傾角平均30°。先后開拓了-10～-580 m 的17 個中段,中段高度20～50 m,-10 m 以上標高作為永久保安礦柱。

Ⅴ號礦體為中間厚上下兩邊薄的透鏡狀礦體,平均厚度80 m,上部走向長約240 m,下部走向300 m,傾角30°。先后開拓了-480～-730 m 六個階段,階段高度50 m。

根據該礦的實際開采情況,確定研究對象為Ⅰ號和Ⅴ號整個礦體,模型簡化Ⅰ號礦體賦存高度-10～-580 m,走向長300 m,平均厚度30 m,傾角30°。Ⅴ號礦體賦存高度-480～-730 m,走向長240 m,厚度80 m,傾角30°。

2 礦體的模型建立與參數的選取

2.1 礦體模型的建立

根據礦體賦存條件計算方案中計算區域取1 800 m×1 500 m×800 m。模型采用巖土工程數值模擬軟件建立模型及劃分網格,共劃分802 872 單元,831 480 個節點,計算模型示意圖如圖1 所示。

圖1 計算模型示意圖

2.2 模擬參數的選取

2.2.1 參數與本構模型的選取

根據工程實際情況,選定模擬計算采用的巖體力學參數如表1 所示,數值計算中主要采用摩爾-庫倫本構模型。

表1 巖體力學參數

2.2.2 邊界條件

邊界應力與約束條件是計算模型的重要內容,直接影響計算結果的可靠性及精度。應用應力解除法測量某金礦地應力得出的最大水平主應力、最小水平主應力、垂直主應力關系式:

式中σh.max—最大水平主應力;

σh.min—最小水平主應力;

σv—垂直主應力,MPa;

H—測點埋深,m。

采用基于應力邊界的快速地應力場擬合方法,應用此法形成的初始應力場如圖2 所示,與實測地應力場有較高的吻合度。

圖2 初始應力場云圖

3 開采過程地表沉降規律分析

3.1 地表豎直位移特征及分析

(1)地表垂直位移特征分析

圖3、4 分別為I #礦體和V #礦體-730 階段開采結束地表沉降位等直線圖。從圖中可以看出,I #礦體開采結束后,地表沉陷位置比較集中,地表最大垂直位移達到了148 mm,影響范圍沿走向方向約900 m,垂直于走向方向約1 100 m。最大下沉量出現在礦山的上盤位置,約為礦體走向方向的中心位置。隨著開采向深部進行,沉降中心逐漸向V #礦體上部移動,V #礦體-730 階段開采結束后地表下沉幅度和影響范圍進一步增大;地表沉陷中心略向V#礦體上方偏移,沉陷盆地形狀基本沒有什么變化,最大下沉量達到166.5 mm,增大幅度達到18.5 mm,影響范圍沿走向方向約900 m、垂直走向方向約1 300 m。表明地表沉降主要由于淺部資源開采引起,深部開采對地表沉降的影響不大。通過等高線圖還可以很清楚地看到地表沉降形成凹陷區域的平面圖,隨V #礦體開采深度的增加,地表沉降逐漸擴大,并且沉降中心逐漸向V #礦體上方移動。等值線在垂直于走向方向由淺部開采到深部開采過程中逐漸變疏,A 較B 位置較密,可知在A 位置地表傾斜變化率較大,此位置地標建筑物更容易產生破壞。

圖3 I#礦體開采結束后地表沉降等值線圖

圖4 V#礦體-730 階段開采結束地表沉降等直線圖

3.2 地表水平位移特征及分析

通過圖5、圖6 可以看出,隨著各礦體相續開采,地表水平移動范圍不斷擴大,該金礦最大水平主應力垂直于走向方向,形成了以采空區為中心,并且地表水平位移向空區兩側不斷減小的水平移動帶,水平移動范圍基本呈現對稱分布。I #礦體開采結束水平移動最大值約為170 mm,V #礦體-730 階段開采結束后正向移動最大值約為180 mm,增加10 mm,幅度不大。A～C 區間水平移動值逐漸減小,以拉伸變形為主,A 區域水平移動值最大,A～B 區域逐漸減小,表明在拉伸變形到達一定平衡時轉而以壓縮變形為主。B 區域等值線較密表明此處地表建筑受壓縮變形破壞影響較大。

圖5 I#礦體開采結束后地表x 水平位移等值線圖

圖6 V#礦體-730 階段開采結束地表x 水平位移等值線圖

3.3 地表斜率與曲率的數值微分計算

通過在地表坐標設置歷史記錄點,將數值計算后地表的垂直位移和水平位移自動記錄并保存,才能進一步計算地表變形曲率K、傾斜i、水平變形值ε。記錄點的間距跟礦體的平開采深度呈正相關。開采深度越大,記錄點的間距也越大。在開采深度大于500 m 時,選取的點間距應為50 m,所以在地表每50 m 設置一個記錄點,記錄水平和垂直位移值。通過仿真軟件計算結束后生成的x、z兩個方向位移值,得出的地表所有點的結果如表2 所示。

表2 地表監測點的水平和豎直位移值

數值微分問題指已知函數y=f(x)的一組離散值(xk,f(xk))(k=1,2,…,n),求f(x)在節點處的導數或微分值。對離散點求導數,使用數值微分方法比較方便。

根據向前差商或向后差商,分別得到微分方程:

對式(4)(5)取算術平均值,并用二階中心差商近似f″(x0),便得到f″(x0)的數值微分公式(6)。

通過求式(6)誤差漸進展開式。根據Taylor 公式可得:

令G(h)為上面兩式等號左端的第二項,對于一階和二階導數可建立外推公式:

計算時,步長分別取100 m 和50 m,即可按照所得的Gm+1(h)值計算出相應的一階導數和二階導數,即地表變形的傾斜和曲率。

計算后所得地表最大傾斜0.6 mm/m,最大曲率0.005 2 mm/m2,最大水平變形0.404 mm/m。對比《建筑物、水體、鐵路及主要井巷煤柱留設與壓煤開采規程》中建筑物保護等級及允許變形值分類。整個地表出現的變形值都沒有超過一般磚木結構允許的臨界變形值。因此,充填開采后地表的變形不會對地表建(構)筑物的安全造成較大的危害。

4 結論

(1)通過FLAC3D對地表豎直位移分析得到I#礦體開采結束后,地表沉陷位置比較集中,沉陷盆地形狀較規整,地表最大下沉量達到了144.9 mm,影響范圍沿走向方向約900 m,垂直于走向方向約1 100 m。V#礦體-730 階段開采結束后,最大下沉量達到166.5 mm,增大幅度達到21.6 mm,影響范圍沿走向方向約900 m、垂直走向方向約1 300 m。

(2)從地表水平位移得到,I #礦體開采結束正向移動最大值約162 mm,V #礦體-730 階段開采結束后正向移動最大值約175 mm。

(3)通過地表斜率與曲率的數值微分計算得到,地表最大傾斜0.6 mm/m,最大曲率0.005 2 mm/m2,最大水平變形0.404 mm/m,均沒有超過一般磚木結構允許的臨界變形值。因此,開采后對采空區進行及時充填,地表的變形不會對地表建(構)筑物的安全造成較大的危害。