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

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

0 引言

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

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

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

1 崩落過程動態(tài)模擬技術(shù)流程

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

1.1 模型動態(tài)更新實現(xiàn)思想

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

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

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

1.2 模型動態(tài)更新基本原則

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

1)質(zhì)量守恒原則

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

m崩=m出+m堆

(1)

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

2)體積平衡原則

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

V空+V崩≥V堆

(2)

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

3)拱形發(fā)育原則

底部結(jié)構(gòu)沉積的礦堆呈散體狀態(tài)，礦堆按照一定的自然安息角呈拱形發(fā)育。程序根據(jù)現(xiàn)有空區(qū)平面形態(tài)和位置，按設(shè)定的安息角在一定拱形范圍內(nèi)重現(xiàn)單元，其形態(tài)見圖1。

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

1.3 動態(tài)過程模擬流程

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

圖2 顯式崩落實現(xiàn)概圖Fig.2 The overview of realization of explicit caving

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

2 應(yīng)用研究

2.1 模型構(gòu)建及地應(yīng)力反演

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

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

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

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

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

表3 接觸面參數(shù)Table 3 The mechanical parameters of interface

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

表4 初始地應(yīng)力反演結(jié)果Table 4 The result of back analysis of origin geostatic stress

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

表5 人工設(shè)定參數(shù)Table 5 Manual setting parameter

2.2 崩落區(qū)及礦堆演化特征

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

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

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

圖5 礦堆形態(tài)演化特征Fig.5 Evolution characteristics of ore heap

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

圖6 體積對比統(tǒng)計Fig.6 The plot of volume contrast statistics

圖6為各項體積統(tǒng)計圖，理論礦堆體積按初始出礦率換算所得，底部空隙體積為累計可填充礦堆體積，實際崩落體積為最終在底部結(jié)構(gòu)實際生成的礦堆體積。工程模擬過程中，由于碎脹系數(shù)與出礦量的相對關(guān)系，理論礦堆體積始終大于崩落體積；當?shù)撞靠障扼w積大于理論礦堆體積時，崩落下的散體礦石可按計劃出礦，此時理論礦堆體積與實際礦堆體積相當；當模擬至15步時，底部空隙體積已小于理論礦堆體積，此時認為礦堆已將底部填滿，故將第15步的崩落礦量在原有出礦率基礎(chǔ)上多出礦10%，所以實際礦堆體積在15步中較理論礦堆小。由圖6可知所設(shè)計流程能夠自動實現(xiàn)對礦堆是否完全填充空隙的判斷，進而動態(tài)調(diào)整出礦率。

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

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

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

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

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

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

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

4 結(jié)論

1)設(shè)計并優(yōu)化了顯式崩落算法流程并利用FLAC3D內(nèi)嵌的編程語言，實現(xiàn)了自然崩落法中崩落體與圍巖的脫離，并完成了礦堆在礦柱上方的負載施加，有利于分析底部結(jié)構(gòu)的穩(wěn)定性。

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

3)設(shè)計流程能自動監(jiān)測礦堆體積與底部空隙體積的相對關(guān)系，進而實現(xiàn)出礦量的動態(tài)調(diào)整。

4)傳統(tǒng)與顯式2種模擬結(jié)果與監(jiān)測數(shù)據(jù)對比，突出了崩落區(qū)消除后的空區(qū)對下一步崩落的促進效果和顯式崩落模擬的準確性。

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

[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] 薛東杰, 周宏偉, 任偉光, 等. 淺埋煤層超大采高開采柱式崩塌模型及失穩(wěn)[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] 王濤, 盛謙, 熊將. 基于顆粒流方法自然崩落法數(shù)值模擬研究[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] 王連慶, 高謙, 王建國, 等. 自然崩落采礦法的顆粒流數(shù)值模擬[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] 孫闖, 張向東, 張濤,等. 深部大跨度泥質(zhì)頂板剪切冒落失穩(wěn)區(qū)預測研究[J]. 中國安全生產(chǎn)科學技術(shù), 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] 曾慶田,劉科偉,嚴體,等.基于多數(shù)值模擬方法聯(lián)合的自然崩落法開采研究[J].黃金科學技術(shù),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數(shù)值模擬方法及工程應(yīng)用[M]. 中國建筑工業(yè)出版社, 2015:170-171.

[13] 陸永青, 計三有, 戴詩亮,等. 散體材料彈性參數(shù)的估計[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] 賈林剛,劉卓然.充填開采充填率與地表移動規(guī)律的數(shù)值模擬研究[J].遼寧工程技術(shù)大學學報(自然科學版),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] 曹帥,杜翠鳳,母昌平,等.崩落法轉(zhuǎn)充填法采礦地表移動二維離散元程序數(shù)值模擬及其規(guī)律驗證[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.