康育鵬，孟海梅

(河南理工大學(xué)測(cè)繪與國(guó)土信息工程學(xué)院，河南焦作454003)

煤炭是我國(guó)的主要能源，在我國(guó)的一次能源消耗中占65%左右，我國(guó)的煤炭資源雖然豐富，但經(jīng)過(guò)長(zhǎng)期大量開采正在逐年減少，尤其是在我國(guó)的東部地區(qū)，煤炭資源正在逐步枯竭，剩余儲(chǔ)量的一半以上屬于“三下”(建筑物下、水體下、鐵路公路下)壓煤。據(jù)不完全統(tǒng)計(jì)，全國(guó)生產(chǎn)礦井“三下”壓煤量達(dá)137.9億t，其中建筑物下壓煤87.6億t，村莊下壓煤占建筑物下壓煤的60%以上［1］。煤礦“三下”壓煤?jiǎn)栴}是當(dāng)前和今后一段時(shí)期內(nèi)煤炭行業(yè)發(fā)展面臨的主要問題之一，關(guān)系到煤炭資源的合理開發(fā)與礦井的可持續(xù)發(fā)展。

固體充填開采是近些年發(fā)展起來(lái)的一種綠色開采技術(shù)，它是利用矸石等固體充填物充入采空區(qū)，占據(jù)上覆巖層的下沉空間，以阻止頂板的下沉［2－3］。相比于傳統(tǒng)的垮落法開采，它可以有效減小采動(dòng)損害，控制上覆巖層與地表的移動(dòng)變形，大大減少矸石等固體廢棄物的排放，進(jìn)而解決煤礦的“三下”壓煤?jiǎn)栴}。

由于使用FLAC3D等數(shù)值模擬軟件進(jìn)行開采沉陷的模擬與概率積分法得到的預(yù)計(jì)結(jié)果比較相近，所以FLAC3D等數(shù)值模擬軟件可以被廣泛應(yīng)用于開采沉陷的預(yù)計(jì)與可視化方面［4］。固體充填開采雖然可以有效減小采動(dòng)損害，但仍不可避免地造成一定程度的地表沉陷，因此需要對(duì)固體充填開采所引起的地表沉陷進(jìn)行預(yù)測(cè)。本文嘗試使用FLAC3D軟件對(duì)固體充填開采進(jìn)行數(shù)值模擬研究，并應(yīng)用其內(nèi)嵌的FISH程序語(yǔ)言，編制相關(guān)程序以提取出固體充填開采地表的沉陷數(shù)據(jù)，結(jié)合MATLAB強(qiáng)大的數(shù)據(jù)可視化功能完成對(duì)固體充填開采地表沉陷的可視化操作。

1 數(shù)值模擬與模型建立

1.1 數(shù)值模擬及FLAC3D簡(jiǎn)介

數(shù)值模擬也叫計(jì)算機(jī)模擬，它是通過(guò)數(shù)值計(jì)算和圖像顯示的方法，借助計(jì)算機(jī)來(lái)獲得滿足工程要求的數(shù)值解，對(duì)工程問題和自然界各類問題進(jìn)行研究。數(shù)值模擬現(xiàn)在已經(jīng)成為礦山開采沉陷研究的重要手段之一。

FLAC3D是三維快速拉格朗日分析(Fast Lagrangian Analysis of Continua)的簡(jiǎn)稱，F(xiàn)LAC3D程序自美國(guó)ITASCA咨詢集團(tuán)推出后，已經(jīng)成為目前巖土力學(xué)計(jì)算中的重要數(shù)值方法之一。該程序是FLAC二維計(jì)算程序在三維空間的擴(kuò)展，用于模擬三維土體、巖體或其他地質(zhì)材料達(dá)到強(qiáng)度極限或屈服極限時(shí)發(fā)生的破壞或塑性流動(dòng)的力學(xué)特征，特別適用于分析漸進(jìn)破壞失穩(wěn)以及模擬大變形［5］。FLAC3D是采用顯式算法來(lái)獲得模型全部運(yùn)動(dòng)方程的時(shí)間步長(zhǎng)解，從而可以追蹤材料的漸進(jìn)破壞和垮落，這對(duì)研究開采的時(shí)間效應(yīng)和空間效應(yīng)是十分有效的。此外，該程序提供了多種模擬材料的本構(gòu)模型，可以輸入和改變整體或局部的模擬材料類型和參數(shù)，增強(qiáng)了程序使用的靈活性，用來(lái)提供采動(dòng)區(qū)域的垮落過(guò)程和開采中的充填過(guò)程［6］。FLAC3D擁有自己的內(nèi)嵌程序語(yǔ)言FISH，使用戶能夠自定義變量和函數(shù)，擴(kuò)大了FLAC3D的應(yīng)用及用戶自有的特色。FLAC3D還具備比較完善的后處理系統(tǒng)，用戶可以提取和查看所需要的計(jì)算信息，并且可以生成和輸出圖形、切片、動(dòng)畫等多種形式的后處理結(jié)果。

