地表塌陷區(qū)膏體回填控制顆粒流數(shù)值模擬

張友志 吳愛(ài)祥 王洪江 李 濤 孫 偉 閆其盼

(北京科技大學(xué)土木與環(huán)境工程學(xué)院，北京100083)

采礦塌陷區(qū)易誘發(fā)泥石流等地質(zhì)災(zāi)害并造成生態(tài)破壞，嚴(yán)重威脅地表建筑穩(wěn)定性和人員安全［1］。據(jù)不完全統(tǒng)計(jì)，全國(guó)因采礦引起的塌陷區(qū)有180多處，塌陷坑1 600多個(gè)，塌陷面積達(dá)1 150 km2之多。礦產(chǎn)資源開(kāi)發(fā)帶來(lái)的塌陷及其環(huán)境地質(zhì)災(zāi)害，已經(jīng)給礦區(qū)及礦業(yè)城市的可持續(xù)發(fā)展帶來(lái)嚴(yán)重的影響［2］。尾砂作為一種廉價(jià)的礦山廢料，大量用于礦山充填，既可以解決礦山充填材料的來(lái)源問(wèn)題，又解決了礦山工業(yè)廢料造成環(huán)境污染和占用耕地的問(wèn)題［3］。目前處理礦山地表塌陷區(qū)的通常做法為采用礦山廢棄物(煤矸石、粉煤灰、廢石等)直接對(duì)塌陷區(qū)進(jìn)行回填處置，但易造成回填體大規(guī)模塌陷、井下泥石流、地下水污染等次生災(zāi)害，引起礦石損失貧化。針對(duì)以上塌陷區(qū)處置方式的弊端，提出將礦山廢棄物(尾砂、廢石等)制備成混合處置體，采用柱塞泵通過(guò)管道泵送至塌陷區(qū)進(jìn)行回填，以期達(dá)到有效控制活動(dòng)塌陷區(qū)、實(shí)現(xiàn)井下安全、高效開(kāi)采的目的。

顆粒流法是一種新興的數(shù)值分析方法，在離散元理論基礎(chǔ)上由 P.A.Cundall和 O.Strack［4］提出，廣泛應(yīng)用于巖石類材料基本特性、顆粒物質(zhì)動(dòng)力響應(yīng)、巖石類介質(zhì)破裂及其發(fā)展等基礎(chǔ)性問(wèn)題的研究［5］，適用于礦區(qū)采空區(qū)中的巖體斷裂、坍塌、破碎和流動(dòng)問(wèn)題。該方法從微觀角度出發(fā)，通過(guò)采用離散單元法來(lái)模擬圓形顆粒介質(zhì)的運(yùn)動(dòng)及其相互作用，由平面內(nèi)的平動(dòng)和轉(zhuǎn)動(dòng)運(yùn)動(dòng)方程來(lái)確定每一時(shí)刻顆粒的位置和速度，能夠自動(dòng)模擬介質(zhì)基本特性隨應(yīng)力環(huán)境的變化，實(shí)現(xiàn)巖土體對(duì)歷史應(yīng)力－應(yīng)變記憶特性的模擬(屈服面變化、Kaiser效應(yīng)等)，自動(dòng)反映介質(zhì)的連續(xù)非現(xiàn)行應(yīng)力－應(yīng)變關(guān)系、屈服強(qiáng)度和此后的應(yīng)變軟化或硬化過(guò)程等。目前，顆粒流法已廣泛應(yīng)用于巖石力學(xué)行為研究［6-10］，但在地表塌陷區(qū)回填處置方面的應(yīng)用并不多見(jiàn)。

本研究以某礦地表塌陷區(qū)處置為工程背景，并基于顆粒流理論和PFC2D程序，采用Fish語(yǔ)言編制命令流，對(duì)塌陷區(qū)擬采用的回填工藝進(jìn)行了數(shù)值模擬，得出不同回填工藝下混合處置體在塌陷坑底部采礦活動(dòng)過(guò)程中的運(yùn)移變化規(guī)律，研究結(jié)果可作為地表塌陷區(qū)回填處置的參考依據(jù)。

1 工程概況

