徐 超，曹明月，李小芳，孫浩石，閆志銘

（1.中國(guó)礦業(yè)大學(xué)（北京）共伴生能源精準(zhǔn)開(kāi)采北京市重點(diǎn)實(shí)驗(yàn)室，北京100083；2.中國(guó)礦業(yè)大學(xué)（北京）應(yīng)急管理與安全工程學(xué)院，北京100083；3.山東科技大學(xué) 礦山災(zāi)害預(yù)防控制重點(diǎn)實(shí)驗(yàn)室，山東 青島266590；4.陽(yáng)泉煤業(yè)（集團(tuán)）有限責(zé)任公司，山西 陽(yáng)泉045000）

煤與瓦斯突出等瓦斯事故嚴(yán)重威脅煤礦安全生產(chǎn)，而且礦井瓦斯事故多發(fā)生在綜采工作面[1-2]。因此基于FLUENT模擬軟件研究采空區(qū)內(nèi)瓦斯運(yùn)移及富集規(guī)律對(duì)預(yù)防礦井瓦斯事故有指導(dǎo)性作用。采空區(qū)是由垮落的巖石和遺煤構(gòu)成的多孔介質(zhì)，多孔介質(zhì)的空隙率和滲透率分布對(duì)流體流動(dòng)有較大影響。李樹剛等通過(guò)研究垮落巖石的碎脹特性，提出了基于碎脹系數(shù)的采空區(qū)空隙率計(jì)算公式，并基于空隙率模型研究了采空區(qū)瓦斯?jié)B流規(guī)律[3]；高建良等研究了在采空區(qū)滲透率均勻分布、分段均勻分布和非均勻連續(xù)分布下的采空區(qū)流體運(yùn)移規(guī)律，其中非均勻連續(xù)分布最符合現(xiàn)實(shí)情況[4]；陳鵬等研究采空區(qū)及其上覆巖層的巖石碎脹性質(zhì)，基于“O”型圈理論提出了更符合實(shí)際的采空區(qū)空隙率三維分布模型[5]。在參考前人研究結(jié)果的基礎(chǔ)上，以平舒煤礦15111綜采工作面為研究對(duì)象，基于“O”型圈和砌體梁理論建立采空區(qū)三維空隙率和滲透率模型，運(yùn)用FLUENT軟件，模擬在有無(wú)重力2種情況下的采空區(qū)瓦斯流動(dòng)狀態(tài)，此外還對(duì)的Blake-Kozeny公式涉及到的多孔介質(zhì)平均粒子直徑的取值進(jìn)行了模擬，得到了在不同粒徑下采空區(qū)瓦斯運(yùn)移和富集規(guī)律，并對(duì)比模擬結(jié)果和實(shí)測(cè)數(shù)據(jù)，分析模擬結(jié)果的準(zhǔn)確性。

1 采空區(qū)三維流場(chǎng)數(shù)學(xué)模型

煤層開(kāi)采后采場(chǎng)中形成矩形采動(dòng)空間，受二次應(yīng)力影響，圍巖發(fā)生垮落、斷裂和變形，形成了采場(chǎng)上覆巖體結(jié)構(gòu)的“砌體梁”模型[6]?；凇捌鲶w梁”理論和“O”型圈理論，提出如下模型：取進(jìn)風(fēng)巷側(cè)工作面與采空區(qū)的交點(diǎn)為坐標(biāo)原點(diǎn)O，x軸方向?yàn)檠夭煽諈^(qū)走向方向，y軸方向?yàn)檠夭煽諈^(qū)傾向方向，z軸方向?yàn)檠夭煽諈^(qū)豎直方向。

1.1 采空區(qū)三維空隙率模型

沿采空區(qū)走向方向，破碎煤體碎脹系數(shù)引起多孔介質(zhì)空隙率變化范圍較大。隨著距工作面越遠(yuǎn)，空隙率呈負(fù)指數(shù)關(guān)系遞減[7]。因此基于以上關(guān)系，建立采空區(qū)空隙率沿x軸方向的分布函數(shù)。

由煤巖體碎脹系數(shù)定義可知采空區(qū)內(nèi)破碎巖體的碎脹系數(shù)沿x軸的分布為：

