崩落法開采上覆崩落體降雨入滲及突水模擬研究

王 平 曾文旭 程愛平 許夢(mèng)國(guó) 張玉山

（1.武漢科技大學(xué)資源與環(huán)境工程學(xué)院，湖北武漢430081；2.冶金礦產(chǎn)資源高效利用與造塊湖北省重點(diǎn)實(shí)驗(yàn)室，湖北武漢430081）

強(qiáng)降雨期間塌陷坑內(nèi)匯集的雨水通過上覆崩落體和圍巖裂隙通道向采場(chǎng)入滲，滲流水導(dǎo)致碎石間的摩擦系數(shù)下降從而使其接觸力降低，增加了突水通道形成的可能性；雨水或碎屑物質(zhì)還可能在采場(chǎng)上方的較大縫隙或空區(qū)內(nèi)滯留，使局部壓力增大，同時(shí)在重力的共同作用下形成突水通道的可能性極高。例如程潮鐵礦使用無(wú)底柱分段崩落法開采多年來，在地表形成約93.4 hm2的塌陷區(qū)，在強(qiáng)降雨時(shí)共引發(fā)過4次大規(guī)模井下突水突泥事故，水裹挾大量泥土和礦巖碎屑從采場(chǎng)進(jìn)路涌出，不僅使大量巷道和生產(chǎn)設(shè)備毀壞，還造成人員傷亡，給礦山造成巨大損失。

許多學(xué)者對(duì)井下突水突泥進(jìn)行了相關(guān)研究，水從源頭到達(dá)井下引發(fā)突水突泥是一個(gè)復(fù)雜的過程。文獻(xiàn)［1—2］研究了具體礦山突水災(zāi)害的形成機(jī)理，認(rèn)為井下突水的形成需要具備流動(dòng)通道和充足的水源這2個(gè)必要條件；文獻(xiàn)［3—4］通過數(shù)值模擬和相似實(shí)驗(yàn)研究了由上覆崩落體導(dǎo)致采場(chǎng)突水突泥時(shí)關(guān)鍵部位的突變過程；文獻(xiàn)［5］進(jìn)行了地表黃土層降雨入滲機(jī)理研究；文獻(xiàn)［6］建立了崩落法礦山井下采場(chǎng)突水突泥形成影響因素的指標(biāo)體系。

綜上所述，對(duì)于降水在上覆崩落體中滲流的研究還較少，本研究結(jié)合松散礦巖的滲流特性及無(wú)底柱分段崩落法礦山采場(chǎng)突水機(jī)理，構(gòu)建數(shù)值計(jì)算模型，分析不同工況下上覆崩落體中滲流速度和壓力分布狀況，研究降水在上覆崩落體中的入滲規(guī)律，對(duì)于預(yù)測(cè)和預(yù)防采場(chǎng)突水突泥具有重要意義。

1 入滲突水機(jī)制

上覆巖層的崩落方式有自然崩落或通過人工爆破崩落，崩落過程中形成的巖塊塊度較大，隨著井下后續(xù)的開采和出礦，崩落范圍越來越大，崩落范圍內(nèi)的巖塊也逐漸下降，致使上覆崩落體得以充分松散，可以將其看成是一種具有較強(qiáng)透水性的多孔介質(zhì)。流體在上覆崩落體內(nèi)滲流，具有流速快、雷諾數(shù)高等特點(diǎn)，流體在其中不僅要受到滲透壓力和孔隙水壓力的作用，還需要考慮黏性流體的剪切應(yīng)力作用，這是以黏性力為主導(dǎo)、適合滲流速度較慢的Darcy定律所無(wú)法描述的［7］。Darcy方程描述的是牛頓流體在某種多孔介質(zhì)中滲流，滲流速度與水勢(shì)梯度正相關(guān)的滲流基本定律，該方法適應(yīng)于地下水在含水層中的滲流計(jì)算；Navier-Stokes方程是不可壓縮黏性流體基于牛頓第二定律的表達(dá)式，它是在動(dòng)能、勢(shì)能和靜壓能平衡同時(shí)考慮流體在重力、壓力和黏性阻力作用下的運(yùn)動(dòng)規(guī)律，該方程適用于管道內(nèi)流體的計(jì)算。Darcy方程和Navier-Stokes方程顯然不適合上覆崩落體內(nèi)的滲流計(jì)算。H C Brinkman在Darcy方程的基礎(chǔ)上考慮Navier-Stokes方程中的黏性流體剪切應(yīng)力項(xiàng)，提出了Brinkman方程，此方程是基于牛頓第二定律，描述了在剪切應(yīng)力和滲透壓力作用下快速滲流體的運(yùn)動(dòng)規(guī)律，對(duì)流體在礦巖松散體內(nèi)的滲流提供了較好的表達(dá)［8］。本研究將采用此方程來描述上覆崩落體內(nèi)的滲流。

