張 村，宋子玉，趙毅鑫,，韓鵬華，滕 騰

( 1. 中國(guó)礦業(yè)大學(xué)( 北京 ) 共伴生能源精準(zhǔn)開(kāi)采北京市重點(diǎn)實(shí)驗(yàn)室，北京 100083；2. 中國(guó)礦業(yè)大學(xué)( 北京 ) 能源與礦業(yè)學(xué)院，北京 100083；3. 安徽理工大學(xué) 深部煤礦采動(dòng)響應(yīng)與災(zāi)害防控國(guó)家重點(diǎn)實(shí)驗(yàn)室，安徽 淮南 232001)

近年來(lái)我國(guó)煤炭去產(chǎn)能戰(zhàn)略持續(xù)推進(jìn)，煤炭主產(chǎn)區(qū)由東部地區(qū)逐漸向西部轉(zhuǎn)移，西部礦區(qū)千萬(wàn)噸級(jí)別高強(qiáng)度開(kāi)采礦井大量涌現(xiàn)。但西部礦區(qū)氣候干燥，蒸發(fā)量是降雨量的6倍左右。再加之西部礦區(qū)高強(qiáng)度的開(kāi)采，對(duì)本已脆弱的西部礦區(qū)生態(tài)水環(huán)境造成嚴(yán)重影響[1-2]。據(jù)統(tǒng)計(jì)，中國(guó)每年因煤炭開(kāi)采破壞地下水約80億t，而利用率僅25%左右，損失的礦井水資源相當(dāng)于中國(guó)每年工業(yè)和生活缺水量( 100億t )的60%[3]。因此，如何實(shí)現(xiàn)煤炭開(kāi)發(fā)與水資源保護(hù)相協(xié)調(diào)，是西部煤炭科學(xué)開(kāi)發(fā)的重大難題之一，也是礦區(qū)生態(tài)文明建設(shè)的核心內(nèi)容。據(jù)此，顧大釗院士[3]提出了煤礦地下水庫(kù)的概念和礦井水井下儲(chǔ)存利用的理念與技術(shù)。目前該技術(shù)在神東礦區(qū)應(yīng)用，已建成35座煤礦地下水庫(kù)，儲(chǔ)水量2 500萬(wàn)m3，供應(yīng)了礦區(qū)95%以上的用水，并為周邊產(chǎn)業(yè)供水，有助于西部礦區(qū)煤炭開(kāi)采水資源保護(hù)利用。

礦井地下水庫(kù)建設(shè)過(guò)程中，礦井采空區(qū)儲(chǔ)水能力和礦井地下水庫(kù)穩(wěn)定運(yùn)行是關(guān)鍵科學(xué)問(wèn)題[3]。由于垮落帶孔隙率高、滲透率大，是礦井地下水庫(kù)最主要的庫(kù)容體，垮落帶的孔隙演化規(guī)律直接影響礦井水庫(kù)的儲(chǔ)水能力[4-7]?？迓鋷в肿鳛榈V井地下水庫(kù)的主要承載體，與殘留煤柱協(xié)同承載著上覆巖層的載荷確保礦井水庫(kù)穩(wěn)定運(yùn)行[8-9]( 圖1 )，垮落帶破碎巖體遇水弱化導(dǎo)致覆巖下沉同樣影響著礦井水庫(kù)的儲(chǔ)水能力?？迓鋷?chǔ)水能力的減小，將導(dǎo)致礦井地下水庫(kù)水位大幅度提升，威脅礦井地下水庫(kù)的穩(wěn)定性。因此，垮落帶破碎巖體在礦井水庫(kù)運(yùn)行過(guò)程中孔隙演化規(guī)律是影響礦井水庫(kù)儲(chǔ)水能力的重要因素。

圖1 礦井地下水庫(kù)垮落帶+殘留煤柱協(xié)同承載結(jié)構(gòu) Fig. 1 Collaborative bearing structure of caving zone of mine underground reservoir and residual coal pillar

