趙文彬，石新巖，張培偉，李振武，張延松，劉方順

(1.山東科技大學(xué) 安全與環(huán)境工程學(xué)院，山東 青島 266590;2.山東濟(jì)寧運(yùn)河煤礦有限責(zé)任公司，山東 濟(jì)寧 272000)

0 引言

隨煤礦開(kāi)采深度增加，受地應(yīng)力、地溫、熱環(huán)境及采空區(qū)空隙率結(jié)構(gòu)變化等因素影響，采空區(qū)自燃帶范圍及遺煤自燃規(guī)律發(fā)生變化。

學(xué)者針對(duì)采空區(qū)自燃區(qū)域劃分展開(kāi)研究：徐精彩等[1]通過(guò)測(cè)定采空區(qū)氧濃度分布狀況，提出采空區(qū)遺煤自燃極限參數(shù)計(jì)算方法，構(gòu)建自燃危險(xiǎn)區(qū)域判定充分條件；李宗翔等[2]利用模擬軟件分析漏風(fēng)對(duì)采空區(qū)溫度場(chǎng)影響，指出自燃“三帶”呈非對(duì)稱性分布；李樹(shù)剛等[3]得出采空區(qū)上覆巖層空隙結(jié)構(gòu)變化規(guī)律理論公式；司俊鴻等[4]建立采空區(qū)空隙率及滲透率三維分布模型，提出采空區(qū)空隙結(jié)構(gòu)在立體空間服從線性規(guī)律；王少鋒等[5]研究采空區(qū)及覆巖空隙率三維分布特點(diǎn)，提出考慮煤層傾角的空隙率變化公式。

學(xué)者一般采用數(shù)值模擬與現(xiàn)場(chǎng)實(shí)測(cè)結(jié)合的方法，研究采空區(qū)自燃區(qū)域劃分影響因素：車強(qiáng)[6]利用CFD模擬得出遺煤自燃前后采空區(qū)氣體組分三維滲流場(chǎng)、濃度場(chǎng)和溫度場(chǎng)多場(chǎng)耦合變化規(guī)律；李宗翔等[7]采用煤樣封閉耗氧實(shí)驗(yàn)得到煤樣隨時(shí)間變化耗氧速率，基于Fluent模擬軟件研究曉南礦通風(fēng)方式對(duì)采空區(qū)自燃發(fā)火的影響。

煤層自燃受多種因素影響，煤層深部受地溫、壓力、空隙結(jié)構(gòu)及煤氧結(jié)合能力影響，采空區(qū)自燃帶分布發(fā)生改變。本文結(jié)合魯西南千米深井、厚煤層、傾斜工作面采空區(qū)特點(diǎn)，基于采空區(qū)空間空隙結(jié)構(gòu)、溫度場(chǎng)及耗氧速率變化規(guī)律，在理論分析基礎(chǔ)上利用Fluent軟件重點(diǎn)模擬不同通風(fēng)方式下采空區(qū)熱環(huán)境及耗氧速率對(duì)工作面流場(chǎng)及煤自燃區(qū)域分布的影響，并結(jié)合現(xiàn)場(chǎng)實(shí)測(cè)結(jié)果對(duì)模擬數(shù)據(jù)進(jìn)行對(duì)比分析。

1 傾斜工作面采空區(qū)空隙率分布規(guī)律

角度及空隙結(jié)構(gòu)影響采空區(qū)風(fēng)流擴(kuò)散，進(jìn)而影響采空區(qū)自燃區(qū)域劃分。

受地應(yīng)力影響，采空區(qū)上部圍巖垮落,破碎巖體散落到采空區(qū)形成冒落松散結(jié)構(gòu)[8-10]。巖石間裂隙分布在整個(gè)采空區(qū)，采空區(qū)空隙結(jié)構(gòu)變化如圖1所示。隨工作面推進(jìn)，采空區(qū)上覆圍巖形成裂隙帶、彎曲變形帶；冒落帶自然垮落，在采空區(qū)形成c-自然堆積區(qū)、b-載荷影響區(qū)、a-壓實(shí)區(qū)。沿工作面推進(jìn)方向，采空區(qū)受垮落影響，工作面附近空隙結(jié)構(gòu)較大；隨采空區(qū)距離增加，采空區(qū)上覆巖層冒落壓實(shí)，空隙結(jié)構(gòu)變小，并在后部a區(qū)逐漸穩(wěn)定。