式中：kp為破碎巖石的碎脹系數(shù)，無(wú)量綱；∑h為直接頂厚度，m；m為煤層開(kāi)采厚度，m；W為采空區(qū)上覆巖層下沉量，m。

式中：kp1為直接頂破碎后巖塊碎脹系數(shù)，無(wú)量綱；x為采空區(qū)內(nèi)某一點(diǎn)距工作面的距離，m；l為基本頂破斷后巖體長(zhǎng)度，m。

根據(jù)破碎煤巖體空隙率與碎脹系數(shù)之間的關(guān)系，空隙率沿x軸方向的分布函數(shù)如下：

式中：n（x）為采空區(qū)底板處y=0時(shí)沿走方向方向空隙率，無(wú)量綱。

根據(jù)“O”型圈理論，可得沿工作面傾向方向上偏離y軸原點(diǎn)的空隙率變化系數(shù)變化函數(shù)[8]：

式中：n（y）為采空區(qū)地板處空隙率沿y軸方向變化系數(shù)，無(wú)量綱；y1為工作面傾向長(zhǎng)度，m；y為采空區(qū)內(nèi)某一點(diǎn)距進(jìn)風(fēng)巷的傾向距離，m，取值范圍［0，y1］。

在xy平面上采空區(qū)空隙率變化函數(shù)為：

在豎直方向上，將關(guān)鍵層視為臨界線，在關(guān)鍵層以下的層位，從垮落帶到裂隙帶，破碎巖體空隙率沿豎直方向呈現(xiàn)逐漸減小的分布特征；在關(guān)鍵層以上層位，由于關(guān)鍵層支撐作用，受采動(dòng)影響較小，空隙率趨近于0。由于數(shù)值模擬模型主要位于底部層位，關(guān)鍵層以上層位對(duì)于數(shù)值模擬結(jié)果影響較小，因此，只關(guān)注關(guān)鍵層以下的區(qū)域[9]。設(shè)該區(qū)域內(nèi)采空區(qū)空隙率在z方向上服從線性遞減規(guī)律，得到空隙率在xyz空間上的分布函數(shù)為：

式中：a、b為待定系數(shù)，無(wú)量綱；z為采空區(qū)內(nèi)某一點(diǎn)豎直高度，m。

式中：H為關(guān)鍵層高度，m。

將15111綜采工作面相關(guān)參數(shù)（∑h=5.9 m、m=3.5 m、kp1=1.3、l=30 m、y1=186 m、H=50 m）代入式（7）得到采空區(qū)三維空隙率分布函數(shù)：

運(yùn)用MATLAB軟件，得到采空區(qū)底板z=0處空隙率變化曲面，采空區(qū)底板垮落巖石空隙率分布如圖1。從圖1可以看出，沿x軸方向，空隙率呈遞減趨勢(shì)，這是由于采空區(qū)深部垮落巖石逐漸被壓實(shí)，從而導(dǎo)致空隙率減??；沿y軸方向，空隙率呈先減小后增大變化趨勢(shì)，這是由于兩道附近頂板受支護(hù)巷道支撐作用，使得兩側(cè)采空區(qū)的空隙率比采空區(qū)中部大。在采空區(qū)底板處，采空區(qū)垮落巖石的空隙率基本呈“滑梯”狀。

圖1 采空區(qū)底板垮落巖石空隙率分布Fig.1 Distribution of caving rock porosity in bottom plate of gob

1.2 采空區(qū)三維滲透率模型

根據(jù)Blake-Kozeny方程，得到采空區(qū)多孔介質(zhì)滲透率k與空隙率n之間關(guān)系，即：

式中：k為多孔介質(zhì)滲透率，10-15m2；Dm為多孔介質(zhì)平均粒子直徑，m；n為多孔介質(zhì)空隙率，無(wú)量綱。

1.3 采空區(qū)非均質(zhì)多孔介質(zhì)動(dòng)量損失模型

在FLUENT軟件中模擬定義采空區(qū)多孔介質(zhì)時(shí)，需要定義其黏性和慣性阻力系數(shù)作為模擬采空區(qū)對(duì)氣體的流動(dòng)阻力[10]。