1.2 地質(zhì)采礦條件

以某礦區(qū)地質(zhì)采礦條件為模型，模擬開采區(qū)域位于某村莊下，所開采的煤層為5煤層，煤層開采的平均厚度為5.0 m，平均采深為695 m，煤層傾角為0°～2°。由于煤層的傾角很小，可以將該煤層看作水平煤層。煤層的上覆巖層以砂巖和砂質(zhì)泥巖為主，其中第四系松散層厚度約為186 m。因?yàn)槊簩娱_采的區(qū)域位于村莊下，而該村莊建筑物較為密集并且抗變形能力不強(qiáng)，容易受到采動(dòng)影響而發(fā)生破壞，所以采用固體充填開采的方法進(jìn)行開采，并且在開采之前需要進(jìn)行開采沉陷預(yù)計(jì)。各巖層的物理力學(xué)參數(shù)如表1所示。

1.3 數(shù)值模型建立

假設(shè)試采工作面長(zhǎng)度為400 m，工作面推進(jìn)距離為1 000 m。工作面采用走向長(zhǎng)壁固體充填開采的方法，充填為即采即充，充填率為100%。結(jié)合礦區(qū)的地質(zhì)采礦條件和試采工作面的尺寸等相關(guān)信息，運(yùn)用FLAC3D數(shù)值模擬軟件建立大小為3 000 m×2 400 m×715 m的三維模型，該模型共有24 800個(gè)矩形網(wǎng)格單元和261 792個(gè)單元節(jié)點(diǎn)。三維數(shù)值模型如圖1所示。

圖1 三維數(shù)值模型

2 等效采厚與模擬開挖

2.1 固體充填開采的等效采厚

固體充填開采之所以能夠減小巖層及地表的移動(dòng)變形，就是因?yàn)轫肥裙腆w充填物占據(jù)了上覆巖體的下沉空間，阻止了頂板的垮落式下沉，使得頂板巖層最終下沉較小，這就相當(dāng)于減小了煤層的開采厚度［7］。固體充填開采的等效開采厚度就是充填工作面煤層頂板從開始移動(dòng)到最終穩(wěn)定時(shí)所移動(dòng)的距離。

根據(jù)固體充填開采的工藝過(guò)程以及頂板巖層的移動(dòng)下沉機(jī)理可知，等效采厚與充填前頂板的移近量、充填的未接頂量、充填率以及固體充填體的壓縮率密切相關(guān)。固體充填開采的等效采厚的計(jì)算公式為:

式中:Me為煤層的等效開采厚度;S1為充填前頂板的移近量;S2為充填的未接頂量;S3為充填體最終的壓縮量。

其中S1、S2和S3的計(jì)算公式為:

式中:M為煤層的實(shí)際開采厚度;η為充填前頂板的下沉率(%);ρ為充填率(%);μ為充填體的壓縮率(%)。

所以固體充填開采的等效采厚的計(jì)算公式可寫為:

由于固體充填開采所引起的巖層與地表移動(dòng)變形是充填工作面的等效采厚引起的，而固體充填開采的等效采厚一般都明顯小于煤層的實(shí)際采厚，所以固體充填開采可以有效控制巖層與地表的移動(dòng)變形。在對(duì)固體充填開采所引起的地表沉陷進(jìn)行預(yù)計(jì)時(shí)，可以把它看作是開采厚度為等效采厚Me的煤層在垮落法開采條件下的情況。沉陷預(yù)計(jì)的方法可以采用概率積分法，數(shù)值模擬分析法等。