采空區(qū)垮落帶可以視為由破碎巖體組成的多孔介質(zhì)，隨著工作面的持續(xù)推進(jìn)，造成上覆巖層的不斷下沉壓實(shí)采空區(qū)垮落帶，垮落帶自身的應(yīng)力、密度、孔隙率、滲透率等參數(shù)均會(huì)發(fā)生變化[10]。在采空區(qū)作為礦井水庫(kù)過(guò)程中，垮落帶破碎巖體將長(zhǎng)期處于水浸受載環(huán)境中，水巖作用下破碎巖體受載內(nèi)部結(jié)構(gòu)劣化、損傷演化規(guī)律復(fù)雜多變[11-13]。在礦井地下水庫(kù)運(yùn)行過(guò)程中，破碎巖體主要發(fā)生受載破碎( 再次破碎和結(jié)構(gòu)調(diào)整 )、水浸軟化( 侵蝕溶解和遇水膨脹 )以及水滲沖蝕( 沖蝕擴(kuò)孔和小破碎巖體攜帶 )等耦合損傷( 圖2 )，進(jìn)而造成破碎巖體再次破碎、孔隙結(jié)構(gòu)調(diào)整和壓實(shí)變形，使得垮落帶孔隙空間大幅度減小，降低礦井水庫(kù)儲(chǔ)水能力。

圖2 破碎煤巖體耦合損傷與孔隙空間變化情況 Fig. 2 Coupling damage and pore space changes of broken coal and rock mass

對(duì)于垮落帶破碎煤巖體的壓實(shí)破碎特征，學(xué)者們[14-20]采用實(shí)驗(yàn)室測(cè)試和數(shù)值模擬進(jìn)行了大量的研 究，主要考慮破碎煤巖體粒徑大小、形狀、強(qiáng)度、水分、風(fēng)化程度等因素對(duì)壓實(shí)及滲流特征的影響，主要表現(xiàn)為：① 在壓實(shí)過(guò)程中破碎巖體的干密度逐漸升高；② 破碎巖體彈性模量隨著壓實(shí)應(yīng)力的增加而線(xiàn)性增加，飽和含水狀態(tài)下彈性模量的增加速度更快，但不會(huì)超過(guò)其完整狀態(tài)下的彈性模量；③ 隨著應(yīng)力的不斷升高，破碎巖體大小的分形維數(shù)將趨于一個(gè)固定值，該值與破碎巖體初始粒徑及強(qiáng)度有關(guān)。在破碎巖體初始粒徑及壓實(shí)應(yīng)力相同的情況下，分形維數(shù)隨著破碎巖體強(qiáng)度的增加逐漸減小。煤礦開(kāi)采過(guò)程中涉及水巖作用的研究主要集中在頂?shù)装逋凰?、巷道圍巖遇水穩(wěn)定性、斷層陷落柱突水以及工作面注水防片幫等方面[21-22]。為了獲得煤巖體遇水強(qiáng)度弱化的細(xì)觀(guān)機(jī)理，李桂臣[23]、鄧華鋒[24]等借助CT和離散元數(shù)值模擬進(jìn)行煤巖體遇水細(xì)觀(guān)劣化機(jī)制的研究。除此之外，很多學(xué)者[25-29]根據(jù)采動(dòng)應(yīng)力和水環(huán)境效應(yīng)研究動(dòng)靜載、加卸載、水化學(xué)侵蝕等條件下浸水煤巖體物理力學(xué)性質(zhì)的演化特征。但是到目前為止，對(duì)于礦井地下水庫(kù)循環(huán)儲(chǔ)放水過(guò)程中水體滲流對(duì)垮落帶孔隙率的影響研究相對(duì)較少。礦井儲(chǔ)放水過(guò)程中由于儲(chǔ)水高度不一致，水流速度存在差異。采空區(qū)垮落帶內(nèi)破碎巖體尺寸范圍大，從粉末顆粒到直徑數(shù)米的破碎巖體均存在。除此之外，采空區(qū)垮落帶由于壓實(shí)程度不一致，采空區(qū)中部孔隙率要小于四周，同樣影響礦井水滲流過(guò)程中對(duì)破碎巖體的攜帶作用，進(jìn)而影響采空區(qū)垮落帶的孔隙率。

