基于修正Mathews 穩定圖法與FLAC3D 的阿舍勒銅礦深部回采方案優化研究

黃 聰，魏超城，丘永富

（1.長沙礦山研究院有限責任公司，湖南 長沙 410012；2.金屬礦山安全技術國家重點實驗室，湖南 長沙 410012；3.新疆哈巴河阿舍勒銅業股份有限公司，新疆 哈巴河 836700；4.新疆金寶礦業有限公司，新疆 富蘊 836199）

0 引言

經過十多年的開采，阿舍勒銅礦淺部資源逐漸枯竭，礦山逐步轉入深部開采階段。受到深部“三高一擾動”的復雜力學環境影響，巖體的力學特征及其工程也相應發生了明顯變化。巖爆、冒頂、采空區失穩等地質災害程度加劇，對礦山深部資源的開采產生了巨大威脅。因此，選擇合理的采場結構參數和回采順序對于控制圍巖變形、預防采空區失穩有著積極作用。

近年來，眾多專家學者對采場結構參數優化與回采順序選擇進行了深入研究。周科平[1]將遺傳算法應用于采場結構參數的優化，并闡述了其基本原理。彭康等[2]通過工業試驗驗證了使用數值模擬來確定采場參數的可行性。趙永等[3]利用Mathews 穩定圖法對赤峰紅嶺鉛鋅礦的采空區進行了穩定性及穩定概率評價。徐荃等[4]采用數值模擬與Mathews穩定圖法相結合的方式對某金礦深部采場布置形式和結構參數進行了優化研究。孫健等[5]針對地下礦山采場結構參數應用Mathews 圖解法對急傾斜采場暴露面冒落問題進行了研究。尹土兵等[6]認為采場暴露面積的確定直接影響采場結構參數的合理布置，采用階段空場嗣后充填采礦法時，需要同時討論采場頂板和充填體側幫的暴露面積。特魯曼等[7]收集了大量工程實例，與加拿大數據庫合并分析后，提出了一種修正的Mathews 圖解法用于預測采場幫壁穩定性。

以上研究對于礦山回采方案優化具有一定的指導意義，但由于阿舍勒礦區構造復雜，并不能直接套用。因此，使用Mathews 穩定圖法與數值模擬相結合的方法，以阿舍勒銅礦深部采場參數、回采順序等為優化對象，分析不同回采順序與采場結構條件下，開采范圍內的礦巖介質應力、位移、塑性區的變化過程，以期為礦山生產提供依據。

1 工程概況

阿舍勒銅礦礦體傾角為55o～75o，礦體產狀總體走向近南北，傾向東，分布范圍為水平方向3#勘探線～11#勘探線，空間范圍+300～+50 m，其礦體賦存條件如圖1 和圖2 所示。目前阿舍勒銅礦的深部生產中段主要為+400～0 m 共9 個中段，中段高度為50 m，采用大直徑深孔空場嗣后充填采礦法開采，采場垂直礦體走向布置，長度為礦體厚度，采場寬度20 m，開采順序為由北向南隔一采一。

圖1 +150 m 中段地質平面圖Fig.1 Geological plane of +150 m middle section

圖2 5#勘探線剖面圖Fig.2 Profile of exploration line 5

阿舍勒銅礦深部開采目前存在貧化損失率居高不下、采場穩定性較差以及地壓控制困難等諸多難題，選取具有代表性的+150 m 中段作為研究對象，旨在通過優化采場結構、選擇合理開采順序，達到控制地壓、降低貧化損失率的目的，實現礦山安全高效生產。

2 基于修正Mathews 穩定圖法采場尺寸確定

Mathews 圖解法是MATHEWS 等基于大量工程實例在1980 年提出的一種礦山開采設計方法，最初應用于對采場穩定性進行預測。經過大量學者的研究與補充，目前Mathews 圖解法已成為一種相對簡單的巖體分級系統，被廣泛應用于國外礦山采場設計。根據阿舍勒銅礦深部采場所使用的采礦方法，采用修正Mathews 穩定性圖解方法分別分析采場頂板和采場上盤圍巖區域的巖體穩定性[8]。

2.1 修正Mathews 穩定圖法

Mathews 穩定性系數的計算公式見式（1）。Mathews穩定性系數與水力半徑的相關關系如圖3 所示。