式中：Si為采空區(qū)多孔介質(zhì)的動(dòng)量損失源，無(wú)量綱；μ為動(dòng)力黏度，Pa·s；ρ為流體密度，kg/m3；Dij和Cij分別為黏性阻力和慣性阻力損失系數(shù)矩陣；vj（j=1，2，3）為流體微元在x、y、z方向上的速度分量，m/s。

在黏性阻力中黏性阻力系數(shù)C1，慣性阻力中內(nèi)部慣性阻力系數(shù)為C2[11]。通過(guò)UDF導(dǎo)入FLUENT模擬軟件中，并在其中調(diào)用。

1.4 采空區(qū)流場(chǎng)控制方程

根據(jù)質(zhì)量及動(dòng)量守恒方程，建立采空區(qū)流場(chǎng)控制方程。

連續(xù)性方程：

式中：ui為在i方向的速度，m/s；xi為i方向上的坐標(biāo)值，m；t為時(shí)間，s；qm為質(zhì)量源項(xiàng)，kg/（m3·s）。

動(dòng)量守恒方程：

式中：uj為在j方向的速度，m/s；p為流體微團(tuán)上的壓力，Pa；τij為應(yīng)力張量，無(wú)量綱；xj為j方向上的坐標(biāo)值，m；fi為在單位質(zhì)量流體微團(tuán)上的i方向上體積力，N；Fi為流體微元上的體力，N[12]。

2 采空區(qū)瓦斯運(yùn)移數(shù)值模擬及結(jié)果

2.1 數(shù)值模型建立

平舒煤礦15111綜采工作面位于陽(yáng)泉煤礦15#煤東翼一采區(qū)，所采的煤層為15#煤層，煤層厚度平均為3.8 m。該綜采工作面走向長(zhǎng)為1 130.7 m，傾斜長(zhǎng)為186 m，預(yù)計(jì)回采期間瓦斯絕對(duì)涌出量為42.2 m3/min，采用“U”型通風(fēng)系統(tǒng)，核定配風(fēng)量為2 500 m3/min。15111綜采工作面采用綜合機(jī)械化一次采全高的采煤工藝。

由于實(shí)際采場(chǎng)的復(fù)雜性與不確定性，在建模時(shí)需要進(jìn)行簡(jiǎn)化。根據(jù)平舒煤礦15111綜采工作面的實(shí)際情況，工作面及進(jìn)、回風(fēng)巷的參數(shù)以專業(yè)人員現(xiàn)場(chǎng)實(shí)測(cè)為主，建立了三維采場(chǎng)物理模型，15111綜采工作面三維采場(chǎng)物理模型參數(shù)表見(jiàn)表1。運(yùn)用ICEM軟件對(duì)15111綜采工作面簡(jiǎn)化物理模型進(jìn)行非結(jié)構(gòu)化網(wǎng)格劃分，并進(jìn)行局部加密，劃分成1 432 500個(gè)網(wǎng)格單元。

表1 15111綜采工作面三維采場(chǎng)物理模型參數(shù)表Table 1 Parameter table of 3D stope physical model in 15111 fully mechanized face

工作面進(jìn)風(fēng)口設(shè)置為速度入口（inlet），進(jìn)風(fēng)風(fēng)速設(shè)置為2.89 m/s；回風(fēng)口設(shè)置為自由出口（outflow），工作面與采空區(qū)之間的面和垮落帶與裂隙帶之間的面，設(shè)為交界面（interior）[13]。在模擬過(guò)程中將采場(chǎng)流體視為不可壓流體、采空區(qū)視為各向同性多孔介質(zhì)并假定采空區(qū)瓦斯涌出源來(lái)自于采空區(qū)，并且瓦斯涌出量為定值。

此次模擬研究重力因素（有無(wú)重力）和Blake-Kozeny公式中采空區(qū)平均粒徑取值（0.08、0.12、0.16 m）對(duì)采空區(qū)瓦斯運(yùn)移及富集規(guī)律的影響[14]。

2.2 模擬結(jié)果

2.2.1 重力環(huán)境設(shè)置的影響

考慮重力環(huán)境設(shè)置，在模擬三維流場(chǎng)中很有必要，在豎直方向上，由于瓦斯密度為0.716 kg/m3，空氣密度是瓦斯密度的1.81倍，瓦斯在采空區(qū)內(nèi)流動(dòng)時(shí)，會(huì)出現(xiàn)由于重力作用的上浮效應(yīng)，使得采空區(qū)瓦斯流動(dòng)規(guī)律有所差別，采空區(qū)瓦斯體積分?jǐn)?shù)分布云圖如圖2。