據(jù)此，筆者通過(guò)理論分析，建立礦井水滲流與破碎巖體耦合作用方程，探討垮落帶孔隙率、破碎巖體尺寸以及礦井水流速對(duì)破碎巖體的攜帶作用的影響。最后將理論模型嵌入至PFC3D中，采用DEM-CFD耦合程序模擬分析了流速對(duì)破碎巖體孔隙率的影響，研究結(jié)果有助于掌握礦井水滲流過(guò)程中垮落帶孔隙率的演化規(guī)律。

1 流固耦合控制方程

1.1 流固耦合作用原理

通過(guò)附加力考慮顆粒與流體的相互作用，根據(jù)礦井地下水庫(kù)液體流動(dòng)情況可知，主要為水平流動(dòng)遷移，假設(shè)研究區(qū)域?yàn)槲挥诹黧w內(nèi)部的破碎巖體，其在礦井地下水庫(kù)的受力條件可以簡(jiǎn)化成如圖3所示的應(yīng)力環(huán)境，則有

圖3 礦井地下水庫(kù)破碎巖體受力條件 Fig. 3 Stress condition of broken rock mass in underground reservoir of mine

式中，為顆粒速度；m為顆粒質(zhì)量；為流體施加在顆粒上的總作用力( 流固耦合相互作用力 )，其由拖曳力和流體壓力梯度導(dǎo)致的力組成為作用在顆粒上的水平外力( 主要為顆粒間的水平接觸力 )之和；為重力加速度；Vm為顆粒體積；ρf為流體密度；為作用在顆粒上的垂直應(yīng)力( 主要由上覆巖層施加的力 )之和；μm為顆粒摩擦因數(shù)，進(jìn)而求出摩擦力ff；為顆粒角速度；I為慣性矩；為作用在顆粒上的力矩。

流體作用于顆粒的拖曳力是基于包含該顆粒的流體單元自身?xiàng)l件為其單獨(dú)定義的。需要說(shuō)明的是，流體-顆粒相互作用力默認(rèn)施加在顆粒形心上，并且沒(méi)有轉(zhuǎn)矩作用。因此拖曳力[30]計(jì)算公式為

這個(gè)校正項(xiàng)使這個(gè)力同時(shí)適用于高孔隙率和低孔隙率系統(tǒng)，并且支持雷諾數(shù)的大范圍取值[31-32]。單個(gè)顆粒所受拖曳力被定義為

式中，Cd為拖曳力系數(shù)，由式( 5 )計(jì)算；r為顆粒半徑；→為流體流速。

式中，Rep為顆粒雷諾數(shù)，進(jìn)而經(jīng)驗(yàn)系數(shù)χ[32]的計(jì)算公式為

顆粒雷諾數(shù)可以表示為

式中，μf為流體的動(dòng)力黏滯系數(shù)。

式中，p為流體壓力。

1.2 滲流求解方程

由于礦井水庫(kù)儲(chǔ)放水過(guò)程中水體體積較大，水體流動(dòng)流速相對(duì)較小，多孔介質(zhì)中的低雷諾數(shù)流體可由達(dá)西定律描述。需要說(shuō)明的是，達(dá)西定律對(duì)應(yīng)的雷諾數(shù)上限在1～10之間，與多孔介質(zhì)的粒徑存在一定的關(guān)系：

式中，K為滲透率。

這里假定流體的壓縮性小到可以忽略，通過(guò)隱式求解可以根據(jù)流速很快得出流體的壓力場(chǎng)。為了考慮流體流動(dòng)受顆粒的影響，滲透系數(shù)由多孔介質(zhì)模型的孔隙度計(jì)算。本文采用Kozeny-Carman方程計(jì)算滲透系數(shù)：