圖1 采空區(qū)空隙結(jié)構(gòu)變化Fig.1 Change of void structure in goaf

據(jù)調(diào)研可知,采空區(qū)傾角空隙率可基本反映現(xiàn)場(chǎng)實(shí)際，結(jié)合采空區(qū)冒落帶空隙率連續(xù)變化方程，得到采空區(qū)空隙率如式(1)所示：

(1)

式中：φG為采空區(qū)底部走向中軸線的空隙率；hd為直接頂高度，m；H為采高，m；kpb為巖體碎脹系數(shù)；l為基本頂巖塊長(zhǎng)度，m；x為采空區(qū)深部走向長(zhǎng)度，m。

受回采工藝影響，傾斜煤層采空區(qū)上下巷端頭空隙高于中間位置，采空區(qū)漏風(fēng)“O”形圈結(jié)構(gòu)愈發(fā)明顯。由“O”形圈理論可知，進(jìn)回風(fēng)巷附近采空區(qū)空隙率大于中部，空隙結(jié)構(gòu)如式(2)所示：

(2)

式中：φG,y為工作面傾向(y軸方向)上偏移坐標(biāo)原點(diǎn)的空隙結(jié)構(gòu)變化系數(shù)；Y為沿工作面走向與模型坐標(biāo)原點(diǎn)偏移距離，m。

沿垂直方向，受重力場(chǎng)作用采空區(qū)空隙率隨冒落高度增加而增大[11-13]。采空區(qū)破碎巖體重力作用示意如圖2所示。由圖2可知，當(dāng)工作面存在傾角時(shí)，采空區(qū)內(nèi)煤矸石受重力作用，底部壓實(shí)上部松散，根據(jù)破碎巖體空隙率軸向應(yīng)力變化規(guī)律及煤巖體自身重力[14]，得到空隙結(jié)構(gòu)如式(3)所示：

圖2 采空區(qū)破碎巖體重力作用示意Fig.2 Schematic diagram for gravity action of broken rock mass in goaf

(3)

式中：φγ為破碎巖體受軸向應(yīng)力后的空隙率；β0為破碎巖體受軸向應(yīng)力前的空隙率；βi為相關(guān)系數(shù)；γ為垮落巖石容重，N/m3；α為工作面傾角,(°)；Y為沿工作面走向與模型坐標(biāo)原點(diǎn)偏移距離，m；ly為工作面長(zhǎng)度，m。

模擬參數(shù)設(shè)置見(jiàn)表1。根據(jù)式(1)～(3)以及表1，利用Python軟件選取距底板高度Z=1，5，7 m時(shí)，各區(qū)域立體空間采空區(qū)空隙率變化如圖3所示。

表1 模擬參數(shù)Table 1 Simulation parameters

圖3 立體空間采空區(qū)空隙率變化Fig.3 Change of porosity in three-dimensional goaf

由圖3可知，空隙結(jié)構(gòu)沿工作面及采深方向呈“凹”形分布，靠近工作面和進(jìn)回風(fēng)巷兩側(cè)空隙結(jié)構(gòu)較大，中部及采空區(qū)深部空隙結(jié)構(gòu)較小。受重力影響，隨采空區(qū)高度Z增加，空隙結(jié)構(gòu)變大；受工作面傾角影響，上巷空隙結(jié)構(gòu)大于下巷。

2 采空區(qū)溫度場(chǎng)對(duì)耗氧速率影響

基于阿倫尼烏斯定律，溫度對(duì)氧氣消耗速率影響如式(4)所示：

Vo2=ACnexp(-E/RT)

(4)

式中：A為指前因子，也稱阿倫尼烏斯常數(shù)；E為活化能，kJ/mol；T為絕對(duì)溫度，K；R為氣體常數(shù)，kJ/(mol·K)；C為氧氣體積分?jǐn)?shù),%；n為常數(shù),取0.5～1。

