基于FLAC3D模型動態更新的自然崩落法開采過程模擬技術*

方傳峰，余傳玉，史存丁，馮興隆

(1.中南大學 資源與安全工程學院，湖南 長沙 410083；2.云南迪慶有色金屬有限責任公司，云南 香格里拉，674400)

0 引言

自然崩落法作為一種低成本，高效率，以及安全性好的采礦方法，適用于大型低品位礦山的開采，自1895年在Pewabic鐵礦首次應用起，技術和裝備不斷發展，在國內外得到了大量應用[1-2]。該方法是通過簡單的拉底和切割工程，在拉底空間上方利用礦體自身的軟弱結構面，在重力及次生構造應力作用下，誘導礦體失穩崩落，從而省卻大量鑿巖爆破工程及費用。

自然崩落法開采過程主要受礦床開采技術條件，以及拉底和出礦過程控制，通常采用經驗推導法、理論計算法、相似物理試驗法以及數值模擬法等分析和確定上述因素對崩落礦巖時空演化特征的影響[3-5]。隨著計算機技術的發展，數值模擬法已成為國內外研究學者用于預測自然崩落力學響應特征的有效途徑。王濤等[6]采用PFC2D研究礦巖在崩落過程中的裂紋擴展及接觸力分布情況；王連慶等[7]以某鎳銅礦的地質條件及礦巖物理力學性質為背景,采用PFC2D數值模擬的方法分析了自然崩落法的崩落規律；孫闖等[8]采用FLAC3D軟件，在獲取泥巖的粘結力軟化和內摩擦角硬化參數的基礎上預測大跨度的泥質頂板冒落失穩； Vyazmensky等[9]采用實例統計、經驗圖分析與有限—離散元軟件數值模擬研究地表沉降的影響因素；Woo等[10]采用FLAC3D與干涉合成孔徑雷達(InSAR)分析、評價與預測礦塊崩落法引起的地表變形，并修正數值模擬參數。

PFC離散元方法不能進行大規模開采模型的計算和分析[11]，FLAC3D雖可以分析大規模問題，但在以往的應用研究中，未引入網格動態生成機制，不能真實反映礦巖崩落和出礦過程的動態性，也不能及時生成底部結構承重礦堆，從而影響了對自然崩落過程力學響應研究的準確性。本文基于FLAC3D軟件設計工作流程，并利用軟件內嵌FISH邏輯語言進行編程與模擬，及時消除已崩落單元且同時在底部結構生成相應質量礦堆，以實現開挖過程的動態模擬。以西南某銅礦為例，通過建立二維礦山模型，對其進行拉底模擬，以分析顯式崩落流程設計的有效性。

1 崩落過程動態模擬技術流程

以實際礦山崩落為動態破壞過程。傳統連續介質算法模擬自然崩落，雖可識別各步拉底對應崩落區，但未做破壞礦巖單元的動態更新，既不能體現崩落過程的動態性，更重要的是導致模擬結果偏于保守，與實際工程相比，無論在空間范圍，還是在崩落量上，均出現較大偏差，制約了生產計劃和監測方案的準確制定。在對自然崩落法開采過程特征進行分析的基礎上，提出和設計了基于FLAC3D的過程動態模擬算法，編制了相應的代碼。

1.1 模型動態更新實現思想

首先，按照傳統FLAC3D力學分析步驟和方法建立模擬模型。

其次，在設定每一步驟拉底范圍后，采用摩爾-庫倫模型模擬計算其力學性態，并按照考慮位移、最大主應力或塑性狀態的單指標判斷準則，或同時考慮多個因素的綜合指標判斷準則，確定潛在破壞和需要動態處理的單元。

最后，利用FLAC3D編寫FISH函數修改單元的本構模型，實現破壞單元的消除和崩落面和底部結構之間單元的顯現，即將發生破壞的單元類型修改為空模型(null)，原屬于空模型的單元修改為其他本構類型，并修改其力學參數以表征松散礦堆。

1.2 模型動態更新基本原則

基于FLAC3D單元的消除與再現[12]，程序主要遵從3個基本原則，實現模型動態更新和崩落過程模擬。

1)質量守恒原則

礦巖破壞后發生冒落，若忽略礦巖的吸水性且考慮出礦作用，則質量恒定，滿足式(1)：

m崩=m出+m堆

(1)

