基于Flac3D的某金礦回采方案優化數值模擬研究

王磊 閆剛 付佳杰 鐘健

收稿日期：2023-12-26； 修回日期：2024-02-15

基金項目：國家自然科學基金青年項目（52004053）；中央高?；究蒲袠I務專項資金資助-國防重大培育項目（N2101040）

作者簡介：王? 磊（1980—），男，工程師，從事礦山安全管理工作；E-mail：hdjkahk@163.com

*通信作者：鐘? ?。?999—），男，碩士研究生，從事地表巖層移動機理研究工作；E-mail：1184729590@qq.com

摘要：為了探究礦床開采順序對地表沉陷的影響，以某金礦的生產實際問題為工程背景，擬定不同開采方案，并采用Flac3D軟件數值模擬的方法，將地表監測的變形量及塑性區體積作為評價指標，分析不同回采方案下的地表移動規律。研究結果表明：2個中段同時開采產生的擾動會對地表移動有較大影響，第七中段與第八中段同時回采時，地表沉陷量最小；塑性區破壞體積與開采順序有較大的關系，同時開采的中段總體積越小，礦山整體開采結束后被破壞的塑性區體積越小。研究結果為礦山安全高效開采提供參考依據。

關鍵詞：充填采礦；回采順序；數值模擬；地表移動；地表沉陷

中圖分類號：TD853.34????????? 文章編號：1001-1277（2024）06-0011-05

文獻標志碼：Adoi：10.11792/hj20240603

引? 言

礦產資源開采引起的地表沉陷給后續的生產與安全帶來了一系列消極影響，在開采過程中合理的回采—充填過程能夠有效減緩地表的沉陷，不僅能夠降低礦石回采的成本，更能提高礦山采礦效率與安全水平。國內外學者利用數值模擬的方式對回采順序優化進行了大量研究，取得了豐碩的成果［1-8］。楊八九等［9］基于數值模擬結果發現部分開采區域有應力集中現象，安全隱患極大，通過多方案的采場結構參數及回采順序論證，優化了后續開采方案，提高了開采安全系數。楊海彬等［10］針對某金礦深部開采地壓顯現問題，對不同礦房、礦柱跨度與分層采高的方案進行模擬，分別從應力、位移、塑性區進行分析，確定了最優采場結構參數優化方案。郭進平等［11］為分析四方金礦深部礦體由空場采礦法轉崩落采礦法回采時的地壓活動規律，對不同回采方案的位移、應力、塑性區分布規律進行了研究，確定了地壓活動最小的最優回采方案。張瑞明等［12］在不同采場跨度和不同回采順序條件下對采場進行了穩定性模擬分析，確定了最優的采場結構參數及回采順序。吳振坤等［13］通過正交試驗結果建立了多指標綜合評價模型，在綜合滿意度最優條件下確定了最優回采順序與采場跨度。陳順滿等［14］建立了多指標的響應面回歸模型，基于整體滿意度準則利用非線性目標優化算法對采場結構參數進行優化。GUO等［15］運用理論分析法對復雜傾斜薄礦體的采場結構參數進行了優化。BAI等［16］提出了網絡流的方法，提高了采場結構參數的優化效率。

某金礦主要有Ⅰ、Ⅱ、Ⅲ號礦體，估算礦石資源儲量達30.387 6萬t，平均金品位1.96 g/t。礦體總體走向NE50°，傾向NW，傾角約30°，礦塊平均長度126 m，平均水平厚度18.5 m，礦體總體呈大透鏡狀，賦存標高為-420～-380 m。根據礦體賦存特征、開采技術條件和礦山開采現狀，采用上向水平分層充填采礦法回采礦體。根據礦體賦存情況布置礦房，礦塊長20～40 m，寬為礦體厚度，高為中段高30～40 m，分段高10 m，分層高2.5 m，底柱高5 m，不留頂柱和間柱。爆破通風后進行頂板和上盤排險工作，局部不穩固地段采用錨桿金屬網支護，錨桿長度一般為1.5～2.0 m，安裝網度為1.0 m×1.0 m～0.5 m×0.5 m，與錨桿同時安裝。