2.2 模型的模擬開挖

煤層的實(shí)際開采厚度M=5 m，充填工作面的尺寸為1 000 m×400 m。假定充填前頂板的下沉率η=5%，充填率ρ=100%，固體充填體的壓縮率μ=10%。結(jié)合固體充填開采的等效采厚計(jì)算公式，得到的等效采厚Me=0.725 m。

在模擬工作面開挖時(shí)，只需要將所開挖的區(qū)域定義為空模型(Null Model)。而在模擬工作面充填開采時(shí)，需要給之前開挖的部分區(qū)域賦予新的材料本構(gòu)模型，如摩爾－庫(kù)侖模型，雙屈服模型等，并設(shè)定對(duì)應(yīng)模型的材料參數(shù)和初始條件，然后進(jìn)行解算即可。模擬固體充填開采的FLAC3D程序?yàn)?

在以上的FLAC3D模擬開挖程序中，可以在solve這一行命令之前加上如“hist id=1 unbal”和“hist id=2 gp zdisp 1500，1200，715 ”這樣的命令，用于采樣最大不平衡力和節(jié)點(diǎn)(1500，1200，715)在 Z方向上的位移隨時(shí)間步的迭代變化，并可使用FLAC3D的plot繪圖命令顯示出來(lái)。如果需要采樣的記錄較多時(shí)，用戶可以輸入“print history”命令，顯示出什么樣的記錄號(hào)對(duì)應(yīng)什么變量或參數(shù)。在解算完成后，用戶可以通過(guò)輸入諸如“print gp disp id 257702”或“print gp disp range x 1500 y 1200 z 714 715”這些命令得到id=257702或range命令包含的相應(yīng)區(qū)域內(nèi)的節(jié)點(diǎn)位移。

在模型的模擬開挖解算完成后，可以運(yùn)用apply命令在模型上邊界的部分區(qū)域施加一定的豎向應(yīng)力，用來(lái)模擬建筑物荷載作用下的情況。比如在開挖后模型的上邊界中央?yún)^(qū)域內(nèi)施加100 kPa的一般建筑工程荷載后，就可以得到模型的上邊界增大了約1.5 cm的豎向位移?？紤]到一般村莊建筑物的荷載較小(＜100 kPa)以及巖土體的固結(jié)變形等因素，可以認(rèn)為村莊建筑物自身荷載對(duì)于煤礦開采沉陷區(qū)域內(nèi)地表的移動(dòng)變形影響很小，某些情況下甚至可以忽略不計(jì)。

由于FLAC3D具有比較強(qiáng)大的后處理功能，在模型的模擬開挖和施加建筑物荷載等解算結(jié)束之后，用戶就可以使用一系列的plot繪圖命令將所需要的圖形顯示出來(lái)。分析云圖如圖2～圖6所示。

圖2 模型的沉降位移云圖

圖3 模型的最大位移云圖

圖4 走向觀測(cè)線沉降位移云圖

圖5 傾斜觀測(cè)線沉降位移云圖

3 地表沉陷數(shù)據(jù)提取與可視化

3.1 觀測(cè)站布設(shè)及地表沉陷數(shù)據(jù)的提取

為了全面、準(zhǔn)確地反映固體充填開采對(duì)地表移動(dòng)變形的影響，需要在三維數(shù)值模型上模擬建立等間距的剖面線狀觀測(cè)站或網(wǎng)狀觀測(cè)站。由于煤層的平均采深為695 m，根據(jù)表2所示的開采深度與測(cè)點(diǎn)間距的關(guān)系，設(shè)定測(cè)點(diǎn)間距為30 m［8］。因?yàn)镕LAC3D軟件和其內(nèi)嵌的FISH程序語(yǔ)言只能查詢或輸出網(wǎng)格節(jié)點(diǎn)和單元體的空間位置等相關(guān)信息，所以在先前建立三維數(shù)值模型時(shí)，就將各個(gè)層組的網(wǎng)格劃分為長(zhǎng)寬均為30 m的單元格，以便對(duì)測(cè)點(diǎn)間距為30 m的各地表移動(dòng)觀測(cè)站進(jìn)行沉陷數(shù)據(jù)的提取。