式中：m崩為崩落礦石質量；m出為出礦量；m堆為出礦后礦堆質量。

2)體積平衡原則

由于破碎后礦巖自身的碎脹性，體積會膨脹。為滿足底部結構存在足夠體積的空區順利生成礦堆，則體積應滿足式(2)：

V空+V崩≥V堆

(2)

式中：V空為已存在空區體積；V崩為待崩落礦巖體積；V堆為待沉積礦堆體積。

3)拱形發育原則

底部結構沉積的礦堆呈散體狀態，礦堆按照一定的自然安息角呈拱形發育。程序根據現有空區平面形態和位置，按設定的安息角在一定拱形范圍內重現單元，其形態見圖1。

1.崩落區；2.已開挖拉底層；3.已開挖聚礦槽；4.已開挖進路；5.礦堆。圖1 礦堆示意Fig.1 Mine heap schematic

1.3 動態過程模擬流程

基于上述3項原則，引入人工控制變量，編制Fish程序，具體動態過程模擬流程見圖2。

圖2 顯式崩落實現概圖Fig.2 The overview of realization of explicit caving

圖2為顯式崩落實現概圖，首先確定每步拉底范圍，其次對第N步模擬拉底，運行平衡后，根據破壞判斷準則選定破壞單元，基于三項原則實現破壞單元的消除，并在底部結構生成相應質量礦堆，實現該拉底步顯式崩落操作后，再次平衡模型。隨后進入第N+1步拉底模擬，以此周而復始，最終實現各拉底步的動態崩落。

2 應用研究

2.1 模型構建及地應力反演

依據西南某銅礦地形等高線構建X方向長1 000 m，厚度2 m，地形復雜的二維數值模型，并根據井巷工程布置構建拉底層、聚礦槽、出礦進路等地下工程，最終模型見圖3。

由于研究內容為地下開挖問題，模擬選用的本構模型為摩爾-庫倫模型。顯式崩落模擬需巖體參數、礦堆參數、接觸面參數，其中巖體參數采用胡克-布朗(Hoek-Brown)強度折減準則并利用RocLab軟件，通過輸入完整巖體參數，選擇符合礦山實際的巖體條件GSI(地質強度指標)、mi(完整巖石參數)、D(擾動因數)最終確定；礦堆為散體，故強度參數設置為0，彈性參數參照文獻[13-14]所提供曲線，根據本次模擬特征近似確定；接觸面按照常規取值[15]，具體參數見表1～表3。

圖3 礦山模擬模型Fig.3 The model of mine

表1 礦山巖體參數Table 1 The mechanical parameters of mine rock mass

表2 礦堆巖體參數Table 2 The mechanical parameters of ore heap rock mass

表3 接觸面參數Table 3 The mechanical parameters of interface

該銅礦地應力的測量主要集中在進路水平，將模型X與Y方向滾支固定，底端完全固定，地表作為自由面，施加初始應力，進行地應力反演。最終進路水平地應力反演值與實際監測值對比見表4。

表4 初始地應力反演結果Table 4 The result of back analysis of origin geostatic stress

模擬過程需人為設定若干參數，且模擬采用的出礦模式為相對出礦量，所需參數見表5。

表5 人工設定參數Table 5 Manual setting parameter

2.2 崩落區及礦堆演化特征

平衡完畢后的模型按圖3模擬開挖：開挖始自拉底層中段，跨度15 m，隨后每步開挖均向拉底層兩側各推進15 m，聚礦槽與出礦進路滯后拉底推進線30 m開挖，各步驟開挖后運行前均按圖2所示流程處理前一步模擬結果，實現顯式崩落。共模擬15步，各步模擬結果匯總見圖4～圖6。

圖4 崩落區演化特征Fig.4 Evolution characteristics of caving area

圖4中，隨拉底向兩側推進，崩落區整體呈穹隆狀向上、向四周延伸，該現象與工程實際經驗相符；拉底第13步時崩透地表，隨后的拉底過程中地表出現大面積破壞，基于地表散體間的摩擦效果，將與水平呈45°范圍內的破壞區作為最終崩落區。

圖5 礦堆形態演化特征Fig.5 Evolution characteristics of ore heap

