金屬礦隱患空區動力失穩過程仿真技術

劉曉明 ，羅周全，徐紀成

(1. 中南大學 資源與安全工程學院，湖南 長沙，410083；2. 中國礦業大學 煤炭資源與安全開采國家重點實驗室，江蘇 徐州，221116；3. 中南大學 現代分析測試中心，湖南 長沙，410083)

隱患空區動力失穩過程是一個極其復雜的力學系統，但從局部而言，都是巖石的變形與破裂問題。隨著巖石力學的發展，各種數值計算方法如有限元法、有限差分法、邊界元法、離散元法、拉格朗日元法和流形元法等廣泛應用于各種巖石力學工程問題中[1?4]。而巖層移動與破裂研究的數值計算方法主要包括 2類：一類是建立在連續介質力學基礎上的數值計算方法，如有限元法、邊界單元法等；另一類則是非連續介質力學分析方法，如離散元法。近年來，細觀力學分析方法為巖石破裂過程的數值模擬方法的發展提供了很好的思路。巖石破裂過程分析系統 RFPA(rock failure process analysis)是基于彈性損傷統計模型的數值模擬工具，能較好地模擬巖石破裂失穩過程，目前已廣泛應用于巖石破裂機制以及巖石力學與工程問題的研究，取得了許多成果[5?6]。隱患空區失穩根據失穩類型可以分為靜力失穩和動力失穩。靜力失穩即在原巖應力場作用下進行空區開挖，由于空區周圍應力集中致使空區頂板或兩幫出現拉伸破壞或剪切破壞；動力失穩指在靜力作用的基礎上，受到周圍爆破動荷載的作用，致使空區失穩。與靜力失穩相比，空區動力失穩模式雖然更復雜，仿真難度更大，但是，更能反映實際空區破壞過程。此外，隱患空區形態邊界復雜，如采用矩形、梯形或圓等規則形狀進行簡化，再用簡化的模型進行空區動力失穩仿真，空區的破壞過程與規律往往與實際不符，不能如實反映實際復雜隱患空區動力失穩過程，因此，有必要繼承隱患空區三維信息獲取的結果，在此基礎上進行動力失穩仿真模型構建。華錫集團銅坑礦92號礦體的開采技術條件復雜，井下隱患空區空間分布復雜，而爆破地震波作用下空區發生動力失穩的趨勢日益明顯。為保障92號礦體的安全高效開采，揭示空區動力失穩的內在機制，有效預防空區塌陷事故的發生，在銅坑礦隱患空區 CMS(cavity monitoring system)探測并獲取準確的三維信息的基礎上，借助RFPA提供的巖石破裂過程動力仿真分析功能，同時考慮開挖過程和爆破震動作用對空區穩定性的影響，對空區動力失穩過程進行仿真，以便為開展空區災害預防與控制、空區處理方案的制定、應力監測點布置、資源開采布局等工作提供基礎依據。

1 隱患空區三維信息獲取

1.1 工程概況

柳州華錫集團銅坑礦主要開采對象為細脈帶、91號和92號礦體，這三大礦體在空間位置上重疊產出。目前細脈帶、91號等礦體開采即將結束，正在大量開采92號礦體及細脈帶及91號礦體邊緣薄貧礦體。在長期開采中形成了大量的采空區，局部地段空區的冒落或坍塌已構成威脅礦山生產的重大安全隱患。為能采取有效措施控制因空區引發的地質災害，急需對現有空區進行徹底探查和分析，以掌握空區空間分布、實際邊界、體積等基本情況，并在此基礎上對空區群危險性進行分析，確定不穩定或隨著周圍礦體的開采將不穩定的空區，以便為空區綜合治理與控制、資源有效回收提供技術支撐[7?8]。

1.2 隱患空區三維信息獲取

空區激光探測系統(CMS)因其能準確測定井下危險采空區及其他硐室空間的三維形態而被廣泛應用于國內外金屬地下礦山[9?12]。因此，研究采用CMS對銅坑礦92號礦體隱患空區進行探測，輸出的點云數據經過濾、拼接及坐標、格式轉換，采用三維礦業軟件Surpac生成空區三維實體模型，并以此為基礎開展模型計算、三維可視化及剖面輸出等應用。隱患空區三維信息獲取及可視化過程如圖1所示[13]。

圖1 隱患空區三維信息獲取及可視化過程Fig.1 Process of 3D information acquisition of disaster cavity and its visualization

隱患空區三維空間信息主要包括空區位置形態、體積、頂板面積及邊界等，這些信息將隨時間發生變化，因此，準確、及時掌握隱患空區信息可為礦山進行采礦貧損指標控制、開采設計、空區管理、安全評價及災害預測等提供可靠的基礎性依據。

1.3 計算區域的選取及其特征分析