圖2 采空區(qū)瓦斯體積分?jǐn)?shù)分布云圖Fig.2 Cloud charts of gas volume fraction distribution

由圖2可以看出，采空區(qū)深部和靠近回風(fēng)巷側(cè)瓦斯體積分?jǐn)?shù)較大，可知工作面向采空區(qū)漏風(fēng)主要區(qū)域?yàn)榭拷M(jìn)風(fēng)巷側(cè)工作面處。從圖2（a）、圖2（b）可以看出，由于垮落帶瓦斯涌出量較大且在無(wú)重力操作環(huán)境下，瓦斯在垮落帶深處積聚，造成采空區(qū)頂部瓦斯體積分?jǐn)?shù)比采空區(qū)底部瓦斯體積分?jǐn)?shù)小，說(shuō)明在不考慮重力因素時(shí)采空區(qū)內(nèi)瓦斯不會(huì)出現(xiàn)上浮效應(yīng)。從圖2（c）、圖2（d）可知，上隅角瓦斯積聚效應(yīng)明顯，這是因?yàn)椴煽諈^(qū)內(nèi)積存大量高體積分?jǐn)?shù)瓦斯，漏入采空區(qū)的風(fēng)流通過(guò)回風(fēng)巷帶出大量瓦斯氣體，在重力環(huán)境下，瓦斯比空氣密度低，會(huì)產(chǎn)生上浮效應(yīng)和流體分層現(xiàn)象，使得瓦斯在上隅角處以下扎姿態(tài)流入回風(fēng)巷，并在上隅角上部采空區(qū)形成“倒三角錐”區(qū)域。

在不考慮重力因素下采空區(qū)最大瓦斯體積分?jǐn)?shù)為0.876 4，在考慮重力因素下采空區(qū)最大瓦斯體積分?jǐn)?shù)為0.839 3，發(fā)現(xiàn)在無(wú)重力因素下采空區(qū)瓦斯體積分?jǐn)?shù)較高，這是因?yàn)椴豢紤]重力因素時(shí)風(fēng)流靜壓變化較小，滲流速度較低，漏風(fēng)量較小。通過(guò)以上分析，在不考慮重力時(shí)，瓦斯不會(huì)出現(xiàn)明顯的上浮效應(yīng)，這與實(shí)際采空區(qū)內(nèi)瓦斯在其中分布存在很大差異，因此在模擬采空區(qū)瓦斯運(yùn)移時(shí)考慮重力很有必要。

2.2.2 平均粒子直徑的影響

根據(jù)FLUENT模擬結(jié)果，不同平均粒子直徑下高（低）瓦斯體積分?jǐn)?shù)占比如圖3。

圖3 不同平均粒子直徑下高（低）瓦斯體積分?jǐn)?shù)占比Fig.3 Proportions of high（low）gas volume fraction under different average particle diameters

由圖3可知，當(dāng)瓦斯體積分?jǐn)?shù)為0～5%，平均粒徑為0.16 m時(shí)占總網(wǎng)格數(shù)最高為38.81%，平均粒徑為0.08 m時(shí)占比最低為37.73%，設(shè)置不同平均粒徑值對(duì)采空區(qū)低體積分?jǐn)?shù)瓦斯富集影響較小，總體上隨著平均粒徑值的增加，低體積分?jǐn)?shù)瓦斯所占區(qū)域占比增大。當(dāng)瓦斯體積分?jǐn)?shù)為5%～15%，平均粒徑為0.16 m時(shí)占總網(wǎng)格數(shù)最高為13.89%，平均粒徑為0.08 m時(shí)占比最低為8.72%，當(dāng)采空區(qū)瓦斯體積分?jǐn)?shù)5%～15%區(qū)段內(nèi)，有明火的情況下就會(huì)發(fā)生爆炸，當(dāng)采空區(qū)粒徑越大，也易發(fā)生瓦斯爆炸事故。在瓦斯體積分?jǐn)?shù)為60%～100%之內(nèi)，平均粒子直徑為0.16 m時(shí)占比最低為0。由此可以看出，隨著采空區(qū)平均粒徑的增大，采空區(qū)高體積分?jǐn)?shù)瓦斯富集區(qū)變小，低體積分?jǐn)?shù)瓦斯富集區(qū)增大。采空區(qū)平均粒徑增大使采空區(qū)整體滲透率增大，工作面向采空區(qū)漏風(fēng)量也增大，風(fēng)流從進(jìn)風(fēng)巷附近漏入采空區(qū)并從回風(fēng)巷附近漏出，帶出更多高體積分?jǐn)?shù)瓦斯，從而使得采空區(qū)整體瓦斯體積分?jǐn)?shù)下降。

