基于Rhino-Griddle-Flac3D的礦體建模與地表穩定性分析

石滿生 許 威 張美道

（1.江西天億礦業有限公司，江西 貴溪 335400；2.江西理工大學資源與環境工程學院，江西 贛州 341000）

礦體地下開采后，周圍巖層的地應力會重新分布。隨著采空區暴露面積的增大，若不采取有效措施，開采擾動將會對自身開采范圍內（地表移動帶）的地表建/構筑物以及相臨礦山的安全生產產生一定影響［1］，準確預測和有效預防的技術手段就是開采前進行三維模擬。在運用三維數值模擬計算預測方面，國內外學者進行了較多研究。鄭文棠等［2］基于AutoCAD平臺，借助AutoLisp語言，采用滑動最小二乘法插值擬合，構建三維可視化模型，揭示邊坡的宏觀破壞模式。楊美宏等［3］采用3Dmine軟件對地表及井下工程建立幾何模型，運用Flac3D模擬分析采空區對地表穩定性的影響。李占金等［4］基于FLAC3D對馬城鐵礦深部大采場開采及回填過程中，圍巖與充填體的穩定性進行了數值模擬計算。陳慶發等［5］以3DMine數字化模型為基礎，提出了3DMine-FLAC3D耦合建模方法和3DMine-Surfer-Rhino-ANSYS-FLAC3D多軟件耦合建模方法，構建了鋅多金屬礦三維數值模型，分析了礦體上行開采地表沉陷規律。余璨等［6］基于DIMINE軟件，創建銅廠礦體沿走向、傾向、厚度3個方向上的礦化數學模型，得出Cu品位空間分布模型。程海勇等［7］應用FLAC3D數值分析軟件，探討了龍首礦進路式采礦及充填對地表下沉的影響規律。程立年等［8］采用FLAC3D數值模擬方法對該礦“三下”開采移動范圍內典型剖面的地表沉降進行了分析，得出該礦“三下”開采移動帶內地表建（構）筑物安全穩定。Nengxiong Xu等［9］通過有限差分法（FDM）進行了開采引起的地面沉降預測，以判斷煤層的開采是否會對大壩產生負面影響。Lü Xiangfeng等［10］使用3D有限元代碼對長壁板的地層行為進行了3D建模，并使用CFD代碼對采空區瓦斯氣流進行了模擬，該研究對采礦引起的地層應力變化，裂縫和氣體流型之間的復雜動力相互作用產生了許多新見解。

為確定天億礦業開采對地表穩定性以及銀海礦業的影響，采用“Rhino-Griddle-Flac3D”聯合方法，建立礦體與圍巖的三維數值計算模型，通過計算4種工況和3個變形參數予以分析。

1 工程概況

天億礦業是一座基建期地下礦山，以層狀斑巖型礦床為主，西南方向與銀海礦業、鮑家礦業毗鄰（如圖1）。其中，鮑家礦體處于淺部，南礦帶位于標高+50 m以上、北礦帶位于標高+80 m以上，且通地表。銀海礦業部分礦權與鮑家南礦帶在垂直方向交接，標高+50 m以上為鮑家南礦帶、+50～-70 m為隔離巖體、-70～-420 m為銀海礦權。銀海礦業為生產礦山，當前設計排采計劃至2041年；鮑家礦業也處于生產階段，但資源即將枯竭，故不予考慮。

根據礦山提供的勘探線剖面圖、地表地形圖，地質勘查報告等資料，應用“Rhino-Griddle-Flac3D”三維數值建模與計算方法，建立天億礦業、銀海礦業、鮑家礦業的礦體與圍巖計算模型，探究天億礦業開采對自身開采范圍地表穩定性以及相臨銀海礦業的影響。

2 建立模型

根據圣維南原理，將計算范圍設為（X、Y、Z）從（39 518 500，3 087 800，-850）至（39 521 000，3 090 600，地表），且設（39 518 500，3 087 800，-850）為模型計算相對原點（0，0，0），X軸朝東為正、模型長度2 500 m，Y軸朝北為正、模型寬度2 800 m，Z軸鉛直朝上為正、模型高度從-850 m標高到地表。

2.1 建立礦體模型

將銀海礦業（144～120#勘探線）、鮑家礦業（130～100#勘探線）、天億礦業（100～187#勘探線）勘探線剖面圖導入Rhino，再通過地表XY軸坐標系、標高等作為參考，使用繞軸旋轉、平移等命令調整好各個勘探線剖面圖的位置。在Rhino中，使用“控制點曲線”命令將礦體輪廓圈定出來再使用放樣、曲線、曲面、組合等命令建立礦體相鄰勘探線剖面之間的同一礦體模型，確保每一個礦體模型都為“封閉的實體—多重曲面”。建立的礦體模型如圖2所示。

2.2 建立圍巖模型

將地質地形圖導入Rhino，提取等高線。使用“曲線分段”命令將等高線每隔2m劃分成高程點（圖3（a）所示），選中所有高程點后使用“嵌面”命令生成地表模型。再用通過計算范圍的矩形線框對地表模型進行布爾分割，得到計算范圍內的地表—圍巖模型。