圖3 Mathews 穩定性系數與水力半徑關系Fig.3 Relationship between Mathews stability coefficient and hydraulic radius

式中：N為Mathews 穩定性系數；Q'為修正的Q值；A為巖石應力系數；B為節理方位系數；C為重力調整系數。

應用Mathews 穩定性圖表方法需要計算兩個參數，即穩定性系數和形狀因素，其中，穩定性系數反映了在一定的應力條件下巖體自立的能力；形狀因素考慮了單獨采場暴露表面的尺寸和形狀，形狀系數是暴露面的水力半徑。

Mathews 穩定性圖解方法采用了修正的NGI 隧道質量指標Q，與Q值不同的是，Q'值中的應力折減系數（SRF）和節理滲水折減系數（Jw）都為1.0，Q'值的計算公式見式（2）。

A值考慮高應力影響降低巖體穩定。A值為完整巖體的單軸抗壓強度與平行開挖面的最大誘導應力的比值。A值與 σc/σi成線性關系，變化范圍為0.1～1.0。巖石應力系數A的圖解如圖4 所示。

圖4 巖石應力系數A 的圖解Fig.4 Diagram of rock stress coefficient A

B值要考慮不連續面的方向影響，根據控制性節理與采場表面的相對方位確定，結構面與開挖面的夾角為90°時，B系數被賦值1，不連續結構面與開挖表面的夾角為20°時，B值為0.2，節理方位系數B的圖解如圖5 所示。

圖5 節理方位系數B 的圖解Fig.5 Diagram of joint azimuth coefficient B

受重力影響，采場頂板穩定性小于側幫，重力調整系數C考慮了重力對采場暴露表面崩落、滑落等穩定性影響，重力調整系數C和采場表面傾角β的關系由式（3）確定。重力調整系數C的圖解如圖6所示。水力半徑用于分析表面積除以暴露墻周長的比值，如圖7 所示，計算公式見式（4）。

圖6 重力調整系數C 的圖解Fig.6 Diagram of gravity adjustment coefficient C

圖7 水力半徑確定方法的圖解Fig.7 Diagram of hydraulic radius determination method

2.2 阿舍勒銅礦采場結構參數優化

2.2.1 穩定性系數和水力半徑

修正的Q值的計算結果見表1。巖石應力系數A見表2。礦山深部以厚大、急傾斜礦體為主，礦體傾向東，傾角60°～80°。其上盤圍巖主要為黃鐵礦，礦體為銅硫礦，下盤圍巖主要為凝灰巖。B值根據控制性節理與采場表面的相對方位確定，見表3。

表1 修正的Q 值計算結果Table 1 Calculation results of corrected Q-value

表2 巖石應力系數ATable 2 Rock stress coefficient A

表3 節理方位系數BTable 3 Joint azimuth coefficient B

受重力的影響，采場頂板的穩定性小于側幫，重力調整系數C考慮了重力對采場暴露表面崩落、滑落等穩定性的影響。

1）采場頂板。采場頂板為水平狀即β=0?，重力調整系數C計算見式（5）。

2）采場上盤。采場上盤傾角一般為70°，即β=70?，重力調整系數C計算見式（6）。

根據求得的Mathews 穩定性系數計算參數，按照式（1）計算得出采場頂板和采場上盤巖體的穩定性系數，運用圖解法得出與穩定性系數對應的容許水力半徑取值區間，結果見表4。

表4 穩定性系數計算結果Table 4 Calculation results of stability coefficient

2.2.2 采場頂板穩定性分析

礦山深部主要采用大直徑深孔采礦法，采場垂直礦體走向布置，中段高度為50 m，礦體最大厚度為50 m，回采規劃中采場寬度以12 m 為主，根據水力半徑計算方法，計算得出不同采場寬度和長度條件下頂板的水力半徑，計算結果如圖8 所示。

圖8 不同采場結構下頂板的水力半徑Fig.8 Hydraulic radius of roof under different stope structures

在對深部巖體質量指標Q值計算的基礎上，根據Mathews 圖解法，在中段高度一定的條件下，對不同采場寬度和長度尺寸下采場頂板的容許水力半徑進行計算分析，得出水力半徑值≤4.8 m 對應的采場尺寸，其采場頂板處于無支護臨界穩定區，整體穩定狀態較好，滿足采場安全生產要求。