通過分析銅坑礦 92號礦體探測空區群的三維空間信息，發現92號礦體西南采區內的空區之間互相重疊、貫通，穩定性較差，鄰近頻繁爆破生產作業誘發采空區失穩破壞的趨勢日益明顯。為此，結合目前銅坑礦的開采現狀，選擇92號礦體西南采區內的采空區及其周邊礦巖作為本次動力失穩仿真技術研究的計算區域。

因空區三維形態復雜，構建復雜空區三維數值計算模型難度大，為此，選用二維 RFPA，選擇典型空區剖面進行分析，相比三維數值分析，二維計算結果安全性更高。因空區主要沿勘探線分布，為選取典型空區進行動力失穩過程模擬仿真，分別沿 205-2，203-2，201-2和202-2這4條勘探線截取剖面，分析剖面上空區動力失穩的過程，其中剖面位置如圖2所示。不同剖面上隱患空區的分布特征見表1。

圖2 隱患空區及計算剖面位置Fig.2 Disaster cavity and section position

2 隱患空區RFPA數值計算模型

2.1 計算模型

目前，RFPA-dynamic尚未開發與其他三維建模軟件或數值軟件數據轉換接口，只能在軟件內建立計算模型，軟件提供了完整的用于建立各種直線、曲線、不規則多邊形、矩形和圓等工具，方便建立各種邊界不規則的隱患空區。因此，充分結合205-2，203-2，201-2和202-2這4條勘探線剖面上空區形態邊界，在RFPA內利用不規則多邊形完成復雜邊界隱患空區數值計算模型的建立，綜合考慮模型尺寸及軟件計算能力，針對不同模型選取不同的尺寸及單元尺寸，具體參數見表2，最后建立的RFPA計算模型如圖3所示。

表1 不同剖面上隱患空區的分布特征Table 1 Distribute characteristic in different sections

表2 RFPA數值計算模型參數Table 2 RFPA numerical calculation model parameters

2.2 基元賦值

RFPA通過 Monte-Carlo方法和統計描述相結合對基元進行初始化賦值[14?15]。設網格中所有基元的彈性模量最大閾值為E0，φ(E)為具有某彈性模量E基元的分布值，則彈性模量Weibull分布函數的積分為

其中：φ(E)為具有彈性模量E的基元的統計量。由上式統計分布構成的基元組成一個樣本空間，在閾值E0不變的情況下，由于m的差別，積分空間分布不完全一樣。這些基元構成的巖石類介質在宏觀上性質可能大體一致(E0相同)，但是，由于細觀結構的無序性，使得基元的空間排列方式有顯著不同。這種細觀上的無序性正好體現了巖石類介質獨特的離散性特征。

圖3 銅坑礦隱患空區RFPA數值計算模型Fig.3 RFPA numerical calculation models of Tongkeng Mine disaster cavity

一般物理空間隨機分布的無序性可以通過Monte-Carlo方法來實現，其產生方法是：產生一組在(0, 1)區間上均勻分布的隨機數序列{iγ≤1|i=1, 2, …,n}；對于任何iγ，存在一個與{iγ≤1|i=1, 2, …,n}相對應的隨機數序列{Ei|i=1, 2, …,n}。由隨機數序列{iγ}映射一組彈性模量參數序列{iγ}。這一組基元彈性模量參數隨機序列逐一賦予網絡中的每一個基元，其他力學參數(強度、泊松比、容重)同樣用此種方法賦值。這種既具有統計性、又具有隨機性(無序性)的基元力學參數賦值方法滿足了網格中基元力學參數非均勻性和隨機性的要求。

根據銅坑礦92號礦體的礦巖的基本條件，采用的基元參數見表3。

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

2.3 初始應力與邊界條件

靜力計算邊界條件為：根據銅坑礦的原巖應力場分布，采用上部受壓、側向受原巖應力場作用的平面應變模型，原巖應力場的水平方向應力為25.4 MPa，垂直方向應力根據模型邊界以上巖石自重求取。

動力計算邊界條件為：在靜力計算邊界條件的基礎上，在空區底板設置爆破地震波。

2.4 爆破震動荷載輸入

將現場實測爆破地震波時間?速度數據進行簡化，并轉換為RFPA能夠接受的應力波應力，如圖4所示。

圖4 爆破震動應力時程曲線Fig.4 Schedule curve of blasting vibration stress

在計算軟件中選擇靜力和動力計算模塊，先設置靜力加載條件，再設置動力加載條件。其中動力加載可以通過應力、位移、波速和波加速度這4種方式實現。本文選擇應力加載方式，在空區底部施加動力荷載。

3 隱患空區動力失穩過程模擬仿真