計算范圍內的圍巖主要為上盤花崗斑巖和下盤晶屑凝灰巖。將相鄰的勘探線剖面圖中的地質體分層線使用“線—面”命令生成地質界面。通過地質界面，使用布爾分割命令將計算范圍內的地表—圍巖模型劃分成2種巖性（圖3（b）所示）。

2.3 基于Griddle插件的網格導出

通過Griddle插件將Rhino中網格導入到Flac3D。該插件與陳金龍等［11］應用KUBRIX插件構建大型復雜地質體相比具有以下3個優勢：①Griddle插件中的GInt功能能夠自動生成交互式網格，清理不正確相交的曲面網格，減少因為網格自交而導致導出時的報錯問題（圖4（a）所示）；②如果生成的網格質量差或者未閉合等情況，Griddle插件能夠顯示出出錯的位置，能夠很快對模型進行調整，極大地提高了檢查模型的效率（圖4（b）所示）；③Griddle插件中的Gsurf命令用于將所選的表面網格重新劃分為所需的單元大小和類型，重新劃分后的曲面網格作為輸入數據到Gvol功能中進行體積剖分，使用表面網格來確定實體單元的大小和類型（圖4（c）所示，礦體與圍巖網格長度不同），Gvol命令能填充四面體或六面單元（六面體、棱鏡），生成用于如Flac3D或3DEC的數值模擬程序的文件格式，該功能滿足于本次研究中所需要的“礦體—圍巖網格不同單元大小”情況。通過“Rhino-Griddle-Flac3D”數值建模方法，將導出的GRID文件導入到Flac3D中（圖4d所示）。

注：銀海為黃色模型，鮑家為粉色模型，天億為綠色模型

3 數值模擬計算與結果分析

3.1 邊界條件

假設模型中礦巖體均為理想彈塑性連續介質，且采用摩爾—庫侖（Mohr-Coulomb）屈服準則，其物理力學模型采用彈塑性力學模型，僅考慮自重應力。數值模型初始邊界條件為：前后左右面及底面采用位移約束，即前后左右面限制水平方向位移，底面限制垂直方向位移［12］。

（1）X軸邊界：約束X方向的移動，其邊界位移為0。

（2）Y軸邊界：約束Y方向的移動，其邊界位移為0。

（3）Z軸邊界：約束Z方向下部邊界位移為0，上部設為自由邊界。

選用折減后的工程巖體物理力學參數如表1。

3.2 工況設計

為探究天億礦業開采對計算范圍內地表穩定性及銀海礦業的影響，根據天億礦業和銀海礦業排采計劃，選擇以下4種代表性工況。工況1（最不利情況）：天億礦業和銀海礦業全部礦體開采，不充填；工況2（最有利情況）：天億礦業和銀海礦業全部礦體開采，嗣后全尾砂膠結充填；工況3（天億礦業開采至+10 m，銀海礦業開采至-240 m）：天億礦業和銀海礦業開采到2035年，嗣后全尾砂膠結充填；工況4（天億礦業開采至-40 m，銀海礦業開采結束）：天億礦業和銀海礦業開采到2041年，嗣后全尾砂膠結充填。由于鮑家礦業即將閉坑且礦體通地表，數值模擬中只開采不充填。

在地表40個重點位置設置變形監測點（圖5），具體為：村莊CZ1～CZ3、冷水河HL1～HL4、公路GL1～GL5、銀海工業場地和隔離礦柱YH1～YH17、銀珠山主副井和南北風井等工業場地GC1～GC9、銀珠山斜坡道平硐PD1～PD2。

3.3 計算結果分析

在Flac3D中，根據所設計4種工況進行數值模擬計算。地表沉降情況見圖6，地表變形量匯總見表2。

由表2可知，工況1～工況4的地表40個監測點的最大傾斜i、最大曲率k及最大水平變形ε均小于《有色金屬采礦設計規范》I級保護物、《煤礦測量規程》（2013版）一般磚石結構建筑物規定i=3 mm/m、k=0.2 mm/m/m、ε=2 mm/m的臨界變形極限值。說明天億礦業在工況1～工況4條件下開采，對地表無顯著不利影響。且對比工況1和工況2發現，充填后天億礦業最大沉降減小了316.68 mm、銀海礦業最大沉降減小了65.11 mm，監測點最大傾斜減小了1.138 9 mm/m、最大曲率減小了0.015 8 mm/m/m、最大水平變形減小了0.217 4 mm/m，說明采空區嗣后充填能顯著降低地表變形，提高地表穩定性。

4 結論

（1）Rhino-Griddle-Flac3D數值建模方法具有較多優勢，特別是滿足多種網格單元需求，保證精度的同時提高運算速度，并通過工程實例詳細介紹了這一建模方法。

（2）數值模擬結果表明，天億礦業開采對地表無顯著不利影響，且采空區嗣后充填能顯著降低地表變形，提高地表穩定性。