某礦細(xì)脈帶礦體自開(kāi)采以來(lái)，井下貫通地表形成的塌陷坑已有10余個(gè)，塌陷坑相互嵌套，形狀各異。如果不加處理，地表塌陷區(qū)仍有不斷擴(kuò)大的趨勢(shì)，最大回填量將達(dá)到236.9萬(wàn)m3。同時(shí)該礦年產(chǎn)礦石量達(dá)220萬(wàn)t，尾礦產(chǎn)率96%。大量的尾礦以10%的濃度被排往尾礦庫(kù)，使尾礦庫(kù)的安全與環(huán)保的壓力增大，并消耗大量的生產(chǎn)用水。另外，該礦地表碎石坑尚有35萬(wàn)m3的碎石，還有－3 mm的重介質(zhì)尾礦砂，均可以用于回填塌陷區(qū)。為了充分利用廢石、尾礦等工業(yè)固體廢棄物，擬通過(guò)向尾礦中添加碎石的方式制備膠結(jié)與非膠結(jié)2種回填處置體，對(duì)地表塌陷區(qū)進(jìn)行回填處理，達(dá)到消除塌陷區(qū)的目的，實(shí)現(xiàn)地表安全保障與環(huán)境保護(hù)的雙重目標(biāo)。

2 PFC2D數(shù)值模擬

在PFC2D的數(shù)值模擬中，顆粒和黏結(jié)體的相關(guān)參數(shù)與通常意義上的宏觀參數(shù)存在較大的區(qū)別，巖石試件的力學(xué)參數(shù)不能直接輸入到軟件中，只有通過(guò)模擬力學(xué)實(shí)驗(yàn)，進(jìn)行大量的數(shù)值模擬試驗(yàn)，根據(jù)數(shù)值試驗(yàn)的結(jié)果才能得到顆?；蛘吣z結(jié)材料的微觀參數(shù)。當(dāng)數(shù)值模擬的試件力學(xué)參數(shù)與實(shí)驗(yàn)室的宏觀參數(shù)基本相符時(shí)，即可得到顆粒的微觀參數(shù)［11-14］。只有將微觀參數(shù)輸入程序進(jìn)行計(jì)算，才能使數(shù)值模擬與實(shí)際情況相符。通過(guò)室內(nèi)實(shí)驗(yàn)對(duì)數(shù)值模型中表土層、穩(wěn)定層、擾動(dòng)帶及膠結(jié)和非膠結(jié)處置體的宏觀參數(shù)進(jìn)行了測(cè)定，各類巖層的宏觀、微觀參數(shù)見(jiàn)表1。

表1 巖層宏觀和微觀參數(shù)對(duì)照Table 1 Contrast of rock macro and microscopic parameters

以T110塌陷區(qū)為模擬對(duì)象進(jìn)行建模:模型長(zhǎng)70 m，高84 m，塌陷坑深度為40 m。針對(duì)T110塌陷區(qū)的回填情況，擬采取完全膠結(jié)回填、完全非膠結(jié)回填、10 m厚度分層回填、20 m厚度分層回填及5 m、15 m厚度交替回填這5個(gè)方案對(duì)回填情況進(jìn)行模擬。原始塌陷坑如圖1所示，以完全膠結(jié)回填為例，回填效果模型如圖2所示。

圖1 原始塌陷坑Fig.1 Original collapse pit

圖2 完全膠結(jié)回填效果模型Fig.2 Effect model of fully cemented backfill

不同回填方式的初始模型建成后，開(kāi)始對(duì)塌陷坑底部的采礦活動(dòng)進(jìn)行模擬。按照該礦的年采礦量對(duì)塌陷坑底部礦石進(jìn)行采出模擬，同時(shí)在回填體表層設(shè)置一系列監(jiān)測(cè)點(diǎn)對(duì)回填體的沉陷位移進(jìn)行監(jiān)測(cè)，在塌陷坑底部設(shè)置測(cè)量圓對(duì)回填體與塌陷坑的接觸狀況及應(yīng)力分布狀況進(jìn)行測(cè)量。監(jiān)測(cè)點(diǎn)及測(cè)量圓布置情況見(jiàn)圖3。

圖3 監(jiān)測(cè)點(diǎn)及測(cè)量圓布置Fig.3 Layout schematic of monitoring point and measuring round

3 數(shù)值模擬結(jié)果分析

以完全膠結(jié)回填、完全非膠結(jié)回填和5 m、15 m厚度交替回填為例，塌陷坑底部采礦活動(dòng)結(jié)束后，不同回填方式下的處置體及塌陷坑運(yùn)移情況分別如圖4～圖6所示。

圖4 完全膠結(jié)回填方案運(yùn)移情況Fig.4 Migration situation of fully cemented backfill scheme

圖5 完全非膠結(jié)回填方案運(yùn)移情況Fig.5 Migration situation of completely non-cemented backfill scheme

圖6 5 m、15 m交替回填方案運(yùn)移情況Fig.6 Migration situation of 5 m and 15 m alternative backfill scheme

3.1 監(jiān)測(cè)點(diǎn)沉陷位移測(cè)量