當(dāng)模型上建立的移動(dòng)觀測(cè)站是剖面線狀觀測(cè)站時(shí)，就僅需要提取走向觀測(cè)線和傾斜觀測(cè)線上各觀測(cè)站上的沉陷數(shù)據(jù)。假定走向觀測(cè)線的Y坐標(biāo)為1 200，傾斜觀測(cè)線的X坐標(biāo)為1 500，利用FISH語(yǔ)言編制以下程序:

圖6 施加100 kPa建筑物荷載后的模型沉降位移云圖

表2 地表移動(dòng)觀測(cè)站測(cè)點(diǎn)密度表

在模型開挖的解算完成后，運(yùn)行以上程序，就可以在xyw.log文件中得到走向觀測(cè)線和傾斜觀測(cè)線上各觀測(cè)站的下沉值。

當(dāng)模型上建立的移動(dòng)觀測(cè)站是網(wǎng)狀觀測(cè)站時(shí)，則需要提取地表上開采影響區(qū)域內(nèi)所有間隔為30 m的各移動(dòng)觀測(cè)站的沉陷數(shù)據(jù)。由于FISH語(yǔ)言中的table函數(shù)只能產(chǎn)生含有X、Y這兩個(gè)字段的表，因此在提取網(wǎng)狀觀測(cè)站的沉陷數(shù)據(jù)時(shí)，就不能利用剖面線狀觀測(cè)站提取沉陷數(shù)據(jù)時(shí)采用的table、xtable、ytable等這樣的函數(shù)。針對(duì)這一問題，可以通過(guò)提取模型開挖后地表上間距為30 m的各節(jié)點(diǎn)的三維坐標(biāo)和位移來(lái)解決，對(duì)應(yīng)的FISH程序?yàn)?

運(yùn)行以上程序，就可以得到模型開挖后地表上各節(jié)點(diǎn)(移動(dòng)觀測(cè)站)的三維坐標(biāo)和位移。通過(guò)選擇xyz.log文本中的X、Y坐標(biāo)和對(duì)應(yīng)節(jié)點(diǎn)Z方向上的位移(Z－Dis)，就可以整理得到地表上各移動(dòng)觀測(cè)站的沉陷數(shù)據(jù)。

通過(guò)提取各個(gè)地表移動(dòng)觀測(cè)站的沉陷數(shù)據(jù)，得到移動(dòng)觀測(cè)站的最大沉降值為56 mm，各觀測(cè)站之間的最大傾斜值為0.102 mm/m。即使在施加了100kPa的一般建筑工程荷載后，移動(dòng)觀測(cè)站的最大沉降值僅為72 mm，各觀測(cè)站之間的最大傾斜值也僅為0.232 mm/m。結(jié)合我國(guó)一般磚木結(jié)構(gòu)建筑的臨界變形值［8］和《建筑地基基礎(chǔ)設(shè)計(jì)規(guī)范》［9］(GB50007－2011)中所規(guī)定的建筑物地基的最大沉降量(見表3)可知，各觀測(cè)站的最大沉降量和觀測(cè)站之間的最大傾斜值都遠(yuǎn)遠(yuǎn)小于規(guī)定的最大沉降量和臨界變形值。因此采用固體充填開采的方法可以有效減輕或消除煤炭開采對(duì)于地表建筑物的影響，適用于村莊等建筑物下采煤的情況。

表3 建筑物地基的最大沉降量和地表的臨界變形值表

3.2 利用MATLAB實(shí)現(xiàn)后處理與可視化操作

MATLAB是美國(guó)MathWorks公司研發(fā)的目前世界上使用最廣泛的科學(xué)計(jì)算軟件之一，它具有矩陣計(jì)算、數(shù)值分析、數(shù)據(jù)可視化以及交互式程序設(shè)計(jì)等強(qiáng)大功能［10］。在這里，本文利用MATLAB強(qiáng)大的數(shù)據(jù)可視化功能，結(jié)合提取的地表移動(dòng)觀測(cè)站的沉陷數(shù)據(jù)，運(yùn)用MATLAB的cftool曲線擬合工具箱以及plot3、mesh、contour、clabel等MATLAB繪圖函數(shù)，實(shí)現(xiàn)固體充填開采地表沉陷的可視化，繪制并輸出相關(guān)的圖形(見圖7～圖10)。