氣體在采空區(qū)流動(dòng)時(shí)，氧氣與遺煤發(fā)生化學(xué)反應(yīng)，導(dǎo)致氧氣濃度降低。煤樣封閉耗氧實(shí)驗(yàn)得到不同氧濃度下煤的耗氧速率，如式(5)所示：

V=-λc(cτ-cb)

(5)

式中：cb為下限氧濃度，mol/m3；λc為氧氣濃度衰減率；τ為氧化時(shí)間，s；cτ為空氣氧化后的氧濃度，mol·m-3。

深部煤層因采空區(qū)暴露遺煤、圍巖溫度高，使煤氧結(jié)合能力增大，耗氧量增加。模擬過(guò)程中，基于式(4)～(5)，重點(diǎn)考慮采空區(qū)熱環(huán)境條件下，煤耗氧速率對(duì)采空區(qū)自燃區(qū)域影響。

3 礦井及模型概況

兗州礦區(qū)某礦工作面位于31采區(qū)東北方向，煤層埋深1 050 m，原巖溫度37 ℃，煤層厚度8.6 m。工作面采取U型通風(fēng)方式，走向長(zhǎng)度150 m,設(shè)計(jì)采高3.5 m，工作面風(fēng)量1 200 m3/min。該礦井屬于低瓦斯礦井，有自燃發(fā)火傾向性，屬自燃煤層，自燃發(fā)火期52 d；工作面傾角較大，地溫梯度3.2 ℃/100 m。

深部?jī)A斜采空區(qū)模型如圖4所示。其中，x軸方向?yàn)椴煽諈^(qū)走向、y軸方向?yàn)楣ぷ髅孀呦颉軸與采空區(qū)底部垂直布局建立坐標(biāo)系。模擬主要參數(shù)見(jiàn)表2。

表2 模擬參數(shù)設(shè)置Table 2 Setting of simulation parameters

圖4 兗州礦區(qū)某礦工作面采空區(qū)模型Fig.4 Goaf model of working face in a certain mine of Yanzhou mining area

4 工作面風(fēng)流溫度對(duì)采空區(qū)熱場(chǎng)分布影響分析

根據(jù)兗州礦區(qū)實(shí)測(cè)數(shù)據(jù)，工作面夏季進(jìn)風(fēng)巷風(fēng)溫18 ℃，回風(fēng)巷風(fēng)溫29 ℃。結(jié)合表1～2設(shè)置模擬參數(shù)，利用Fluent軟件用戶自定義(UDF)功能編寫(xiě)空隙率變化方程。模擬不同通風(fēng)方式下工作面風(fēng)流溫度對(duì)采空區(qū)熱場(chǎng)影響范圍，如圖5～6所示。

對(duì)比圖5與圖6可知，傾斜工作面上巷采空區(qū)空隙結(jié)構(gòu)較大，下巷空隙結(jié)構(gòu)小，導(dǎo)致下行風(fēng)進(jìn)風(fēng)側(cè)冷風(fēng)擴(kuò)散范圍大于上行風(fēng)進(jìn)風(fēng)側(cè)。

圖5 下行風(fēng)風(fēng)溫對(duì)采空區(qū)熱場(chǎng)影響Fig.5 Influence of downward wind temperature on thermal field of goaf

圖6 上行風(fēng)風(fēng)溫對(duì)采空區(qū)熱場(chǎng)影響Fig.6 Influence of upward wind temperature on thermal field of goaf

對(duì)比工作面中下部溫度場(chǎng)變化可知，受采空區(qū)內(nèi)部地溫溫差(根據(jù)煤層傾角和地溫梯度計(jì)算，上下巷原巖溫差2.5 ℃)引起的浮升力影響，下行風(fēng)時(shí)，進(jìn)風(fēng)側(cè)新鮮冷空氣與采空區(qū)內(nèi)部浮升氣流相互影響，內(nèi)部熱風(fēng)與工作面冷風(fēng)交匯后導(dǎo)致采空區(qū)熱風(fēng)涌向工作面。由于浮升力較小，2者在工作面中下部產(chǎn)生壓力平衡區(qū)，熱風(fēng)涌出造成工作面中下部溫度升高；上行風(fēng)風(fēng)流方向與采空區(qū)浮升力方向一致，采空區(qū)熱風(fēng)涌向工作面現(xiàn)象不明顯。