不同平均粒子直徑下采空區(qū)瓦斯體積分?jǐn)?shù)分布云圖如圖4。不同平均粒子直徑下z=5 m截面處瓦斯體積分?jǐn)?shù)分布等值線如圖5。

圖4 不同平均粒子直徑下采空區(qū)瓦斯體積分?jǐn)?shù)分布云圖Fig.4 Cloud charts of gas volume fraction distribution in goaf with different average particle diameters

圖5 不同平均粒子直徑下z=5 m截面處瓦斯體積分?jǐn)?shù)分布等值線Fig.5 Contours of gas volume fraction distribution at z=5 m cross section at different mean particle diameters

由圖4可知，不同平均粒徑下采空區(qū)瓦斯分布差別較大。當(dāng)平均粒徑為0.08 m時(shí)，采空區(qū)瓦斯體積分?jǐn)?shù)范圍為0～1；平均粒徑為0.12 m時(shí)，采空區(qū)瓦斯體積分?jǐn)?shù)范圍為0～0.839；平均粒徑為0.16 m時(shí)，采空區(qū)瓦斯體積分?jǐn)?shù)范圍為0～0.569。

由圖5可知，在z=5 m截面處，不同平均粒徑下采空區(qū)整體瓦斯分布規(guī)律一致，在回風(fēng)巷側(cè)采空區(qū)深部形成高體積分?jǐn)?shù)瓦斯富集區(qū)。為了更準(zhǔn)確的看到z=5 m處的瓦斯運(yùn)移情況，布置1條沿走向的觀測(cè)線y=100 m，z=5 m，觀測(cè)采空區(qū)瓦斯體積分?jǐn)?shù)變化和采空區(qū)流體滲流速度變化。

不同平均粒子直徑下觀測(cè)線上采空區(qū)瓦斯體積分?jǐn)?shù)及滲流速度變化如圖6。

圖6 不同平均粒子直徑下觀測(cè)線上采空區(qū)瓦斯體積分?jǐn)?shù)及滲流速度變化Fig.6 Variation of seepage velocity and gas volume fraction at the observation line at different mean particle diameters in the area

從圖6可知，越遠(yuǎn)離工作面，觀測(cè)線上瓦斯體積分?jǐn)?shù)越小。不同平均粒徑下采空區(qū)瓦斯體積分?jǐn)?shù)相差較大，當(dāng)平均粒徑為0.08 m時(shí)，觀測(cè)線上最大瓦斯體積分?jǐn)?shù)為1；當(dāng)平均粒徑為0.16 m時(shí)，觀測(cè)線上最高瓦斯體積分?jǐn)?shù)為0.463，隨著平均粒徑增大，觀測(cè)線上采空區(qū)最大瓦斯體積分?jǐn)?shù)減少了約1/2。由此可見(jiàn)，采空區(qū)平均粒徑的取值對(duì)采空區(qū)瓦斯流動(dòng)影響較大。采空區(qū)初始滲流速度較小，說(shuō)明工作面漏入采空區(qū)風(fēng)量只為整體風(fēng)量的一小部分，不同平均粒子直徑下速度變化不明顯。在距離工作面0～3 m內(nèi)采空區(qū)平均粒徑越大，流體初始滲流速度越小，這是因?yàn)轱L(fēng)流漏風(fēng)位置主要位于靠近進(jìn)風(fēng)巷側(cè)工作面內(nèi)，平均粒徑越大漏入采空區(qū)的風(fēng)量越多，從而導(dǎo)致工作面中部風(fēng)流越小，滲流速度越小。在距離工作面3～32.3 m內(nèi)，采空區(qū)平均粒徑越大，流體滲流速度越小，滲流速度衰減的越慢。在距離工作面32.3～200 m，隨著風(fēng)流向采空區(qū)深部流動(dòng)，流體滲流速度逐漸減小，并在最深處達(dá)到0 m/s。