各開挖步驟礦堆堆積模擬效果見圖5，礦堆于崩透地表前所呈形態為頂端帶有一定自然安息角的弧形，一定高度以下崩落區完全填充，同時崩落下的礦堆未侵占未崩落單元。當第13步崩透地表后，引起地表大范圍破壞，圖4中第13步開挖對應的崩落區將按照與地表輪廓平行的狀態進行填充，最終實現地表下沉效果。因此，通過圖2所示的顯式崩落流程可近似實現礦體的崩落下降效果和礦堆在底部結構的堆積效果。

圖6 體積對比統計Fig.6 The plot of volume contrast statistics

圖6為各項體積統計圖，理論礦堆體積按初始出礦率換算所得，底部空隙體積為累計可填充礦堆體積，實際崩落體積為最終在底部結構實際生成的礦堆體積。工程模擬過程中，由于碎脹系數與出礦量的相對關系，理論礦堆體積始終大于崩落體積；當底部空隙體積大于理論礦堆體積時，崩落下的散體礦石可按計劃出礦，此時理論礦堆體積與實際礦堆體積相當；當模擬至15步時，底部空隙體積已小于理論礦堆體積，此時認為礦堆已將底部填滿，故將第15步的崩落礦量在原有出礦率基礎上多出礦10%，所以實際礦堆體積在15步中較理論礦堆小。由圖6可知所設計流程能夠自動實現對礦堆是否完全填充空隙的判斷，進而動態調整出礦率。

3 2種方式崩落高度與體積對比分析

為更好理解崩落區不做沉降處理的傳統崩落模擬方式與顯式崩落模擬對崩落過程的不同影響，選用同一模型通過傳統崩落模擬方式進行拉底模擬，對比2種方式在崩落高度與崩落體積2方面的差別。

以拉底層上端為基準水平，將每一步崩落區最高值作為崩落高度，圖7為2種模擬方式的崩落高度與實際監測頂板高度對比。

由于工程與地質原因，頂板監測數據有若干缺失，但并不影響2種模擬方式對比。顯式崩落模擬相較傳統模擬，與監測數據吻合度更高，更能反映頂板崩落高度的實際變化趨勢。隨拉底推進，顯式崩落速度明顯快于傳統崩落，其崩透地表時預測拉底步僅比實際崩透地表晚一步。傳統模擬方式頂板崩落高度變化明顯滯后于顯式崩落和實際監測值。原因是顯式崩落移除了每步滿足破壞準則的單元，進而促進了崩落區進一步向上發展，所以與傳統崩落模擬方式相比，崩落速度更快且更符合實際情況。

圖7 崩落高度曲線Fig.7 The plot of caving height

由圖8所示，與傳統崩落模擬方式相比，顯式崩落模擬方法誘發的崩落量更大，原因是在每一步運行時，及時移除了上一步的崩落區單元，促進了與上一步崩落區緊貼的圍巖發生更大范圍的破壞。因此，從顯式崩落的高度與體積均大于傳統方式可知：顯式崩落模擬能夠體現崩落區消除后所產生的空區對下一步崩落的促進效果，該效果有助于提高對崩落預測的準確性。

圖8 2種模擬方式崩落體積Fig.8 The plot of caving volume of two simulation methods

4 結論

1)設計并優化了顯式崩落算法流程并利用FLAC3D內嵌的編程語言，實現了自然崩落法中崩落體與圍巖的脫離，并完成了礦堆在礦柱上方的負載施加，有利于分析底部結構的穩定性。

2)實現了崩落區形態的自動識別，并引入礦堆碎脹系數、出礦量、自然安息角等變量，便于對礦山工程實際問題進行特定分析。

3)設計流程能自動監測礦堆體積與底部空隙體積的相對關系，進而實現出礦量的動態調整。

4)傳統與顯式2種模擬結果與監測數據對比，突出了崩落區消除后的空區對下一步崩落的促進效果和顯式崩落模擬的準確性。

5)雖然該優化算法用于解決自然崩落法模擬過程中存在的問題，但其思路也適用于煤礦開挖等與開挖崩落相關的工程過程的模擬。

[1] WOO K S, EBERHARDT E, ELMOl D, et al. Empirical investigation and characterization of surface subsidence related to block cave mining[J]. International Journal of Rock Mechanics & Mining Sciences, 2013, 61(61):31-42.

[2] ELMO D, STEAD D. An Integrated Numerical Modelling-Discrete Fracture Network Approach Applied to the Characterisation of Rock Mass Strength of Naturally Fractured Pillars[J]. Rock Mechanics & Rock Engineering, 2010, 43(1):3-19.