塌陷區(qū)內(nèi)處置體在回填過(guò)程中會(huì)出現(xiàn)不同程度下陷的情況，通過(guò)對(duì)處置體表層顆粒在重力方向的位移變化情況進(jìn)行監(jiān)測(cè)，可在很大程度上對(duì)處置體整體下陷情況做出評(píng)判，從而對(duì)回填方案的優(yōu)劣情況做出評(píng)判。各方案在圖3中5個(gè)監(jiān)測(cè)點(diǎn)的沉陷位移曲線見(jiàn)圖7。

從圖7中可以得出:

(1)塌陷區(qū)內(nèi)回填處置體在底部開(kāi)挖過(guò)程中出現(xiàn)不均勻下降的現(xiàn)象，其中處置體與塌陷區(qū)邊幫接觸處的沉陷位移量較大，這是由于回填處置體和塌陷區(qū)邊幫接觸處不易形成穩(wěn)定的膠結(jié)層造成的。另外處置體中心有較大沉陷位移量，這與工程實(shí)際也是相吻合的。

(2)完全膠結(jié)回填方案中處置體沉陷位移量最小，中心沉陷位移量為1.1 m，完全非膠結(jié)回填方案的沉陷位移量最大，中心沉陷位移量為2 m。分層回填方案的中心沉陷位移量在1.3～1.35 m，三者相差不大。從經(jīng)濟(jì)角度出發(fā)，優(yōu)先推薦5、15 m分層交替回填方式。

3.2 回填體在底部開(kāi)挖后的應(yīng)力分布

圖8為回填模型在塌陷坑底部礦石開(kāi)挖后的接觸力分布情況。顆粒間接觸力的大小與線條粗細(xì)成正比例關(guān)系。從圖8中可以看出，模型中穩(wěn)定層兩側(cè)的顆粒間接觸力均以豎向壓力作用為主，說(shuō)明在模擬計(jì)算過(guò)程中，擾動(dòng)帶內(nèi)顆粒間存在一定的空隙空間，顆粒間的接觸力較小，而穩(wěn)定層內(nèi)顆粒間接觸力則相對(duì)較大。此外，在圖8中的淺色標(biāo)記區(qū)域內(nèi)，顆粒間的接觸力以橫向作用力為主，形成一個(gè)橫向應(yīng)力拱，并與穩(wěn)定層顆粒共同構(gòu)成了一個(gè)連續(xù)的力學(xué)拱結(jié)構(gòu)，這對(duì)阻止塌陷區(qū)持續(xù)下陷、維持塌陷區(qū)的穩(wěn)定性起到了重要的穩(wěn)固作用。因而在接觸力分布圖中，塌陷區(qū)下方橫向力越明顯，說(shuō)明顆粒間橫向壓力作用越大，形成的橫向應(yīng)力拱也越穩(wěn)固，對(duì)維持塌陷區(qū)的穩(wěn)定性也就越有利。

圖8 回填模型在底部礦石開(kāi)挖后的接觸力分布情況Fig.8 Distribution of contact force after ore being excavated at bottom

3.3 回填體與塌陷坑離層分析

隨著塌陷坑內(nèi)處置體的不斷下移，處置體與塌陷坑底部巖層不一定能夠良好接觸，從而在接觸處可能會(huì)產(chǎn)生一定的離層空間。接觸處離層現(xiàn)象形成的空洞會(huì)對(duì)處置體的均勻下移起到阻礙作用，從而致使塌陷坑內(nèi)處置體不能連續(xù)均衡穩(wěn)定下移，對(duì)塌陷坑的穩(wěn)定性造成影響。離層空間的大小是通過(guò)離層因子來(lái)判定的，在PFC2D中離層因子可通過(guò)對(duì)測(cè)量圓內(nèi)的孔隙率進(jìn)行測(cè)量獲得，從而對(duì)不同方案的離層程度進(jìn)行分析與評(píng)判:離層因子值較小，則說(shuō)明測(cè)量圓內(nèi)部顆粒分布比較均勻，即處置體顆粒與塌陷坑底部接觸較為良好;而離層因子值較大則說(shuō)明塌陷坑底部與處置體接觸不佳。5個(gè)回填方案離層因子的測(cè)量數(shù)據(jù)見(jiàn)表2。

表2 各回填方案離層因子Table 2 Separation factor of each backfill scheme

從表2中可以得出:

(1)完全膠結(jié)回填方案的離層因子最大，離層因子值為0.18，完全非膠結(jié)回填方案的離層因子最小，值為0.12。