結(jié)合圖5、圖6可知，平均粒徑的取值影響著采空區(qū)漏風(fēng)情況，從而影響瓦斯運(yùn)移及富集規(guī)律，采空區(qū)平均粒徑越大，表明采空區(qū)內(nèi)煤巖空隙越大，也易引起采空區(qū)瓦斯事故。因此實(shí)際工程中，使用填充法填充采空區(qū)，減少采空區(qū)煤巖裂隙，從而減少采空區(qū)漏風(fēng)，以達(dá)到防止瓦斯事故的目的。

3 現(xiàn)場(chǎng)實(shí)測(cè)對(duì)比

為了驗(yàn)證模擬得到結(jié)果是否與實(shí)際相吻合，此次研究對(duì)比了15111工作面實(shí)際工作情況下進(jìn)回風(fēng)巷風(fēng)量和模擬得到的進(jìn)回風(fēng)巷風(fēng)量。不同平均粒徑下進(jìn)回風(fēng)巷風(fēng)量分布見(jiàn)表2。

由表2可知，不同平均粒子直徑下進(jìn)回風(fēng)巷風(fēng)量差別不大，整體規(guī)律都呈進(jìn)風(fēng)巷比回風(fēng)巷風(fēng)量大的趨勢(shì)，這是由于風(fēng)流有一部分漏入采空區(qū)所致，隨著平均粒徑的增大回風(fēng)量越來(lái)越小，說(shuō)明漏入采空區(qū)的風(fēng)量變多。對(duì)比實(shí)測(cè)結(jié)果和模擬結(jié)果，可知當(dāng)平均粒子直徑為0.12 m時(shí)誤差最小，為3.981 3%，說(shuō)明當(dāng)設(shè)定初始條件平均粒子直徑為0.12 m時(shí)，模擬結(jié)果最符合實(shí)際情況，也說(shuō)明了此次模擬較為準(zhǔn)確。

表2 不同平均粒徑下進(jìn)回風(fēng)巷風(fēng)量分布Table 2 Air volume distribution of intake and return air roadways with different average particle sizes

4結(jié) 論

1）垮落巖石空隙率在采空區(qū)底板處，空隙率呈“滑梯”狀分布；在豎直方向上，垮落巖石空隙率逐漸減小，并在關(guān)鍵層處達(dá)到最小值0。

2）重力環(huán)境設(shè)置對(duì)采空區(qū)瓦斯流動(dòng)影響較大，在有重力因素下，瓦斯會(huì)在上隅角處以下扎姿態(tài)進(jìn)入回風(fēng)巷，形成一個(gè)“倒三角錐”區(qū)域。

3）不同平均粒子直徑下，采空區(qū)瓦斯整體運(yùn)移規(guī)律較為一致，瓦斯體積分?jǐn)?shù)分布范圍有較大差別。隨著平均粒子直徑的增大，采空區(qū)高體積分?jǐn)?shù)瓦斯富集區(qū)變小，低體積分?jǐn)?shù)瓦斯富集區(qū)增大；在靠近工作面中部的采空區(qū)內(nèi)，平均粒徑增大，流體初始滲流速度減小，在采空區(qū)內(nèi)滲流速度整體呈減小趨勢(shì)。

4）通過(guò)對(duì)比平舒煤礦15111綜采工作面實(shí)測(cè)數(shù)據(jù)，考慮重力因素時(shí)模擬得到的進(jìn)回風(fēng)巷風(fēng)量變化規(guī)律與實(shí)際相符，模擬結(jié)果較為準(zhǔn)確，其中設(shè)置平均粒子直徑為0.12 m時(shí)最符合實(shí)測(cè)數(shù)據(jù)。當(dāng)模擬采空區(qū)瓦斯運(yùn)移情況時(shí)，模擬設(shè)置階段考慮重力因素和平均粒子直徑的取值，會(huì)使得采空區(qū)瓦斯運(yùn)移模擬模型更符合實(shí)際情況。