式中，為了計(jì)算滲透系數(shù)，同時(shí)考慮礦井水庫(kù)條件，孔隙率上限設(shè)置為0.7，孔隙率超過(guò)0.7則滲透系數(shù)取常數(shù)。

上述是針對(duì)低雷諾數(shù)的達(dá)西流滲流方程，當(dāng)雷諾數(shù)超過(guò)10時(shí)，達(dá)西流不再適用，此時(shí)多孔介質(zhì)中的滲流可用Forchheimer公式進(jìn)行描述：

式中，β為非達(dá)西滲流因子。

為了簡(jiǎn)化計(jì)算，本文主要還是以達(dá)西滲流為主進(jìn)行分析。

2 影響因素分析

2.1 垮落帶孔隙率

為了方便計(jì)算，本文采用表1中的數(shù)據(jù)進(jìn)行計(jì)算，同時(shí)假設(shè)顆粒的初始速度為0，將表1中的數(shù)據(jù)代入式( 8 )，依據(jù)本計(jì)算案例獲得的孔隙率與流體對(duì)顆粒的作用力的關(guān)系見(jiàn)式( 12 )，孔隙率由0～0.7下的流體對(duì)顆粒的作用力如圖4所示。為了實(shí)現(xiàn)正負(fù)值的統(tǒng)一，在下文分析中，統(tǒng)一將顆粒運(yùn)動(dòng)方向的力設(shè)為正。

表1 模型計(jì)算參數(shù) Table 1 Model calculation parameters

由圖4可以看出，破碎巖體所受流體壓力大小隨著孔隙率的增加急劇減小，在孔隙率逼近0時(shí)，流體壓力逼近無(wú)窮大。在孔隙率由0.1達(dá)到0.3時(shí)，流體壓力由338 N降至5.9 N，孔隙率升高至0.7時(shí)，流體壓力降至0.26 N。除此之外，由于拖曳力/流體壓力的比值K約等于1，因此，梯度力基本可以忽略不計(jì)。水中矸石顆粒的靜摩擦因數(shù)設(shè)為0.49[33]，在考 慮浮力且不考慮垮落帶上覆壓力和接觸力的情況下，顆粒加速度a的計(jì)算公式為

圖4 多孔介質(zhì)顆粒加速度及流體壓力與孔隙率的相關(guān)關(guān)系 Fig. 4 Relationship between particle acceleration and fluid force and porosity of porous media

由式( 13 )可以看出，孔隙率只影響流體壓力部分導(dǎo)致的加速度，因此，只有當(dāng)流體壓力大于摩擦力時(shí)，破碎巖體顆粒才能運(yùn)移。結(jié)合式( 13 )及圖4可以看出，在表1條件下，加速度與流體壓力變化趨勢(shì)一致，只有當(dāng)垮落帶孔隙率為0.029 6時(shí)，顆粒才能被攜帶移動(dòng)。結(jié)合采空區(qū)垮落帶實(shí)際孔隙率范圍( 隨著壓實(shí)程度的不同，在0.3～0.6之間[10])，表明在表1條件下，破碎巖體不會(huì)被流體攜帶走。在這種情況下，則需要考慮破碎巖體粒徑以及流速對(duì)流體壓力的影響。

2.2 破碎巖體尺寸

采空區(qū)垮落帶破碎巖體分布尺寸差異巨大，從粉末狀的毫米級(jí)到巖體的米級(jí)均存在。因此，有必要探索巖體尺寸對(duì)流體壓力的影響情況。在表1條件下，巖體尺寸分別與流體壓力和顆粒加速度的關(guān)系如式( 14 )和( 15)所示，對(duì)應(yīng)的曲線(xiàn)如圖5所示。

圖5 流體壓力及顆粒加速度與破碎巖體尺寸的相關(guān)關(guān)系 Fig. 5 Relationship between fluid pressure and particle acceleration and particle size