(2)方案3、4、5 的離層因子分別為0.15、0.16 和0.148，三者相差不大。

4 結(jié)論

(1)塌陷坑底部礦石的開(kāi)采會(huì)對(duì)塌陷處回填體造成影響，使之產(chǎn)生位移，完全膠結(jié)回填方案中回填處置體受到的影響最小，沉陷位移也最小，分層回填方案效果次之，完全非膠結(jié)回填方案位移最大，控制處置體位移的能力也最差。

(2)膠結(jié)作用下，回填處置體與周圍巖體形成新的應(yīng)力分布體系，所形成的橫向連續(xù)應(yīng)力拱結(jié)構(gòu)可以阻止塌陷區(qū)下陷，從而對(duì)地表塌陷區(qū)穩(wěn)定性起到較好的控制作用。

(3)回填處置體與塌陷坑底部巖層之間形成的離層空間對(duì)處置體的整體穩(wěn)定均衡下移起到阻礙作用，膠結(jié)回填層厚度不小于5 m、非膠結(jié)回填層厚度5～15 m時(shí)，塌陷坑的穩(wěn)定性最佳。

［1］ Olivier D，Harlalka A.Numerical study of the soil-structure interaction within mining subsidence areas［J］.Computers and Geotechnics，2010，37:802-816.

［2］ Hu R L，Yue Z Q，Wang L C，et al.Review on current status and challenging issues of land subsidence in China［J］.Engineering Geology，2004，76:65-77.

［3］ 劉同有.充填采礦技術(shù)與應(yīng)用［M］.北京:冶金工業(yè)出版社，2001.Liu Tongyou.Filling Mining Technology and Application［M］.Beijing:Metallurgical Industry Press，2001.

［4］ Cundall P A，Strack O.Particle Flow Code in 2D［M］.Minnesota:Itasca Consulting Group，Inc.，1999.

［5］ 吳順川，周 喻，高 斌.卸載巖爆試驗(yàn)及PFC3D數(shù)值模擬研究［J］.巖石力學(xué)與工程學(xué)報(bào)，2010，29(S2):4082-4088.Wu Shunchuan，Zhou Yu，Gao Bin.Study of unloading tests of rock burst and PFC3Dnumerical simulation［J］.Chinese Journal of Rock Mechanics and Engineering，2010，29(S2):4082-4088.

［6］ 王培濤，楊天鴻，柳小波.邊孔角對(duì)無(wú)底柱分段崩落法放礦影響的顆粒流數(shù)值模擬研究［J］.金屬礦山，2010(3):12-16.Wang Peitao，Yang Tianhong，Liu Xiaobo.Particle flow numeral simulation investigation on influence of lateral opening angle on ore drawing with sublevel pillarless caving［J］.Metal Mine，2010(3):12-16.

［7］ 劉增輝，高 謙，李 欣，等.回填廢石層下放礦損貧預(yù)測(cè)及顆粒流數(shù)值分析［J］.金屬礦山，2011(2):5-8.Liu Zenghui，Gao Qian，Li Xin，et al.Ore dilution prediction of ore drawing below backfilling of waste rock layer and particle flow numerical analysis［J］.Metal Mine，2011(2):5-8.

［8］ Potyondy D O，Cundall P A.A bonded-particle model for rock［J］.International Journal of Rock Mechanics ＆ Mining Sciences，2004，41(8):1329-1364.

［9］ Potyondy D O，Cundall P A.Simulating stress corrosion with a bonded-particle model for rock［J］.International Journal of Rock Mechanics ＆ Mining Sciences，2007，44(5):677-691.

［10］ Wang C，Tannant D D，Lilly P A.Numerical analysis of the stability of heavily jointed rock slopes using PFC2D［J］.International Journal of Rock Mechanics ＆ Mining Sciences，2003，40:415-424.

［11］ Jeoungseok Y.Application of experimental design and optimization to PFC model calibration in uniaxial compression simulation［J］.International Journal of Rock Mechanics ＆ Mining Sciences，2007，44:871-889.

［12］ Okan S，Nuri A A.Numerical simulation of rock cutting using the discrete element method［J］.International Journal of Rock Mechanics ＆ Mining Sciences，2011，48:434-442.

［13］ Tawadrous A S，DeGagn D.Prediction of uniaxial compression PFC3Dmodel micro-properties using artificial neural networks［J］.International Journal for Numerical and Analytical Methods in Geomechanics，2009，33:1953-1962.

［14］ Kazerani T，Zhao J.Micromechanical parameters in bonded particle method for modeling of brittle material failure［J］.International Journal for Numerical and Analytical Methods in Geomechanics，2010，34:1877-1895.