Darcy定律是對(duì)飽和砂土中滲流流速和壓力之間線性關(guān)系的描述，通過飽和砂的滲流實(shí)驗(yàn)得出：

式中，V為流體流速，m/s；μ為動(dòng)黏系數(shù)，Pa·s；K為滲透率，m2；p為流體壓力，Pa；Z為位置高度，m；ρ為流體密度，kg/m3；?p為廣義壓力，Pa。

根據(jù)牛頓第二定律，流體微元所受重力與壓力之和等于其加速度和質(zhì)量之乘積，這是對(duì)歐拉方程的描述，Navier-Stokes方程是在歐拉方程的基礎(chǔ)上考慮了黏性力，式（2）為Navier-Stokes方程：

在Darcy方程的基礎(chǔ)上考慮Navier-Stokes方程中的黏性剪切應(yīng)力項(xiàng)，其中壓力?p為壓力項(xiàng)和重力項(xiàng)合并后的廣義壓力，得到Brinkman方程：

針對(duì)巖體破碎帶，文獻(xiàn)［9］中對(duì)斷層活化導(dǎo)水進(jìn)行模擬后發(fā)現(xiàn)發(fā)生突水時(shí)斷層導(dǎo)水性明顯增強(qiáng)，文獻(xiàn)［10］引入滲透率突跳系數(shù)概念，認(rèn)為滲流模型中單元破壞后滲透系數(shù)在原來的基礎(chǔ)上會(huì)增大若干倍；文獻(xiàn)［11］中對(duì)峰后破碎巖石采用立方定律和strain partition技術(shù)建立滲流—體應(yīng)變關(guān)系方程和數(shù)值模型，該模型能夠定量描述峰后破碎巖石滲透性的急劇增大，即破碎巖石散體的突水也符合滲透率突跳系數(shù)的概念。這里將引用滲透率突跳系數(shù)概念，用滲透率急劇增大來模擬大突水通道的形成。

2 入滲突水物理模型

根據(jù)無(wú)底柱分段崩落法采場(chǎng)結(jié)構(gòu)及上覆崩落體形態(tài)，簡(jiǎn)化后建立如圖1所示的降雨入滲及突水模型，程潮鐵礦塌陷坑內(nèi)黃土層長(zhǎng)時(shí)間受塌陷擾動(dòng)和雨水沖刷，不僅嚴(yán)重流失而且存在許多裂隙通道，所以該模型忽略黃土層的存在；模型主要由上覆崩落體和崩落法采場(chǎng)進(jìn)路巷道構(gòu)成，其中分段高度取8 m，采場(chǎng)進(jìn)路巷道高度取3 m；上覆崩落體整體下移時(shí)，不同粒徑散體的移動(dòng)速度不一致，細(xì)小顆粒相對(duì)粗大顆粒移動(dòng)速度較快，發(fā)生上覆巖層的自然分級(jí)現(xiàn)象［12］，其中在礦體上方并與礦體接觸的覆蓋巖層一般由人工爆破放頂形成，它的主要作用是崩礦時(shí)形成擠壓爆破條件和對(duì)冒落圍巖起到緩沖，覆蓋巖層的形成方式和整個(gè)上覆崩落體的分級(jí)現(xiàn)象使這一區(qū)域相對(duì)上覆圍巖自然崩落區(qū)域具有散體塊度小和較密實(shí)的特點(diǎn)，同時(shí)滲透率也相對(duì)較小，模型中將上覆崩落體劃分為覆蓋巖層區(qū)域和自然崩落區(qū)域，覆蓋巖層的厚度取約等于2個(gè)分段高度。

圖1中水從上邊界入口開始向上覆崩落體內(nèi)入滲，在采場(chǎng)流出進(jìn)入回采巷道或流向下一分層。模擬中不考慮兩側(cè)其他水源對(duì)研究區(qū)域的影響，邊界為滑移屬性的隔水邊界，與流體間不存在黏滯效應(yīng)。