圖7 走向觀測(cè)站地表沉陷擬合曲線圖

圖8 傾斜觀測(cè)站地表沉陷擬合曲線圖

圖9 地表沉陷等高線圖

圖10 地表沉陷三維網(wǎng)格圖

如果需要得到移動(dòng)觀測(cè)站之外地表上任意點(diǎn)的下沉值，可以利用MATLAB提供的插值函數(shù)來(lái)解決。對(duì)于剖面線狀觀測(cè)站，可以利用interp1()這樣的一元插值函數(shù)來(lái)求取。而對(duì)于網(wǎng)狀觀測(cè)站的情況，則可運(yùn)用interp2()或griddata()這些二元插值命令來(lái)獲得移動(dòng)觀測(cè)站之外地表上任意點(diǎn)的下沉值。

4 結(jié)論

(1)通過(guò)對(duì)某礦區(qū)固體充填開采進(jìn)行數(shù)值模擬研究，結(jié)合固體充填開采等效采厚的相關(guān)概念，認(rèn)為固體充填開采可以有效控制地表的移動(dòng)變形，該方法對(duì)于地面的減沉效果顯著，可以應(yīng)用于村莊等建筑物下采煤的情況。

(2)運(yùn)用FLAC3D軟件建立了固體充填開采的三維數(shù)值模型，模擬了模型在固體充填開采后的移動(dòng)和變形情況，生成了固體充填開采的沉陷云圖和相關(guān)剖面圖。結(jié)合FLAC3D內(nèi)嵌的FISH程序語(yǔ)言，編制了相關(guān)程序，提取出兩種不同類型的地表移動(dòng)觀測(cè)站的沉陷數(shù)據(jù)。

(3)利用提取出的地表移動(dòng)觀測(cè)站沉陷數(shù)據(jù)，結(jié)合MATLAB強(qiáng)大的數(shù)據(jù)處理和可視化功能，對(duì)固體充填開采進(jìn)行后處理與可視化操作，生成了曲線圖、網(wǎng)格圖、等高線圖等圖形，更加形象直觀地反映出固體充填開采對(duì)于地表移動(dòng)變形的影響。

［1］ 郭增長(zhǎng)，柴華彬.煤礦開采沉陷學(xué)［M］.北京:煤炭工業(yè)出版社，2013.

［2］ 姚寶志.固體充填開采地表沉陷規(guī)律數(shù)值模擬［J］.煤礦安全，2012，43(1):171－173.

［3］ 劉瑞斌.固體充填開采地表沉陷預(yù)測(cè)及控制研究［D］.焦作:河南理工大學(xué)，2014.

［4］ 謝和平，周宏偉，王金安，等.FLAC在煤礦開采沉陷預(yù)測(cè)中的應(yīng)用及對(duì)比分析［J］.巖石力學(xué)與工程學(xué)報(bào)，1999，18(4):397－401.

［5］ 彭文斌.FLAC 3D實(shí)用教程［M］.北京:機(jī)械工業(yè)出版社，2007.

［6］ 栗帥，郭廣禮，徐斗斗，等.基于FLAC3D和SURFER的矸石充填開采沉陷數(shù)值模擬［J］.金屬礦山，2010(7):19－22.

［7］ 姜巖，徐永梅，張倬峰，等.充填開采等效厚度與充填率研究［C］//中國(guó)煤炭學(xué)會(huì).全國(guó)“三下”采煤學(xué)術(shù)會(huì)議論文集，2012.

［8］ 郭文兵，柴華彬.煤礦開采損害與保護(hù)［M］.北京:煤炭工業(yè)出版社，2008.

［9］ GB 50007－2011.建筑地基基礎(chǔ)設(shè)計(jì)規(guī)范［S］.

［10］ 陳勇，郭文兵，文運(yùn)平，等.基于MATLAB求取地表移動(dòng)預(yù)計(jì)參數(shù)的方法研究［J］.河南理工大學(xué)學(xué)報(bào):自然科學(xué)版，2009，28(6):714 －718.