上行風(fēng)回風(fēng)側(cè)由于采空區(qū)空間空隙結(jié)構(gòu)較大，利于氣體流動(dòng)，導(dǎo)致回風(fēng)側(cè)溫度升高；下行風(fēng)回風(fēng)側(cè)空隙結(jié)構(gòu)小，風(fēng)流受采空區(qū)風(fēng)流影響較小，風(fēng)流溫度較低。

5 采空區(qū)熱環(huán)境對(duì)煤自燃區(qū)域影響模擬分析

基于設(shè)定采空區(qū)初始溫度條件，結(jié)合前期熱重實(shí)驗(yàn)結(jié)果，考慮煤低溫氧化在該溫度條件下耗氧速率，通過(guò)UDF編程模擬該溫度條件對(duì)煤自燃區(qū)域影響。模擬中以O(shè)2體積分?jǐn)?shù)的5%、15%作為自燃帶范圍劃分依據(jù)。

5.1 采空區(qū)熱環(huán)境對(duì)自燃帶范圍影響模擬

為觀察多條件采空區(qū)熱環(huán)境對(duì)O2濃度擴(kuò)散影響，各取上下行風(fēng)采空區(qū)O2濃度變化截面，得到不同條件下上行風(fēng)、下行風(fēng)自燃帶范圍模擬結(jié)果如圖7～8所示。受地溫梯度影響，由于工作面傾角為30°，上下巷高差為75 m，需將地溫梯度變化方程編寫(xiě)進(jìn)采空區(qū)溫度計(jì)算區(qū)域，溫差取2.5 ℃。

由圖7可知，隨采空區(qū)熱場(chǎng)條件增加，上行風(fēng)時(shí)采空區(qū)自燃帶范圍整體擴(kuò)大，回風(fēng)側(cè)變化大于進(jìn)風(fēng)側(cè)。對(duì)比圖7(a)和圖7(b)可知，由進(jìn)風(fēng)側(cè)散熱帶明顯突出、工作面中下部及回風(fēng)側(cè)氧化帶后移5～10 m。

圖7 不同條件下上行風(fēng)自燃帶模擬對(duì)比Fig.7 Simulation contrast of spontaneous combustion zone with upward wind under different conditions

由圖8可知，隨采空區(qū)熱場(chǎng)條件疊加增加，采空區(qū)自燃帶范圍整體縮小，回風(fēng)側(cè)變化小于進(jìn)風(fēng)側(cè)；對(duì)比圖8(a)和圖8(b)可知，進(jìn)風(fēng)側(cè)氧化帶范圍繼續(xù)縮小，“前凸”現(xiàn)象明顯、靠近工作面散熱帶范圍明顯擴(kuò)大。

圖8 不同條件下下行風(fēng)自燃帶模擬對(duì)比Fig.8 Simulation contrast of spontaneous combustion zone with downward wind under different conditions

5.2 現(xiàn)象分析

由模擬結(jié)果可知，深部?jī)A斜厚煤層受工作面冷風(fēng)流、采空區(qū)內(nèi)部圍巖散熱及浮升力等綜合影響，導(dǎo)致溫度場(chǎng)發(fā)生改變，進(jìn)而影響各區(qū)域煤的氧化活性(耗氧量)，使采空區(qū)自燃區(qū)域范圍發(fā)生改變。主要規(guī)律概括為以下2點(diǎn)：

1)受深部礦井傾斜采空區(qū)空隙結(jié)構(gòu)影響，上巷空隙結(jié)構(gòu)大，漏風(fēng)多，使上行風(fēng)進(jìn)風(fēng)側(cè)自燃帶變化范圍小于回風(fēng)側(cè)，下行風(fēng)進(jìn)風(fēng)側(cè)自燃帶變化范圍大于回風(fēng)側(cè)。