2.2.3 采場上盤穩定性分析

由于深部開采中段高度固定為50 m，礦體最大厚度為50 m，上盤圍巖的傾角平均為70°，根據水力半徑計算公式，分別計算不同采場寬度尺寸下，上盤圍巖的水力半徑（以上盤斜長計算），計算結果見表5。

表5 不同采場上盤圍巖的水力半徑Table 5 Hydraulic radius of hanging wall surrounding rock in different stopes

結合表4 采場上盤圍巖穩定性系數對應的容許水力半徑值可知，當采場上盤的水力半徑≤4.9 m，上盤圍巖處于穩定區或無支護過渡區，則滿足采場安全生產需要。即在采場高度一定的條件下，采場寬度≤12 m，可以保證礦體最厚處（礦房長度50 m）的采場上盤圍巖整體穩定狀態較好，滿足采場安全生產要求。

3 深部開采回采順序優化

根據分段空場嗣后充填法的特點，綜合利用了CAD、Midas GTS、FLAC3D建立深部礦體三維數值模型，基于深部開采的實際情況、前期開展的工程地質調查和室內巖石力學實驗結果，對不同的開采方案進行了數值模擬。在地質條件、采場結構尺寸一定的條件下，使用FLAC3D軟件對各開采順序方案進行模擬計算，比較采場附近的應力和位移，分析各回采方案的優劣，從而得到合理的回采順序[9-14]。

3.1 數值模擬模型建立

深部+150 m 中段是研究對象，在網格劃分時將150 m 礦體部分網格細化，并適當控制最大單元和最小單元尺寸間的比例以防止計算出現較大的偏差。

計算模型x方向563 m，y方向315 m，z方向250 m，單元大小礦體2m×5m×5m，圍巖4m×5m×5m，模型單元總數31.2 萬單元，計算模型如圖9所示。

圖9 數值模擬模型Fig.9 Numerical simulation model

3.2 模擬方案

阿舍勒銅礦深部采用從上至下的中段開采順序，在中段開采順序和采場結構參數確定的條件下進行采場回采順序的模擬分析。根據礦床開采技術條件和大直徑深孔空場嗣后充填法的特點，采場垂直礦體走向布置，礦房長度為礦體厚度，中段高度50 m，為分析不同的采場回采順序對采場穩定性的影響，以+150 m 中段為模擬對象，共設計了四種不同的回采方案[15]（圖10 和表6），即從中央向兩端推進回采（方案1）、從兩端向中央推進回采（方案2）、從北端向南端推進回采（方案3）、從南端向北端推進回采（方案4），每種回采方案礦房、礦柱的回采方式均為隔三采一。采場編號如圖10 所示。

表6 模擬方案詳情Table 6 Details of simulation scheme

圖10 采場編號示意圖Fig.10 Schematic diagram of stope number

3.3 計算結果及比較分析

對于數值模擬計算的結果，主要從應力分布狀態、位移變化狀態、塑性區分布狀態等三個方面進行對比分析。

3.3.1 采場應力分布狀態的比較分析

根據各方案數值模擬計算的結果，全部采場回采結束后的最大主應力等值線分布圖如圖11 所示。由圖11 可知，在采場回采過程中，最大主應力集中區主要出現在兩側幫中部，最大主應力值-1.62～-0.95 MPa，側幫巖體應力表現為受壓狀態，應力集中區域在采場上盤，采場易產生下盤巖體受壓底鼓及上盤巖體受壓冒落的現象。

圖11 采場最大主應力等值線圖Fig.11 Contour map of maximum principal stress in stope

根據數值模擬計算結果，對各方案中采場在回采過程中產生的最大主應力值進行統計，采場最大主應力值統計見表7。由表7 可知，方案4 中各采場的最大主應力值均小于其他幾種方案，表明方案4中采場的最大主應力處于較低的應力狀態。由各采場最大主應力變化趨勢可以看出，方案1、方案2 中采場最大主應力值變化區間較大，部分采場集中應力值較大，方案3、方案4 中采場最大主應力值變化趨勢相似，變化區間較小。由此表明方案1、方案2中局部采場應力集中突出現象明顯，方案3、方案4中采場壓應力變化較為平緩，應力在采場回采中分布較為均勻。