該金礦Ⅰ、Ⅱ號礦體賦存高度差較大且礦體走向相似，在開采設計時可一同考慮設計開采第一至第七中段。Ⅲ號礦體距Ⅰ、Ⅱ號礦體較遠，賦存高度差較小，在回采時需單獨考慮設計開采第八中段，Ⅲ號礦體中段可與Ⅰ、Ⅱ號礦體的其中一個中段同時回采。若回采過程中采動疊加影響過大，會影響地表建（構）筑物等的安全，而按照中段順序進行逐中段開采會導致礦區生產能力小、成本高等問題。為降低采動疊加效應對地表的影響，保障礦山安全、高效生產，需對中段回采順序進行研究。本文為探究Ⅲ號礦體中段的最佳開采順序，利用Flac3D軟件數值模擬分析礦體不同回采方案的回采—充填過程對地表移動規律的影響，對礦山回采順序進行優化，以提高回采順序的合理性，為安全、高效回采提供保證。

1? 礦區模型建立及回采方案設計

1.1? 計算區域與網格劃分

根據巖體構造特點和礦區地形地貌，建立三維數值模型（見圖1）。根據圣維南原理，圍巖區域應取采空區范圍的3～5倍來消除局部效應的影響［17］。本文三維數值模型研究的尺寸為長×寬×高=1 200 m×1 100 m×700 m。模型共分1 183 137個四面體單元，984 844個節點。

1.2? 邊界條件和初始條件

本次模擬邊界條件采用位移約束，即模型的前后、左右邊界施加水平約束，邊界水平初始位移為零；模型的底面邊界施加水平和垂直約束，邊界水平和垂直初始位移為零；模型頂面邊界為自由面，不施加約束條件。計算初始條件考慮自重應力場，重力加速度取9.81 m/s2，對模型進行計算，平衡后將位移場、速度場、塑性區清零。

1.3? 巖體力學參數和本構模型

巖石是一種脆性材料，當荷載達到屈服強度后將發生破壞、弱化，屬于彈塑性體。本次計算模擬的礦體及圍巖本構模型采用莫爾-庫侖彈塑性模型，其破壞準則為［18］：

fs=σ1-σ31+sin φ1-sin φ-2C1+sin φ1-sin φ （1）

式中：σ1、σ3分別為最大和最小主應力（MPa）；C、φ分別為內聚力（MPa）和內摩擦角（°）。

當fs>0時，材料將發生剪切破壞。在通常應力狀態下，巖體的抗拉強度很低。因此，可根據抗拉強度準則判斷巖體是否產生拉破壞。

數值模擬模型的礦巖及充填體物理力學參數見表1。

1.4? 回采方案設計

為確定Ⅲ號礦體的最佳回采順序，根據礦區礦體實際情況，采用自上而下的下行式開采，模擬Ⅲ號礦體第八中段分別與前7個中段同時回采的7種不同回采方案，見表2。各中段位置關系見圖2。

2? 回采方案優化分析

把地表監測的變形量及塑性區體積作為評價指標對各個方案進行評價，對不同模擬結果進行分析，計算結果見表3。

2.1? 地表位移分析

不同方案地表總位移云圖見圖3。由圖3可知：方案5的地表總位移最大，為44.847 mm，方案7的地表總位移最小，為43.350 mm。這說明，第五中段與第八中段同時開采對地表沉降的影響最大，第七中段與第八中段同時開采對地表沉降的影響最小。在實際生產中，每個中段開采后的采空區體積均不同，在模擬計算中，每個中段的橫截面積設置為相同尺寸，僅采場長度不同。前7個中段所采礦體體積中第五中段體積最大，為8.4×105 m3；第七中段體積最小，為4.5×105 m3。在模擬方案中，2個中段同時開采產生的擾動會對地表移動有較大影響，其中體積最大的第五中段與體積最小的第七中段分別對應了地表沉陷量最大值與最小值。結合實際中段采空區體積，對比不同方案下的地表位移情況，方案7中第八中段與采空區體積最小的第七中段同時開采，產生的擾動影響最低，對地表的保護作用最好。

2.2? 塑性區分析