由圖5可以看出，雖然巖體尺寸的增加使得巖體所受流體壓力呈冪函數(shù)上升，但由于破碎巖體半徑的增加同樣使得破碎巖體質(zhì)量大幅度增加，導(dǎo)致總體上破碎巖體尺寸越大，破碎巖體加速度下降越快。在流體速度為0.01 m/s時(shí)，只有破碎巖體半徑為0.000 48 m的破碎巖體能被攜帶走。如果再考慮上覆巖層的應(yīng)力和破碎巖體之間的接觸力，流體速度為0.01 m/s時(shí)，基本攜帶不了破碎巖體，只能攜帶破碎巖體孔隙空間的粉塵。上述計(jì)算是在孔隙率為0.5的情況下獲得的，垮落帶中部孔隙率相對(duì)于垮落帶邊界要低得多。因此，如果相同條件下，孔隙率為0.3時(shí)，粒徑為0.001 34 m的破碎巖體將被流體攜帶，同樣很小。對(duì)于破碎巖體半徑變化過(guò)程中的拖曳力和梯度力而言，流體壓力中的拖曳力仍占絕對(duì)優(yōu)勢(shì)，拖曳力占流體壓力的比值K維持在1，壓力梯度力同樣可以忽略不計(jì)。

2.3 礦井水流速

除了垮落帶自身粒徑和孔隙率以外，流體流速同樣影響垮落帶內(nèi)能被流體攜帶的巖體尺寸，圖6是在表1條件下，垮落帶內(nèi)流體壓力及加速度與流體流速的關(guān)系( 式( 16 )～( 17 ) )。流速的增加使得流體攜帶破碎巖體的能力加大。半徑為1 m的破碎巖體，在流速達(dá)到1.925 m/s時(shí)，也能達(dá)到啟動(dòng)加速度。而假設(shè)垮落帶破碎巖體最大半徑達(dá)到10 m，則同等條件下流速需要達(dá)到6.1 m/s才能使得不考慮覆巖壓力和破碎巖體接觸力的巖體發(fā)生流動(dòng)。同樣的，由圖6可以看出，流體壓力中拖曳力占主導(dǎo)，梯度力可以忽略不計(jì)。

圖6 流體壓力及顆粒加速度與流體速度的相關(guān)關(guān)系 Fig. 6 Relationship between fluid pressure and particle acceleration and fluid velocity

3 數(shù)值模擬分析

筆者通過(guò)理論分析可以看出，礦井地下水庫(kù)在蓄放水過(guò)程中對(duì)垮落帶巖體將會(huì)產(chǎn)生攜帶作用，但在分析過(guò)程中并沒(méi)有考慮破碎巖體間的接觸力對(duì)破碎巖體運(yùn)移的影響，也無(wú)法掌握整個(gè)垮落帶孔隙率的演化情況。因此，為了更加直觀(guān)地分析流體運(yùn)移對(duì)垮落帶孔隙率的影響特征，筆者采用DEMCFD流固耦合模擬方法進(jìn)行實(shí)驗(yàn)室尺度的流固耦合模擬。

3.1 流固耦合實(shí)現(xiàn)方法

筆者采用PFC3D進(jìn)行模擬，PFC3D中內(nèi)置了計(jì)算流體動(dòng)力學(xué)( CFD )模塊，能夠在PFC3D中分析流固耦合作用問(wèn)題，但PFC3D中不包括CFD求解器，只提供了與CFD軟件連接的命令和腳本函數(shù)，需要與其他軟件進(jìn)行同步分析。

3.1.1 孔隙率計(jì)算

有2種方法可用于計(jì)算流體單元的孔隙率，采用表征顆粒的中心位置或多面體，即中心法和多面體法。本文采用計(jì)算效率更快的中心法，中心法適用于顆粒完全包含于該流體單元中，因此需要對(duì)流體網(wǎng)格和顆粒尺寸進(jìn)行限制。流體單元零孔隙度的產(chǎn)生會(huì)導(dǎo)致流體控制方程中產(chǎn)生奇點(diǎn)( 式( 3 )，拖曳力會(huì)無(wú)窮大 )，為了避免由奇點(diǎn)引起的問(wèn)題，本文限制流體單元的孔隙率不小于0.005。需要注意的是，流體單元孔隙率計(jì)算采用整個(gè)顆粒的體積，未對(duì)重疊部分進(jìn)行修正，這就造成低孔隙率多孔介質(zhì)的流體孔隙率要略大于實(shí)際孔隙率。