表7 采場最大主應力值Table 7 Maximum principal stress value of stope

3.3.2 采場位移變化狀態的比較分析

根據數值模擬的計算結果，各方案采場回采結束后所產生的位移等值線分布圖如圖12 所示。由圖12 可知，各采場的最大位移均位于采場的兩側幫和上盤巖體，側幫位移值為19.9～44.9 mm，上盤巖體位移值為16.2～51.9 mm，頂板產生的位移則較小，位移值為4.6～37.8 mm（表8～表10）。

表8 采場側幫最大位移值Table 8 Maximum displacement value of stope side slope

表9 采場上盤巖體最大位移值Table 9 Maximum displacement value of hanging wall rock mass in stope

表10 采場頂板最大位移值Table 10 Maximum displacement of stope roof

圖12 采場位移等值線圖Fig.12 Contour map of stope displacement

由表8～表10 可知，方案1 中回采采場在中間時側幫、上盤和頂板位移值均較小，但隨著采場向兩端的推進，位移值明顯增加，尤其端部采場各區域位移值較大，方案1 中各采場上盤、頂板的位移較其他方案均較大。方案2、方案3 各采場側幫、上盤和頂板位移值較接近，在中部和端部的采場中位移變化較方案4 偏大。方案4 中采場位移較其他方案均較小，表明方案4 中采場位移變化較小，有利于采場的穩定，從位移變化的角度考慮，方案4 明顯優于其他方案。

3.3.3 采場塑性區分布狀態的比較分析

根據數值模擬的計算結果，各方案采場回采結束后所產生的塑性區分布圖如圖13～圖16 所示。由圖13～圖16 可知，方案1 和方案2 中采場充填體的塑性區面積較接近，呈現大面積的貫通性拉張破壞，上盤和頂板巖體表現為剪切破壞，且方案2 中上盤巖體剪切破壞面積明顯大于方案1，采場側幫以拉張破壞為主，這與充填體強度較小有較大的關系，采場充填體較易受拉垮落。

圖13 方案1 采場回采結束后塑性區分布圖Fig.13 Distribution map of plastic zone after stope mining in scheme 1

圖14 方案2 采場回采結束后塑性區分布圖Fig.14 Distribution map of plastic zone after stope mining in scheme 2

圖15 方案3 采場回采結束后塑性區分布圖Fig.15 Distribution map of plastic zone after stope mining in scheme 3

圖16 方案4 采場回采結束后塑性區分布圖Fig.16 Distribution map of plastic zone after stope mining in scheme 4

方案3、方案4 中采場充填體塑性區面積相近，采場充填體的塑性區面積小于方案1、方案2，且塑性破壞區的也未完全貫通，這一結果有利于改善充填體內的受力狀態，減少充填體的破壞程度，增強礦體二步驟回采時充填體的自穩能力。所以從塑性區分布的角度考慮，方案3、方案4 明顯優于其他方案。

通過從應力分布狀態、位移變化狀態和塑性區分布狀態三個方面對上述四種模擬方案進行對比分析可知，方案4 最優，其次為方案3。因此，建議深部采場中段回采順序為由南端向北端推進回采。

4 結語

采用Mathews 穩定性圖解法，根據深部采場的巖體特征，對采場頂板和上盤巖體的穩定性進行了計算分析。當深部采場的中段高度為50 m，采場長度為礦體最大厚度時，采場寬度為12 m 即可處于穩定狀態。根據水力半徑計算結果，只要采場長度小于50 m，采場頂板和上盤圍巖處于無支護臨界穩定區，整體穩定狀態較好，滿足采場安全生產要求。

通過三維數值模擬計算方法，對深部采場的回采順序進行了模擬對比分析，共設計了四種回采方案，即從中央向兩端推進回采、從兩端向中央推進回采、從南端向北端推進回采、從北端向南端推進回采，對采場回采過程中的應力分布狀態、位移變化特點和塑性區特征進行了對比分析，方案4（從南端向北端推進回采）最優，方案3 次之。因此，建議深部采場中段回采順序為由礦體南端向北端推進依次回采。