隱患空區動力失穩仿真的整個計算過程分為靜力計算和動力計算。動力分析以靜力計算為前提，即在完成靜力計算的基礎上，再施加動力時程進行計算。靜力分析時，首先進行初始應力平衡計算，然后，根據空區形成次序分步模擬空區形成過程，每步形成一個空區并計算至平衡狀態。動力計算在靜力計算結束后進行，通過在空區底板施加爆破震動波，再進行動力計算，直至空區破壞，最終獲得空區從形成到失穩破壞的過程模擬仿真。

3.1 勘探線205-2上的空區

圖5所示為勘探線205-2空區動力擾動下失穩破壞過程的剪應力分布。從圖5可見：隨著T506空區的開挖，空區圍巖應力發生重分布，T505空區形成后應力進一步增大，尤其表現在間柱上；施加爆破動載后，空區圍巖受力狀態有所改變，首先在 T506空區右幫出現裂紋；隨著應力波的逐步向上傳播，裂紋不斷往 T505空區方向擴展。由于空區間柱最小寬度僅為5.5 m，裂紋最先在該位置貫通，致使空區間柱完全破壞。

圖5 勘探線205-2剖面空區動力失穩過程仿真Fig.5 Dynamic destabilization numerical simulations of cavity in section 205-2

3.2 勘探線201-2上的空區

圖6所示為勘探線201-2空區動力擾動下失穩破壞過程的剪應力分布。從圖6可見：

(1) 由于S1空區與其他空區間距高達82 m，因此，其他空區的開挖對其影響較小。動荷載作用后，空區左右兩幫周圍圍巖出現局部裂紋，但并未擴展。

(2) T111與T112空間上呈現上下斜交，間柱最小寬度僅為6 m，靜力作用下隨著空區的逐步開挖，空區圍巖應力發生重分布，間柱出現應力集中。

(3) 動荷載作用后，T112空區右幫首先出現局部片幫跨落現象，隨著爆破應力波的進一步作用，片幫范圍逐步加大，并有往下擴展的趨勢。

(4) 動荷載進一步作用后，T110-111頂板靠T112處出現裂紋，擴展與 T112片幫區域連通，最后，整個間柱出現斷裂破壞。

圖6 勘探線201-2剖面空區動力失穩過程仿真Fig.6 Dynamic destabilization numerical simulation of cavity in section 201-2

3.3 勘探線202-2上的空區

圖7所示為勘探線202-2空區動力擾動下失穩破壞過程的剪應力分布。從圖7可見：

(1) 隨著空區的逐步開挖，空區圍巖應力發生重分布，尤其在空區開挖第3步后，由于空區分布密集，空區間柱跨度小，使得空區間柱應力集中。

(2) 空區開挖第 4步后，空區間柱應力集中進一步加強，并出現局部破壞，如T214與T213，T215與T216的間柱均不同程度出現裂紋，此時，裂紋并未擴散。

(3) 動荷載作用后，應力集中程度得到加強，間柱裂紋出現擴展，并沿著間柱薄弱地帶發展，致使T214與T213間柱首先發生破壞。

(4) 動荷載進一步作用后，T215與T216以及T211與 T212的間柱依次出現破壞，致使整個空區向下垮塌。

圖7 勘探線202-2剖面空區動力失穩過程仿真Fig.7 Dynamic destabilization numerical simulation of cavity in section 202-2

3.4 勘探線203-2上的空區

圖8所示為勘探線203-2空區動力擾動下失穩破壞過程的剪應力分布。從圖8可見：

圖8 勘探線203-2剖面空區動力失穩過程仿真Fig.8 Dynamic destabilization numerical simulation of cavity in section 203-2

(1) 模型中T214-316與T312呈上下疊加式，即T314，T315和T316空區的頂板為T312空區的底板，相交處為最危險處，容易出現斷裂。

(2) 空區開挖第 3步后，空區間柱出現明顯的應力集中現象，T314，T315和T316空區頂板出現局部破壞。

(3) 動荷載作用后，空區圍巖應力出現二次分布。T312空區靠T314底板出現裂紋，隨著動荷載的作用，T312空區右幫出現局部破壞。

(4) 動荷載進一步作用后，原有裂紋繼續擴展，最后空區間柱全部垮塌破壞。

4 結論

(1) 在研究 RFPA基本思想和巖石破裂過程仿真原理的基礎上，提出了運用RFPA開展隱患空區動力失穩模擬仿真的方法。

(2) 在對銅坑礦隱患空區現場實測并準確獲取相關信息的基礎上，結合隱患空區危險性可視化分級結果，對處于205-2，203-2，201-2和202-2這4條勘探線上的空區空間分布特征進行了分析，為空區動力失穩過程模擬仿真提供便利。

(3) 模擬了空區從形成到動力擾動作用下的破壞過程，獲得了空區破壞過程的裂紋擴展規律，形成了金屬礦隱患空區動力失穩過程模擬方法。

