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

（1.中國礦業(yè)大學（北京）共伴生能源精準開采北京市重點實驗室，北京100083；2.中國礦業(yè)大學（北京）應急管理與安全工程學院，北京100083；3.山東科技大學 礦山災害預防控制重點實驗室，山東 青島266590；4.陽泉煤業(yè)（集團）有限責任公司，山西 陽泉045000）

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

1 采空區(qū)三維流場數(shù)學模型

煤層開采后采場中形成矩形采動空間，受二次應力影響，圍巖發(fā)生垮落、斷裂和變形，形成了采場上覆巖體結構的“砌體梁”模型[6]?；凇捌鲶w梁”理論和“O”型圈理論，提出如下模型：取進風巷側(cè)工作面與采空區(qū)的交點為坐標原點O，x軸方向為沿采空區(qū)走向方向，y軸方向為沿采空區(qū)傾向方向，z軸方向為沿采空區(qū)豎直方向。

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

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

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

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

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

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

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

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

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

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

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

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

式中：H為關鍵層高度，m。

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

運用MATLAB軟件，得到采空區(qū)底板z=0處空隙率變化曲面，采空區(qū)底板垮落巖石空隙率分布如圖1。從圖1可以看出，沿x軸方向，空隙率呈遞減趨勢，這是由于采空區(qū)深部垮落巖石逐漸被壓實，從而導致空隙率減小；沿y軸方向，空隙率呈先減小后增大變化趨勢，這是由于兩道附近頂板受支護巷道支撐作用，使得兩側(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之間關系，即：

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

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

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

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

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

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

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

連續(xù)性方程：

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

動量守恒方程：

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

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

2.1 數(shù)值模型建立

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

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

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

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

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

2.2 模擬結果

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

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

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

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

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

2.2.2 平均粒子直徑的影響

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

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

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

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

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

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

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

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

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

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

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

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

3 現(xiàn)場實測對比

為了驗證模擬得到結果是否與實際相吻合，此次研究對比了15111工作面實際工作情況下進回風巷風量和模擬得到的進回風巷風量。不同平均粒徑下進回風巷風量分布見表2。

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

表2 不同平均粒徑下進回風巷風量分布Table 2 Air volume distribution of intake and return air roadways with different average particle sizes

4結 論

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

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

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

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