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

( 1. 中國礦業(yè)大學( 北京 ) 共伴生能源精準開采北京市重點實驗室，北京 100083；2. 中國礦業(yè)大學( 北京 ) 能源與礦業(yè)學院，北京 100083；3. 安徽理工大學 深部煤礦采動響應與災害防控國家重點實驗室，安徽 淮南 232001)

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

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

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

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

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

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

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

1 流固耦合控制方程

1.1 流固耦合作用原理

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

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

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

流體作用于顆粒的拖曳力是基于包含該顆粒的流體單元自身條件為其單獨定義的。需要說明的是，流體-顆粒相互作用力默認施加在顆粒形心上，并且沒有轉矩作用。因此拖曳力[30]計算公式為

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

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

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

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

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

式中，p為流體壓力。

1.2 滲流求解方程

由于礦井水庫儲放水過程中水體體積較大，水體流動流速相對較小，多孔介質中的低雷諾數(shù)流體可由達西定律描述。需要說明的是，達西定律對應的雷諾數(shù)上限在1～10之間，與多孔介質的粒徑存在一定的關系：

式中，K為滲透率。

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

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

上述是針對低雷諾數(shù)的達西流滲流方程，當雷諾數(shù)超過10時，達西流不再適用，此時多孔介質中的滲流可用Forchheimer公式進行描述：

式中，β為非達西滲流因子。

為了簡化計算，本文主要還是以達西滲流為主進行分析。

2 影響因素分析

2.1 垮落帶孔隙率

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

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

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

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

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

2.2 破碎巖體尺寸

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

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

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

2.3 礦井水流速

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

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

3 數(shù)值模擬分析

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

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

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

3.1.1 孔隙率計算

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

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

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

式中，dc為流域最小寬度；Δxcfd為流體單元長度。

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

3.1.3 耦合實現(xiàn)

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

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

3.2 數(shù)值模擬構建

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

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

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

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

3.3 模擬結果分析

根據(jù)模擬結果分別繪制了孔隙率監(jiān)測單元第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半徑的破碎巖體加速度分別能夠達到272.4，1 002.7和3 767.1 m2/s；3 mm半徑的破碎巖體加速度分別能夠達到77.6，302.0和1 166.3 m2/s。可以看出在模擬條件下，如果單個破碎巖體被流體攜帶速度會很大。但模型內(nèi)充滿了破碎顆粒，且模型末端設有墻體，存在反向作用力，這就導致顆粒并不能大規(guī)模地被流體攜帶，具體如圖9所示。由圖8孔隙率的演化規(guī)律同樣可以看出，流體速度為0.5和1 m/s時，模型的孔隙率變化幅度非常小，結合圖9顆粒的位移情況，只存在與大顆粒孔隙間的小顆粒發(fā)生位移。對于大顆粒而言，一方面同等條件下大顆粒的加速度本來就要小于小顆粒；另一方面，大顆粒由于直徑相對較大，顆粒配位數(shù)大，成為模型的骨架，接觸應力相對較大( 圖8 )。這就導致在流體速度為0.5和1.0 m/s時，大顆粒的位移基本為0。只有在流體速度達到2.0 m/s時，由于拖曳力明顯增加( 圖9 )，使得外側大小顆粒均向內(nèi)擠壓，外側孔隙率大幅度增加，這由圖8( c )可以看出。

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

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

4 展 望

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

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

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

5 結 論

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

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

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