3.1.2 流體網(wǎng)格尺寸

PFC3D中CFD模塊支持全部由四面體或全部由六面體組成的網(wǎng)格。六面體單元表面必須是平面。為了解決流動(dòng)結(jié)構(gòu)問(wèn)題，流體網(wǎng)格需要滿(mǎn)足以下不等式：

式中，dc為流域最小寬度；Δxcfd為流體單元長(zhǎng)度。

由于上述耦合方法用于描述在一個(gè)流體單元中發(fā)生的平均耦合力。因此假定局部的孔隙度均勻分布在1個(gè)單元內(nèi)，顆粒周?chē)牧鲃?dòng)并沒(méi)有被明確地表示出來(lái)。為了得到好的結(jié)果，1個(gè)CFD單元中應(yīng)該包含若干個(gè)PFC3D顆粒，即

3.1.3 耦合實(shí)現(xiàn)

當(dāng)CFD模塊激活時(shí)，流體-顆粒相互作用力在PFC周期序列中被施加到顆粒上。PFC3D循環(huán)計(jì)算過(guò)程中，流體-顆粒相互作用力和每個(gè)流體單元的孔隙度依據(jù)給定的時(shí)間間隔重復(fù)計(jì)算。在模擬過(guò)程中可以強(qiáng)制更新上述變量，也可以開(kāi)啟和關(guān)閉耦合模式。流體力學(xué)雙向耦合是通過(guò)流體求解器和PFC3D之間進(jìn)行一系列數(shù)據(jù)交互實(shí)現(xiàn)的。每個(gè)流體單元的孔隙度取決于PFC3D。每個(gè)流體單元體積的體積力由PFC3D決定。每個(gè)流體單元中的流體速度、流體壓力、流體壓力梯度、流體密度、流體動(dòng)力黏滯系數(shù)都由流體軟件決定。

CFD模塊自動(dòng)為PFC顆粒施加流體-顆粒相互作用力。雙向耦合通過(guò)在流體流動(dòng)模型中更新孔隙率和滲透率，在PFC的CFD模塊更新流體速度場(chǎng)來(lái)實(shí)現(xiàn)。模型中設(shè)定每隔一定步數(shù)的力學(xué)周期更新穩(wěn)態(tài)流場(chǎng)。與流體求解軟件同步和數(shù)據(jù)交換可以通過(guò)FISH或Python腳本中的TCP套接字通信實(shí)現(xiàn)。

3.2 數(shù)值模擬構(gòu)建

為了提高計(jì)算效率，本文的模擬尺度為實(shí)驗(yàn)尺度，模型橫截面尺寸為0.1 m×0.1 m，長(zhǎng)度為0.2 m，模型四周為不滲透墻體。模型入口處滲流速度分別設(shè)為0.5，1和2 m/s，出口處設(shè)置壓力為0。模型中破碎巖體半徑分為2組，分別為0.001和0.003 m，占比分別是10%和90%，共31 830個(gè)球體。流體單元采用立方體模型，邊長(zhǎng)為0.02 m，符合上述尺寸要求。具體模型如圖7所示。

圖7 DEM-CFD滲流模型及孔隙率監(jiān)測(cè)點(diǎn)布置 Fig. 7 DEM-CFD seepage model and arrangement of porosity monitoring points

模擬過(guò)程中，為了考慮破碎巖體形狀對(duì)破碎巖體流動(dòng)的影響，破碎巖體接觸模型采用滾動(dòng)阻力線(xiàn)性模型( Rolling Resistance Linear Model )，模擬參數(shù)見(jiàn)表2。為了量化分析巖體不同位置處破碎巖體孔隙率的演化特征，在破碎巖體內(nèi)部布置了10組，每組9個(gè)的監(jiān)測(cè)單元，具體如圖7所示。