3 入滲及突水模擬

COMSOL Multiphysics系統(tǒng)也稱用FEMLAB系統(tǒng)，是以有限元為基礎(chǔ)，通過求解偏微分方程或偏微分方程組來實(shí)現(xiàn)對(duì)科學(xué)和工程領(lǐng)域的各種物理過程的模擬。構(gòu)建計(jì)算模型，運(yùn)用COMSOL Multiphysics系統(tǒng)中自由和多孔介質(zhì)流動(dòng)模塊，對(duì)場(chǎng)變量不隨時(shí)間變化的穩(wěn)態(tài)滲流及突水進(jìn)行模擬。

3.1 計(jì)算模型及模擬過程

根據(jù)圖1所示的物理模型，在COMSOL數(shù)值軟件中建立上覆崩落體滲流及突水計(jì)算模型（圖2）。模型長(zhǎng)×高約等于80 m×100 m，頂部邊界為進(jìn)水邊界。模擬時(shí)入口根據(jù)不同工況可以設(shè)置為壓力邊界條件和流量邊界條件，雨水在靜水壓力和重力的驅(qū)動(dòng)下從地表塌陷坑開始依據(jù)Brinkman滲流方程向上覆崩落體內(nèi)入滲，出口為自由邊界，其壓力值設(shè)為0 Pa，位于模型下方，分別對(duì)應(yīng)回采巷道和下個(gè)分層。

程潮鐵礦礦區(qū)某天降雨強(qiáng)度達(dá)到63.7 mm，塌陷區(qū)降雨入滲當(dāng)量速度（單位面積的入滲水量）為1.63×10-6m/s，根據(jù)參數(shù)反演法得到程潮鐵礦上覆崩落體的平均滲透系數(shù)為8.03×10-4m/s（即平均滲透率為1.07×10-10m2），自然崩落區(qū)和覆蓋巖層區(qū)滲透率根據(jù)不同孔隙率來選取，具體參數(shù)見表1所示［13］。假設(shè)流體的密度和黏度都是常數(shù)；水流從采場(chǎng)出口流出后在回采進(jìn)路內(nèi)自由流動(dòng)；模型中流體在各區(qū)域邊界兩側(cè)的壓力和流速相等，即流速和壓力保持連續(xù)性。

首先模擬發(fā)生突水前降雨入滲階段，入口設(shè)置為恒定流量邊界，保持流速當(dāng)量為1.63×10-6m/s，在此模型入口相當(dāng)于約3.74 kPa的靜水壓力，研究突水前上覆崩落體內(nèi)滲流流速和壓力分布情況；第二步假設(shè)覆蓋巖層之上的自然崩落巖體中存在不規(guī)則空隙，黏性流體容易在其中滯留使壓力增大，在出礦和重力作用下從采場(chǎng)爆堆涌出，形成突水大通道，所以把覆蓋巖層上界面到出口這部分區(qū)域定為突水通道的形成區(qū)域，引用滲透率突跳系數(shù)概念，將此區(qū)域的滲透率提高10倍，來模擬突水通道形成后階段，入口邊界條件分別保持恒定流量和恒定水壓，分別研究采場(chǎng)發(fā)生突水后入滲流量不變時(shí)和保持初始水壓不變時(shí)滲流場(chǎng)流速和壓力的分布狀況。

3.2 結(jié)果分析

模擬結(jié)果如圖3～圖6所示，圖3和圖4分別是工況一中流體在上覆崩落體內(nèi)的流速和壓力分布云圖，圖5和圖6分別是地表塌陷坑到井下采場(chǎng)（沿圖3中折線A-B-C）的流速和壓力分布曲線。其中工況一是突水之前降雨入滲模擬，入口保持流量恒定；工況二表示入口邊界條件依舊為流量恒定時(shí)的突水工況；工況三表示維持入口邊界條件為恒定水壓時(shí)的突水工況。