不同方案中段采空區與塑性區體積分布見圖4。由圖4可知：塑性區體積與開采順序有較大的關系。其中，方案5開采后的塑性區體積最大，為11.46×105 m3，方案7開采后的塑性區體積最小，為7.9×105 m3。同時，開采的中段體積越小，礦山整體開采結束后被破壞的塑性區體積越小。結合塑性區分布圖（見圖5），從整體上看，被破壞的塑性區體積較小，僅在第三中段產生了較多的塑性區域，沒有出現大面積的塑性貫通現象。因此，可以認為礦山開采后的圍巖穩定。為了減少第三中段被破壞的塑性區體積，在第三中段充填時，可以適當提高充填體強度，提高第三中段的充填體承載能力，保障后續中段的安全生產。

2.3? 地表移動安全性綜合評價

對于地表移動的安全性評價，反映地表變形沉降和移動特征的指標除了地表垂直沉降位移和水平移動位移，一般還有地表的傾斜變形、曲率、水平變形［19］。

1）地表傾斜變形。相鄰點在豎直方向的相對移動量與相鄰兩點間水平距離的比值。它反映地表沿某一方向的坡度，通常以i表示。

iAB=WB-WAlAB=ΔWABlAB（2）

式中：iAB為傾斜值（mm/m）；WA、WB分別為地表觀測點A、B的下沉值（mm）；lAB為地表觀測點A、B點間的水平距離（m）。

2）地表曲率。兩相鄰線段的傾斜差和兩線段中點間水平距離的比值。它反映了觀測線斷面上的彎曲程度，以K表示，地表下沉曲線上凸為正、下凹為負。

KB=iBC-iAB12（lBC+lAB）=2ΔiA-B-ClBC+lAB（3）

式中：KB為曲率（10-3/m）；iAB、iBC分別為地表A、B點間和B、C點間的平均斜率（mm/m）；lAB、lBC分別為地表A、B點間和B、C點間的水平距離（m）。

3）地表水平變形。相鄰兩點的水平移動差值與兩點間水平距離的比值。它反映相鄰兩觀測點間單位長度的水平移動差值，通常以ε表示。水平變形正值表示拉伸，分布在移動盆地拐點和邊界；負值表示壓縮，分布在移動盆地兩拐點之間。

εAB=UB-UAlAB=ΔUABlAB（4）

式中：εAB為水平變形值（mm/m）；UA、UB分別為地表觀測點A、B的水平移動值（mm）。

對地表沉陷量最小及被破壞塑性區體積最小的方案7的地表監測點位移規律進行分析，結果見圖6。由圖6可知：開采充填后的地表最大傾斜變形i=-0.735 mm/m，最大水平變形ε=-0.316 mm/m，最大曲率K=0.037×10-3/m，均符合GB 50771—2012 《有色金屬采礦設計規范》允許的地表變形值。

3? 結? 論

1）在模擬方案中，兩個中段同時開采產生的擾動會對地表移動有較大影響，第七中段與第八中段同時回采時，地表沉陷量最小，對地表的保護作用最好。

2）塑性區體積與開采順序有較大關系，同時開采的中段體積越小，礦山整體開采結束后被破壞的塑性區體積越小。為了減少被破壞的塑性區體積，在充填時可以適當提高充填體強度及承載能力。

3）基于最優回采順序的模擬結果，得出地表最大傾斜變形i=-0.735 mm/m，最大水平變形ε=-0.316 mm/m，最大曲率K=0.037×10-3/m，均符合GB 50771—2012 《有色金屬采礦設計規范》允許的地表變形值，地下開采活動對地表河流及建（構）筑物影響在安全范圍內。

4）礦山采用優選出的第七、第八中段回采順序方案，極大地降低了礦山開采過程中的安全風險，并提高了回采率，為礦山企業帶來顯著的經濟效益。

［參 考 文 獻］

［1］? 姜永恒.鎮沅金礦無底柱分段崩落采礦法采場結構參數的確定［J］.黃金，2022，43（5）：37-42.

［2］? 劉濱，李振理，尹旭巖，等.上向水平分層充填采礦法在新城金礦的應用［J］.黃金，2018，39（1）：48-51.

［3］? 崔曉林.三鑫金銅礦深部緩傾斜薄礦體回采工藝的優化選擇［J］.黃金，2023，44（1）：33-37.

