基于Rhino-Griddle-Flac3D的礦體建模與地表穩(wěn)定性分析

石滿生 許 威 張美道

（1.江西天億礦業(yè)有限公司，江西 貴溪 335400；2.江西理工大學(xué)資源與環(huán)境工程學(xué)院，江西 贛州 341000）

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

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

1 工程概況

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

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

2 建立模型

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

2.1 建立礦體模型

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

2.2 建立圍巖模型

將地質(zhì)地形圖導(dǎo)入Rhino，提取等高線。使用“曲線分段”命令將等高線每隔2m劃分成高程點(diǎn)（圖3（a）所示），選中所有高程點(diǎn)后使用“嵌面”命令生成地表模型。再用通過計(jì)算范圍的矩形線框?qū)Φ乇砟Ｐ瓦M(jìn)行布爾分割，得到計(jì)算范圍內(nèi)的地表—圍巖模型。

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

2.3 基于Griddle插件的網(wǎng)格導(dǎo)出

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

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

3 數(shù)值模擬計(jì)算與結(jié)果分析

3.1 邊界條件

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

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

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

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

選用折減后的工程巖體物理力學(xué)參數(shù)如表1。

3.2 工況設(shè)計(jì)

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

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

3.3 計(jì)算結(jié)果分析

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

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

4 結(jié)論

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

（2）數(shù)值模擬結(jié)果表明，天億礦業(yè)開采對地表無顯著不利影響，且采空區(qū)嗣后充填能顯著降低地表變形，提高地表穩(wěn)定性。