從圖5中可以看出，工況一中滲流場(chǎng)流速在接近出口時(shí)發(fā)生大幅度變化，工況二、三中在突水通道區(qū)域開始發(fā)生大幅度變化，在出口方向上呈指數(shù)型增長(zhǎng)，在此之前速度保持小幅度增長(zhǎng)；工況一中流速突變之前平均值約為2×10-6m/s，出口處平均流速約1.5×10-5m/s，這是由于滲流面積總體縮小造成的；滲流壓力表現(xiàn)為自然崩落區(qū)域卸壓較慢，而靠近出口的卸壓速度快速提高，這是由于滲流場(chǎng)下半?yún)^(qū)域的滲透率較上半部分低，對(duì)上半?yún)^(qū)域流體的卸壓起到阻礙作用，而靠近出口時(shí)由于流體速度極大提高，使壓力得以迅速釋放。

工況二中采場(chǎng)上方突水通道打開后，保持入口邊界補(bǔ)給水量不變，入口邊界水壓從3.74 kPa降到了2.3 kPa；滲流速度和工況一相比，突水通道打開使流速增速加大，在采場(chǎng)出口處平均流速達(dá)到2.5×10-5m/s，這表明采場(chǎng)上方突水通道打開之后會(huì)起到大幅度的卸壓作用，流速也會(huì)極大提高，滲流場(chǎng)中流向突水采場(chǎng)的水量會(huì)增大。

工況三滲流場(chǎng)入口為恒定壓力的邊界條件，突水通道打開之后，對(duì)自然崩落區(qū)域的滲流失去阻礙作用，而且此時(shí)流體已獲得較快的速度，使上半?yún)^(qū)域流體壓力得以迅速釋放（如圖6）；滲流速度變化和工況二相似，但每個(gè)階段都明顯增大，出口處流速最高達(dá)到6×10-5m/s，約是工況一的4倍；這說明當(dāng)突水水源保持恒定高水壓力時(shí)，一旦發(fā)生突水災(zāi)害后所產(chǎn)生的破壞力也越大。

4 礦山記錄數(shù)據(jù)對(duì)比

程潮鐵礦2001年進(jìn)行強(qiáng)含水帶疏干作業(yè)之后，地下水對(duì)礦山開采的影響很微弱，礦區(qū)在2005年之前平均年降雨量1 218.73 mm，4—7月雨季時(shí)日降雨量30～216.6 mm，礦區(qū)地形特點(diǎn)有利于周邊地表徑流向陷落區(qū)匯集，大氣降雨和地表徑流沿陷落區(qū)或移動(dòng)區(qū)進(jìn)入地下是開采期間最主要的充水因素。2005年時(shí)統(tǒng)計(jì)近幾年井下雨季平均涌水量，其中西區(qū)平均涌水量4 000 m3/d，當(dāng)時(shí)西區(qū)在-307.5 m和-325 m水平同時(shí)采礦，2個(gè)分層共有約120個(gè)回采巷道，計(jì)算得單個(gè)回采巷道日平均涌水流量約33.3 m3/d；工況一為日降雨強(qiáng)度達(dá)到63.7 mm時(shí)上覆崩落體內(nèi)雨水入滲模擬，采場(chǎng)滲流出口處平均流速約1.5×10-5m/s，掌子面爆堆豎直投影面積約13.35 m2，計(jì)算得單個(gè)回采巷道涌水量為17.3 m3/d。

模擬結(jié)果與實(shí)際數(shù)據(jù)有一定偏差，主要是源于模型尺寸、邊界條件、參數(shù)分布等與現(xiàn)場(chǎng)存在差別，另外也沒有考慮崩落塌陷區(qū)之外的圍巖裂隙滲流。雖然有不可避免的誤差存在，但模擬結(jié)果與實(shí)際記錄數(shù)據(jù)在同一數(shù)量等級(jí)，模擬計(jì)算得到的滲流流速和壓力變化規(guī)律與實(shí)際情況符合，今后可以考慮更多條件，建立更完善的立體模型進(jìn)行模擬。

5 結(jié)論

（1）采用Brinkman方程比較適合上覆崩落體內(nèi)滲流的表征，根據(jù)上覆礦巖散體滲透率上大下小的特點(diǎn)，水向下滲流過程中，卸壓速度總體上呈先慢后快的趨勢(shì)，崩落體上部分滲流流體壓力較高，是促發(fā)突水通道形成的一個(gè)重要因素；

（2）地表塌陷坑匯水深度越大，井下采場(chǎng)一旦發(fā)生突水所造成的損失程度越大，甚至形成井下泥石流，在降雨天氣做好塌陷坑內(nèi)的防排水工作能有效減少井下突水事故發(fā)生幾率或降低事故危害。