表2 破碎巖體接觸參數(shù)[33] Table 2 Particle contact parameters[33]

3.3 模擬結(jié)果分析

根據(jù)模擬結(jié)果分別繪制了孔隙率監(jiān)測(cè)單元第1層，中間第5層以及最后1層孔隙率的演化規(guī)律，具體如圖8所示。

圖8 不同流速模型孔隙率的演化規(guī)律 Fig. 8 Evolution of porosity in different flow velocity models

不同滲流速度條件下的破碎巖體位移、接觸力以及拖曳力如圖9所示。

根據(jù)理論分析可知，在流體速度分別為0.5，1.0和2.0 m/s，破碎巖體在不考慮上覆壓力和破碎巖體之間接觸力的情況下，1 mm半徑的破碎巖體加速度分別能夠達(dá)到272.4，1 002.7和3 767.1 m2/s；3 mm半徑的破碎巖體加速度分別能夠達(dá)到77.6，302.0和1 166.3 m2/s?？梢钥闯鲈谀M條件下，如果單個(gè)破碎巖體被流體攜帶速度會(huì)很大。但模型內(nèi)充滿(mǎn)了破碎顆粒，且模型末端設(shè)有墻體，存在反向作用力，這就導(dǎo)致顆粒并不能大規(guī)模地被流體攜帶，具體如圖9所示。由圖8孔隙率的演化規(guī)律同樣可以看出，流體速度為0.5和1 m/s時(shí)，模型的孔隙率變化幅度非常小，結(jié)合圖9顆粒的位移情況，只存在與大顆?？紫堕g的小顆粒發(fā)生位移。對(duì)于大顆粒而言，一方面同等條件下大顆粒的加速度本來(lái)就要小于小顆粒；另一方面，大顆粒由于直徑相對(duì)較大，顆粒配位數(shù)大，成為模型的骨架，接觸應(yīng)力相對(duì)較大( 圖8 )。這就導(dǎo)致在流體速度為0.5和1.0 m/s時(shí)，大顆粒的位移基本為0。只有在流體速度達(dá)到2.0 m/s時(shí)，由于拖曳力明顯增加( 圖9 )，使得外側(cè)大小顆粒均向內(nèi)擠壓，外側(cè)孔隙率大幅度增加，這由圖8( c )可以看出。

圖9 不同流速破碎巖體位移、拖曳力以及接觸力情況 Fig. 9 Particle displacement，drag force and contact force at different flow rates

對(duì)于不同滲流速度而言，主要改變的是流體對(duì)于破碎巖體的拖曳力( 梯度力可以忽略不計(jì) )，在流速由0.5 m/s上升到1.0 m/s時(shí)，拖曳力大幅度增加，但拖曳力增加的同時(shí)，破碎巖體之間的接觸應(yīng)力同樣增加，并不能打破原來(lái)的平衡結(jié)構(gòu)，只能使得游離于結(jié)構(gòu)之外的小顆粒巖體向內(nèi)部運(yùn)移，造成外部的孔隙率稍微增大，且流體速度越大，小顆粒巖體運(yùn)移越多，孔隙率增加越快。在流速上升到2.0 m/s時(shí)，大顆粒巖體結(jié)構(gòu)被沖破，向內(nèi)部擠壓，導(dǎo)致模型內(nèi)部接觸應(yīng)力大幅度提升。這使得外側(cè)的孔隙率大幅度提升，內(nèi)側(cè)孔隙率逐漸減小，進(jìn)而影響模型的整體滲透率。模擬結(jié)果揭示了破碎煤巖體水滲流過(guò)程中內(nèi)部破碎巖體遷移的過(guò)程以及孔隙率的演化情況。

4 展 望

