基于三維空隙率模型的采空區瓦斯運移及富集規律

2021-06-03 09:32:46曹明月李小芳孫浩石閆志銘

煤礦安全 2021年5期

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

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

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

1 采空區三維流場數學模型

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

1.1 采空區三維空隙率模型

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

由煤巖體碎脹系數定義可知采空區內破碎巖體的碎脹系數沿x軸的分布為：

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

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

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

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

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

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

在xy平面上采空區空隙率變化函數為：

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

式中：a、b為待定系數，無量綱；z為采空區內某一點豎直高度，m。

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

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

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

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

1.2 采空區三維滲透率模型

根據Blake-Kozeny方程，得到采空區多孔介質滲透率k與空隙率n之間關系，即：

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

1.3 采空區非均質多孔介質動量損失模型

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

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

在黏性阻力中黏性阻力系數C1，慣性阻力中內部慣性阻力系數為C2[11]。通過UDF導入FLUENT模擬軟件中，并在其中調用。

1.4 采空區流場控制方程

根據質量及動量守恒方程，建立采空區流場控制方程。

連續性方程：

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

動量守恒方程：

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

2 采空區瓦斯運移數值模擬及結果

2.1 數值模型建立

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

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

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

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

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

2.2 模擬結果

2.2.1 重力環境設置的影響

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

圖2 采空區瓦斯體積分數分布云圖Fig.2 Cloud charts of gas volume fraction distribution

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

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

2.2.2 平均粒子直徑的影響

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

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

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

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

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

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

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

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

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

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

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

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

3 現場實測對比

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

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