2)受深部礦井采空區(qū)內(nèi)部熱源影響，采空區(qū)熱向工作面涌出，高溫影響耗氧速率的變化。

上行風(fēng)時(shí)，浮升力與風(fēng)流方向相同，有利于新鮮風(fēng)流在采空區(qū)擴(kuò)散，采空區(qū)自燃帶范圍擴(kuò)大、工作面中部及下巷附近“前凸”現(xiàn)象明顯。下行風(fēng)時(shí)，浮升力與風(fēng)流方向相反，抑制工作面風(fēng)流流向采空區(qū)，使采空區(qū)內(nèi)部流場(chǎng)紊亂，導(dǎo)致采空區(qū)自燃帶范圍縮小。

5.3 浮升力對(duì)采空區(qū)立體自燃帶影響研究

為觀察該區(qū)域氧氣濃度變化，沿采空區(qū)深部走向，在距進(jìn)回風(fēng)巷道5 m處各取1個(gè)截面，得到不同通風(fēng)條件下進(jìn)回風(fēng)側(cè)O2體積分?jǐn)?shù)變化，如圖9所示。沿截面各取3條監(jiān)測(cè)直線，得到不同通風(fēng)條件下氧氣濃度變化如圖10～11所示。

圖9 不同通風(fēng)方式下進(jìn)回風(fēng)側(cè)O2體積分?jǐn)?shù)變化Fig.9 Change of O2 volume fraction in inlet and outlet sides under different ventilation modes

圖10 上行風(fēng)進(jìn)回風(fēng)巷道采空區(qū)O2體積分?jǐn)?shù)變化折線Fig.10 Variation broken line of O2 volume fraction in goaf at inlet and return air duct with upward wind

由圖9(a)監(jiān)測(cè)面1、2可知，進(jìn)風(fēng)下巷周圍O2濃度較低，采空區(qū)自燃帶靠近工作面且范圍較??；回風(fēng)側(cè)自燃帶范圍后移且區(qū)域增大；工作面附近O2濃度較大。O2濃度變化邊界線呈“橢圓狀”，采空區(qū)O2濃度分布呈中間低，頂部及底部O2濃度高的趨勢(shì)；遠(yuǎn)離工作面方向，O2濃度較低，且剖面變化弧度方向相反。

由圖9(b)監(jiān)測(cè)面1、2可知，進(jìn)風(fēng)上巷O2分布范圍大于回風(fēng)側(cè)，回風(fēng)下巷O2分布范圍前移；進(jìn)風(fēng)側(cè)前部O2濃度值小于回風(fēng)側(cè)前部。下行風(fēng)進(jìn)風(fēng)側(cè)空隙結(jié)構(gòu)大，氧氣擴(kuò)散范圍大；受耗氧速率影響，氧氣濃度值較低?；仫L(fēng)側(cè)空隙結(jié)構(gòu)小，氧氣擴(kuò)散范圍??；受內(nèi)部浮升力作用，回風(fēng)側(cè)氧氣濃度值升高。

由圖10可知，受工作面風(fēng)流影響，進(jìn)風(fēng)側(cè)在18 m處進(jìn)入氧化帶，回風(fēng)側(cè)在23 m處進(jìn)入氧化帶；回風(fēng)側(cè)空隙結(jié)構(gòu)較大，O2消耗速率較快，自燃帶處于35～55 m范圍；進(jìn)風(fēng)側(cè)空隙結(jié)構(gòu)較小，氣體不易流動(dòng)，O2消耗速率緩慢，自燃帶處于25～80 m。

由圖11可知，進(jìn)回風(fēng)側(cè)均在18 m處進(jìn)入氧化帶；進(jìn)風(fēng)側(cè)空隙結(jié)構(gòu)較大，自燃帶處于18～45 m；由于回風(fēng)側(cè)空隙結(jié)構(gòu)較小，且受浮升力作用，自燃帶范圍縮小為18～40 m。

圖11 下行風(fēng)進(jìn)回風(fēng)巷道采空區(qū)O2體積分?jǐn)?shù)變化折線Fig.11 Variation broken line of O2 volume fraction in goaf at inlet and return air duct with downward wind