[3] 趙文. 巖石力學[M]. 長沙：中南大學出版社, 2014:187-188.

[4] BRADVB H G, BROWNE T. 地下采礦巖石力學[M]. 佘詩剛, 朱萬成, 趙文，等譯. 北京: 科學出版社, 2011.

[5] 薛東杰, 周宏偉, 任偉光, 等. 淺埋煤層超大采高開采柱式崩塌模型及失穩[J]. 煤炭學報, 2015, 40(4):760-765.

XUE Dongjie, ZHOU Hongwei, REN Weiguang, et al. Instability of pillar collapse model and generation of cracks in thick coal seam mining at shallow depth[J]. Journal of China Coal Society, 2015,40(4):760-765.

[6] 王濤, 盛謙, 熊將. 基于顆粒流方法自然崩落法數值模擬研究[J]. 巖石力學與工程學報, 2007, 26(S2):4202-4207.

WANG Tao，SHENG Qian，XIONG Jiang．Research on numerical simulation based on particle flow method[J]．Chinese Journal of Rock Mechanics and Engineering，2007,26(S2):4202-4207.

[7] 王連慶, 高謙, 王建國, 等. 自然崩落采礦法的顆粒流數值模擬[J]. 北京科技大學學報, 2007, 29(6):557-561.

WANG Lianqing, GAO Qian, WANG Jiaguo, et al. Particle flow code numerical simulation of natural caving mining method[J]. University of Science & Technology Beijing, 2007, 29(6):557-561.

[8] 孫闖, 張向東, 張濤,等. 深部大跨度泥質頂板剪切冒落失穩區預測研究[J]. 中國安全生產科學技術, 2016, 12(1):23-27.

SUN Chuang, ZHANG Xiangdong, ZHANG Tao, et al. Prediction on shear caving instability region of deep and large span muddy roof[J]. Journal of Safety Science and Technology, 2016,12(1):23-27.

[9] VYAZMENSKY A, ELMO D, STEA D. Role of Rock Mass Fabric and Faulting in the Development of Block Caving Induced Surface Subsidence[J]. Rock Mechanics & Rock Engineering, 2010, 43(5):533-556.

[10] WOO K S, EBERHARDT E, RABUS B, et al. Integration of field characterisation, mine production and InSAR monitoring data to constrain and calibrate 3-D numerical modelling of block caving-induced subsidence[J]. International Journal of Rock Mechanics & Mining Sciences, 2012, 53(9):166-178.

[11] 曾慶田,劉科偉,嚴體,等.基于多數值模擬方法聯合的自然崩落法開采研究[J].黃金科學技術,2015,23(1):66-73.

ZENG Qingtian, LIU Kewei, YAN Ti, et al. Study on natural caving mining method based on multi-numerical simulation method[J]. Gold Science and Technology, 2015, 23(1): 66-73.

[12] 王濤, 韓煊, 趙先宇, 等. FLAC3D數值模擬方法及工程應用[M]. 中國建筑工業出版社, 2015:170-171.

[13] 陸永青, 計三有, 戴詩亮,等. 散體材料彈性參數的估計[J]. 武漢理工大學學報(交通科學與工程版), 1997(2):178-185.

LU Yongqing, JI Sanyou, DAI Shiliang, et al. Estimate of elastic parameters for granular material[J]. Journal of Wuhan Transportation University, 1997, 21(2): 178-185.

[14] 賈林剛,劉卓然.充填開采充填率與地表移動規律的數值模擬研究[J].遼寧工程技術大學學報(自然科學版),2015,34(9):1010-1015.

JIA Lingang，LIU Zhuoran．Study on numerical simulation of filling rate and surface subsidence in backfill mining[J]．Journal of Liaoning Technical University (Traffic Science and Engineering Edition)，2015，34(9)：1010-1015.

[15] 曹帥,杜翠鳳,母昌平,等.崩落法轉充填法采礦地表移動二維離散元程序數值模擬及其規律驗證[J].巖土力學,2015,36(6):1737-1751.

CAO Shuai，DU Cuifeng，MU Changping，et al．UDEC based modelling of mining surface movement due to transforming from block caving to sublevel filling and its law verification[J]．Rock and Soil Mechanics, 2015, 35(6) :1737-1751.