［4］? 國紹林，宋召法，劉濱，等.踏步式上向水平分層充填采礦法回采步距優化及應用［J］.黃金，2019，40（1）：26-29.

［5］? 周亞博.基于Flac3D數值模擬的深部采場結構參數優化及應用［J］.黃金，2022，43（5）：51-55.

［6］? 高慶偉，孫嘉，郝顯福.基于3D-σ有限元法的采場結構參數優化［J］.黃金，2017，38（6）：29-32.

［7］? 王亞軍，李向東，盛佳.前河金礦Ⅳ S3礦體280 m中段回采順序優化［J］.礦業研究與開發，2018，38（5）：13-15.

［8］? 關凱，朱萬成，張洪訓，等.新城金礦中段頂柱盤區回采順序優化［J］.金屬礦山，2013（9）：44-49.

［9］? 楊八九，侯克鵬，孔學偉，等.基于數值模擬分析的采場應力阻隔及地壓控制研究［J］.有色金屬工程，2021，11（4）：118-124.

［10］? 楊海彬，李文超，秦秀合，等.某金礦深部采場結構參數與回采順序優化研究［J］.礦業研究與開發，2020，40（12）：18-22.

［11］? 郭進平，強浩，劉少青，等.深部礦體無底柱分段崩落法回采順序優化研究［J］.礦業研究與開發，2020，40（9）：11-17.

［12］? 張瑞明，魏丁一，杜翠鳳，等.階段空場嗣后充填采礦法采場跨度及回采順序優化研究［J］.礦業研究與開發，2022，42（3）：10-14.

［13］? 吳振坤，金愛兵，陳帥軍.基于多指標綜合評價的采場結構參數及回采順序優選［J］.礦業研究與開發，2020，40（2）：1-7.

［14］? 陳順滿，吳愛祥，王貽明，等.基于響應面法的破碎圍巖條件下采場結構參數優化研究［J］.巖石力學與工程學報，2017，36（增刊1）：3 499-3 508.

［15］? GUO Y，HOU K，LI W.Numerical optimization of stope structural parameters of complex inclined thin orebody in Mengnuo lead-zinc mine［J］.Electronic Journal of Geotechnical Engineering，2014，19：9 479-9 490.

［16］? BAI X，MARCOTTE D，SIMON R.Underground stope optimization with network flow method［J］.Computers & Geosciences，2013，52（3）：361-371.

［17］? 杜坤，李夕兵，劉科偉，等.采空區危險性評價的綜合方法及工程應用［J］.中南大學學報（自然科學版），2011，42（9）：2 802-2 811.

［18］? 韓斌，吳愛祥，劉同有，等.金川二礦區多中段機械化盤區回采順序的數值模擬優化研究［J］.礦冶工程，2004（2）：4-7.

［19］? LI H，GUO G，ZHAI S C.Mining scheme design for super-high water backfill strip mining under buildings：A Chinese case study［J］.Environmental Earth Sciences，2016，75（12）：1 017.

Numerical experimental study on optimization

of mining scheme in a gold mine based on Flac3D

Wang Lei1，Yan Gang1，Fu Jiajie1，Zhong Jian2

（1.Canzhuang Gold Mine，Zhaojin Mining Industry Co.，Ltd.;

2.School of Resources and Civil Engineering，Northeastern University）

Abstract：In order to explore the effect of mining sequence on surface subsidence，taking the practical production problems of a gold mine as the engineering background，different mining schemes were designed.The Flac3D software was used for numerical simulation，with surface deformation magnitude and plastic zone volume monitored as evaluation indicators to analyze the movement patterns of the surface under different mining schemes.The research results indicate that disturbances generated by simultaneous mining of two levels have a significant impact on surface movement，and the smallest surface subsidence occurs when the seventh and eighth levels are simultaneously mined.The volume of plastic zone damaged is closely related to the mining sequence，and the smaller the total volume of levels mined simul-taneously，the smaller the volume of plastic zones damaged after the completion of overall mining in the mine.The research findings provide a reference basis for safe and efficient mining in the mine.

Keywords：filling mining;mining sequence;numerical simulation;surface movement;surface subsidence