在采空區(qū)垮落帶內(nèi)，破碎巖體不僅受到自身重力與水中浮力的影響，最主要的是受到上覆巖層的壓實(shí)應(yīng)力的影響，而在壓實(shí)應(yīng)力的作用下，破碎巖 體啟動(dòng)的拖曳力大幅度上升。因此，后期需要研究考慮壓實(shí)應(yīng)力條件下的流體滲流特征。同時(shí)，采空區(qū)垮落帶破碎巖體的粒徑分布范圍更廣，后期由實(shí)驗(yàn)室尺度擴(kuò)展至工程尺度( 擴(kuò)尺度 )的模擬是更好地服務(wù)工程應(yīng)用的基礎(chǔ)。而目前采空區(qū)擴(kuò)尺度模擬首先需要解決的是垮落帶破碎巖體的級(jí)配設(shè)置以及模擬的計(jì)算效率問(wèn)題。

此外，在模擬過(guò)程中并未考慮水流攜帶走的部分小顆粒，也沒(méi)有考慮蓄放水過(guò)程中水壓的變化。但本文的模擬方法能夠根據(jù)垮落帶真實(shí)巖體尺寸分布特征和覆巖應(yīng)力大小構(gòu)建相應(yīng)的工程尺度模型，同時(shí)根據(jù)現(xiàn)場(chǎng)實(shí)際條件下出水口的位置及允許攜帶顆粒的大小進(jìn)行模擬。在這種情況下就能模擬獲得真實(shí)條件下不同水位放水過(guò)程中孔隙率的演化情況。同時(shí)可以根據(jù)真實(shí)放水速度計(jì)算出入水口水壓的變化及覆蓋范圍，進(jìn)而可以在模擬過(guò)程中評(píng)估水位高度和流速的變化，具體工程尺度的模擬平面截圖如圖10所示。

圖10 工程尺度流固耦合模擬方案 Fig. 10 Flow-solid coupling model of caving zone under engineering-scale conditions

5 結(jié) 論

( 1 ) 針對(duì)礦井地下水庫(kù)蓄放水的工程背景，構(gòu)建了礦井水滲流的破碎巖體流固耦合模型，給出了 基于礦井水流速、垮落帶孔隙率、破碎巖體半徑等因素的破碎巖體所受流體壓力的計(jì)算方程。量化分析了垮落帶孔隙率、破碎巖體半徑以及流體流速對(duì)破碎巖體所受流體壓力的影響。

( 2 ) 流體壓力主要由拖曳力和梯度力組成，但在礦井地下水庫(kù)條件下，梯度力可以忽略不計(jì)?？迓鋷Э紫堵实脑黾?，促使流體壓力急劇減小。考慮到采空區(qū)垮落帶孔隙一般為0.3～0.6，在流體流速為0.01 m/s時(shí)無(wú)法攜帶垮落帶破碎巖體，而當(dāng)流體速度升高至1.925 m/s時(shí)，在孔隙率0.5時(shí)能攜帶半徑1 m的破碎巖體。流體壓力隨著破碎巖體粒徑的增加呈冪函數(shù)式增加，但破碎巖體的增加同時(shí)造成質(zhì)量的升高，使得破碎巖體移動(dòng)所需的流體壓力大幅度增加，進(jìn)而造成破碎巖體越小越容易被流體攜帶。

( 3 ) 基于本文構(gòu)建的破碎巖體流固耦合模型，借助PFC3D與Python求解器進(jìn)行了破碎巖體模型不同流速的數(shù)值模擬。由于破碎巖體間接觸應(yīng)力的存在，使得流體很難打破原有的平衡結(jié)構(gòu)。只有游離于破碎巖體結(jié)構(gòu)孔隙內(nèi)的小顆粒巖體才能被流體攜帶，這對(duì)整體模型的孔隙率影響較小。只有在流體流速達(dá)到一定值時(shí)，破碎巖體結(jié)構(gòu)被擠壓，模型外部孔隙率大幅度降低，內(nèi)部孔隙率則逐層出現(xiàn)遞減，進(jìn)而影響垮落帶整體的滲透率分布情況。