由圖10～11可知，當(dāng)處于采空區(qū)同一位置時(shí)，中部監(jiān)測(cè)直線O2體積分?jǐn)?shù)變化值大于其他監(jiān)測(cè)直線，表明受溫度影響下采空區(qū)浮升力及耗氧速率影響，采空區(qū)濕熱氣流升騰作用使自燃帶范圍在垂直方向產(chǎn)生差異。

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

兗州礦區(qū)某礦工作面采用上行通風(fēng)方式已推進(jìn)一段距離。為了解工作面正常推進(jìn)過(guò)程中自燃發(fā)火區(qū)域基本狀況，分別在工作面進(jìn)回風(fēng)巷道布置A，B2個(gè)測(cè)點(diǎn)，如圖12所示。由于安設(shè)在工作面軌道巷測(cè)點(diǎn)無(wú)法正常保護(hù)，主要觀測(cè)回風(fēng)巷道指標(biāo)參數(shù)。束管監(jiān)測(cè)系統(tǒng)觀測(cè)時(shí)間從2019年7月24日至2019年9月17日，時(shí)間為56 d。獲取井下同一時(shí)期內(nèi)日推進(jìn)距離及采空區(qū)內(nèi)不同監(jiān)測(cè)日期O2濃度值，如圖13所示。

圖12 工作面推進(jìn)過(guò)程中測(cè)點(diǎn)布置Fig.12 Arrangement of measuring points in advancing process of working face

圖13 不同監(jiān)測(cè)日期的O2體積分?jǐn)?shù)Fig.13 O2 volume fractions at different monitoring dates

由圖13可知，正常推進(jìn)過(guò)程中，2019年7月28日之前O2濃體積分?jǐn)?shù)在15%以上；2019年8月7日以后，降至5%以下。

結(jié)合該時(shí)期內(nèi)推進(jìn)距離，計(jì)算回風(fēng)側(cè)自燃帶實(shí)測(cè)范圍19～47 m，回風(fēng)側(cè)自燃帶模擬范圍21～52 m。該工作面已實(shí)現(xiàn)智能化開(kāi)采，受巷道支護(hù)加強(qiáng)影響，窒息帶實(shí)測(cè)范圍后移；由于采空區(qū)遺煤量、工作面風(fēng)流速度、深部礦井周圍環(huán)境等因素對(duì)O2擴(kuò)散影響，實(shí)測(cè)回風(fēng)側(cè)氧化帶范圍大于模擬范圍，測(cè)量誤差小于等于3 m。

7 結(jié)論

1)深部?jī)A斜厚煤層在采動(dòng)過(guò)程中，巷道及采空區(qū)空隙結(jié)構(gòu)隨采深增加而減小，隨高度增加而變大；通過(guò)對(duì)煤自燃理論分析，指出采空區(qū)內(nèi)部耗氧速率受深部煤層溫度場(chǎng)影響，呈不均衡性。

2)受工作面風(fēng)流溫度影響，下行風(fēng)在進(jìn)風(fēng)巷附近對(duì)采空區(qū)熱場(chǎng)影響大于上行風(fēng)；上行風(fēng)工作面中后部受采空區(qū)熱場(chǎng)影響大于下行風(fēng)；下行冷風(fēng)與浮升力攜帶的熱空氣在工作面中下部匯合，將攜帶熱風(fēng)涌向工作面。

3)受采空區(qū)空隙結(jié)構(gòu)、耗氧速率及內(nèi)部溫度場(chǎng)作用，不同通風(fēng)方式對(duì)采空區(qū)自燃帶分布規(guī)律影響較大，其中上行風(fēng)自燃帶范圍大于下行風(fēng)。

4)上行通風(fēng)自燃帶進(jìn)風(fēng)側(cè)范圍25～40 m，回風(fēng)側(cè)范圍21～52 m；下行通風(fēng)自燃帶進(jìn)風(fēng)側(cè)范圍15～40 m，回風(fēng)側(cè)范圍15～24 m。結(jié)合工作面實(shí)測(cè)參數(shù)得出模擬與實(shí)際誤差較小，結(jié)果較為吻合。