[1] 何滿潮, 李學元, 劉斌, 等. 深部開采巖體力學研究[J]. 巖石力學與工程學報, 2005, 24(16): 2803?2813.HE Manchao, LI Xueyuan, LIU Bin, et al. Study on rock mechanics in deep mining engineering[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(16): 2803?2813.

[2] 林杭, 曹平, 李江騰, 等. 基于Surpac的Flac3D三維模型自動構建[J]. 中國礦業大學學報, 2008, 37(3): 339?342.LIN Hang, CAO Ping, LI Jiangteng, et al. Automatic generation of Flac3Dmodel based on Surpac[J]. Journal of China University of Mining & Technology, 2008, 37(3): 339?342.

[3] 王濤, 陳曉玲, 楊建. 基于3DGIS和3DEC的地下洞室圍巖穩定性研究[J]. 巖石力學與工程學報, 2005, 24(19): 3476?3481.WANG Tao, CHEN Xiaoling, YANG Jian. Study on stability of underground cavern based on 3DGIS and 3DEC[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(19):3476?3481.

[4] Tang C A, Tang S B. Applications of rock failure process analysis (RFPA) method[J]. Journal of Rock Mechanics and Geotechnical Engineering, 2011, 3(4): 352?372.

[5] 朱萬成, 唐春安, 黃志平, 等. 靜態和動態載荷作用下巖石劈裂破壞模式的數值模擬[J]. 巖石力學與工程學報, 2005, 24(1):1?7.ZHU Wancheng, TANG Chunan, HUANG Zhiping, et al.Numerical simulation on splitting failure mode of rock under static and dynamic loadings[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(1): 1?7.

[6] ZHU Wancheng, TANG Chunan. Numerical simulation on shear fracture process of concrete using mesoscopic mechanical model[J]. Constructions and Building Materials, 2002, 16(8):453?463.

[7] 高峰. 頂板誘導崩落機理及次生災變鏈式效應控制研究[D].長沙: 中南大學資源與安全工程學院, 2009: 20?24.GAO Feng. Study on mechanism of induced caving roof and control of its secondary catastrophe chain effect[D]. Changsha:Central South University. School of Resources and Safety Engineering, 2009: 20?24.

[8] 蘇家紅. 大范圍隱患區資源開采安全控制技術研究[D]. 長沙:中南大學資源與安全工程學院, 2008: 1?2.SU Jiahong. Study on safety control to mining under large scale hidden trouble zone[D]. Changsha: Central South University.School of Resources and Safety Engineering, 2008: 1?2.

[9] LUO Zhouquan, LIU Xiaoming, ZHANG Bao, et al. Cavity 3D modeling and correlative techniques based on cavity monitoring[J]. Journal of Central South University of Technology, 2008, 15(5): 639?644.

[10] LIU Xiling, LI Xibing, LI Faben, et al. 3D cavity detection technique and its application based on cavity auto scanning laser system[J]. Journal of Central South University of Technology,2008, 15(4): 285?288.

[11] LIU Xiaoming, LUO Zhouquan, YUAN Wenni, et al. CMS applications in metal mines of China[C]// Proceedings of Pacific-Asia Workshop on Computational Intelligence and Industrial Application. Wuhan, China, 2009: 353?356.

[12] 劉曉明, 羅周全, 孟穩權, 等. 深井采場大規模垮塌三維探測及可視化計算[J]. 中南大學學報: 自然科學版, 2011, 42(1):158?163.LIU Xiaoming, LUO Zhouquan, MENG Wenquan, et al. 3D monitoring and visible calculation of large-scale collapse area in deep mine[J]. Journal of Central South University: Science and Technology, 2011, 42(1): 158?163.

[13] 劉曉明, 羅周全, 袁雯妮, 等. 基于 CMS的隱患空區三維特征信息獲取[J]. 科技導報, 2011, 29(5): 32?36.LIU Xiaoming, LUO Zhouquan, YUAN Wenni, et al. 3D information acquisition of disaster cavity based on CMS[J].Science & Technology Review, 2011, 29(5): 32?36.

[14] 謝林茂, 朱萬成, 王述紅, 等. 含孔洞巖石試樣三維破裂過程的并行計算分析[J]. 巖土工程學報, 2011, 33(9): 1447?1455.XIE Linmao, ZHU Wancheng, WANG Shuhong, et al.Three-dimensional parallel computing on failure process of rock specimen with a pre-existing circular opening[J]. Chinese Journal of Geotechnical Engineering, 2011, 33(9): 1447?1455.

[15] 黃志平, 唐春安, 馬天輝, 等. 卸載巖爆過程數值試驗研究[J].巖石力學與工程學報, 2011, 30 (S1): 3120?3128.HUANG Zhiping, TANG Chunan, MA Tianhui, et al. Numerical test investigation on unloading rockburst processes[J]. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(S1):3120?3128.