張鵬宇，陳曉坤，趙 亮，魏高明，王 奇

（1.西安科技大學(xué) 安全科學(xué)與工程學(xué)院，陜西 西安 710054；2.西安科技大學(xué) 煤火災(zāi)害防治陜西省重點實驗室，陜西 西安 710054；3.兗礦能源集團公司 鮑店煤礦，山東 鄒城 273500）

隨著現(xiàn)代化礦井開采強度與機械化程度不斷提升，礦井采空區(qū)面積和范圍隨之增加，形成大面積采空區(qū)[1]。在工作面正常推采期間，采掘活動以及頂板周期來壓使得沿空側(cè)預(yù)留煤柱壓酥破碎，形成大量裂縫，導(dǎo)致采空區(qū)之間相互貫通，為煤的氧化升溫提供漏風通道[2]。隨著大面積采空區(qū)持續(xù)供氧，采空區(qū)遺煤長時間被氧化，增大了采空區(qū)煤自燃危險性。采空區(qū)煤自燃產(chǎn)生如CH4、CO2、CO、C2H6、C2H4、C3H8、H2、SO2、H2S 等有害氣體，會在漏風通道內(nèi)運移并涌出[3-5]。因此，采空區(qū)煤自燃環(huán)境作用下氣體非控運移有可能引發(fā)瓦斯、煤塵爆炸和礦井通風系統(tǒng)紊亂等一系列嚴重的次生災(zāi)害[6-7]。采空區(qū)煤-氣-熱共生問題日益嚴峻，大面積采空區(qū)流場特性和有害氣體運移機制亟待進一步深入研究。

國內(nèi)外學(xué)者對采空區(qū)流場特征進行了大量研究，主要采用現(xiàn)場觀測、相似模擬實驗和數(shù)值模擬方法[8-10]。由于采空區(qū)條件復(fù)雜，人員無法進入其內(nèi)部，采空區(qū)隱蔽災(zāi)害信息判斷和氣體流場探測極為困難[11]。物理相似模擬實驗雖然在一定程度上解決了現(xiàn)場對于采空區(qū)垮落環(huán)境下束管監(jiān)測系統(tǒng)監(jiān)測點位破壞和數(shù)據(jù)局限性等問題，但在實驗中不可避免產(chǎn)生的誤差，無法在全域條件下研究采空區(qū)時空演化規(guī)律。而數(shù)值模擬可以彌補這一缺陷，并進一步為物理相似模擬實驗的研究提供驗證依據(jù)，為現(xiàn)場工程應(yīng)用的可行性提供指導(dǎo)[12-13]。許多研究人員已經(jīng)建立了多物理場耦合模型，旨在分析采空區(qū)中的氣體滲流[14-15]、氣體濃度場[16-17]、以及溫度場[18-19]變化規(guī)律，這些研究為采空區(qū)模擬復(fù)雜流場和氣體運移提供了必要的理論基礎(chǔ)。由于已建立數(shù)值模型通?；趩我还ぷ髅娌煽諈^(qū)的漏風規(guī)律進行研究，缺乏考慮采空區(qū)連通后氣體流場的運移規(guī)律；此外，忽略了由于采空區(qū)內(nèi)部氣體密度不同，煤自燃溫度導(dǎo)致采空區(qū)壓力以及氣體滲流和擴散能力的改變。為此，針對大面積采空區(qū)建立了一個涉及多孔介質(zhì)的氣體運輸、擴散以及能量傳遞之間復(fù)雜相互作用的多物理場耦合模型?？紤]在重力影響下，煤自燃溫度會使大面積采空區(qū)局部氣體運動變得更加復(fù)雜，根據(jù)大面積采空區(qū)壓力、溫度和氣體濃度分布規(guī)律，對采空區(qū)氣體流場形成過程進行分析，為大面積采空區(qū)的漏風控制和煤自燃防治提供科學(xué)指導(dǎo)。

1 大面積采空區(qū)多物理場耦合模型

1.1 基本守恒方程

流體運動必須遵循3 個規(guī)則，即質(zhì)量守恒、動量守恒和能量守恒[20]。

質(zhì)量守恒方程又稱為連續(xù)性方程，是指單位時間內(nèi)微元體中流體質(zhì)量的增加等于同一時間間隔內(nèi)流入該微元體的凈質(zhì)量，如式（1）：

式中：ρ 為混合氣體的密度，kg/m3；ε 為采空區(qū)孔隙度；u 為氣體滲流速度，m/s；Sm為采空區(qū)遺煤耗氧源項，kg/（m3·s）；t 為時間，s，。

動量守恒方程描述了流體系統(tǒng)所受外力與流體系統(tǒng)動量變化之間的關(guān)系，即微元體中流體動量的增加率等于作用在微元體上各種力之和，如式（2）[21]：

式中：p 為氣體靜壓，Pa；τ 為黏性應(yīng)力張量；ρgn為重力體積力，kN/m3；α 為多孔介質(zhì)滲透率，m2；μ 為氣體動力黏度，Pa·s。

能量守恒方程建立在物質(zhì)運動變化過程中控制體內(nèi)能量轉(zhuǎn)換的等量關(guān)系，決定了采空區(qū)高溫點的溫度分布和變化[22]。采空區(qū)氣體在運移過程中，溫度的變化同時也會引起采空區(qū)氣體動能和內(nèi)能改變，如式（3）：

式中：ρf、ρs為氣體和固體介質(zhì)的密度，kg/m3；Cpf、Cps為氣體和固體介質(zhì)的定壓比熱容，J/（kg·K）；kf、ks為氣體和固體介質(zhì)的導(dǎo)熱系數(shù)，W/（m·K）；T 為溫度，K；hi為組分i 的焓，kJ/mol；Ji為組分i 的擴散通量，kg/（m2·s）；Sf為遺煤的耗氧放熱源項，W/m3。

煤氧化反應(yīng)速率r 與氧氣濃度和環(huán)境溫度有關(guān)，可根據(jù)Arrhenius 方程定義為式（4）：

式中：A 為指前因子，s-1；E 反應(yīng)活化能，kJ/mol；R 為普適氣體常數(shù)，kJ/（mol·K）；CO2為氧氣濃度，kmol/m3；n 為反應(yīng)級數(shù)。

指前因子A 和反應(yīng)活化能E 可根據(jù)煤低溫氧化實驗獲得，利用Fluent 二次開發(fā)工具UDF 編寫遺煤耗氧及放熱源項。

1.2 粒子輸送方程

粒子輸送方程是指任一瞬時系統(tǒng)內(nèi)物理量（質(zhì)量、動量和能量）隨時間的變化率等于該時間控制體內(nèi)物理量的變化率與通過控制體表面凈通量之和，是流體質(zhì)點物理量在流場中在時間空間上輸運傳送變化過程[23]。粒子輸送方程如式（5）：

式中：ci為氣體組分i 的質(zhì)量分數(shù)；Qi為組分i的增減源項，kg/（m3·s）。

1.3 大面積采空區(qū)孔隙率和滲透率方程

在采空區(qū)形成過程中，采空區(qū)上覆巖層垮落破碎，垮落巖石填充采空區(qū)形成多孔介質(zhì)空間。在采空區(qū)邊界煤巖層的支撐作用和地層的垮落壓實作用下，上覆巖層的變形呈現(xiàn)出一定垮落規(guī)律，用碎脹系數(shù)Kp表征煤巖塊堆積與壓實狀況，采空區(qū)走向上碎脹系數(shù)Kp（x）如式（6）：

采空區(qū)內(nèi)部的孔隙度分布與巖石的破碎程度相關(guān)，當巖石破碎程度較小時，巖塊尺寸較大，孔隙度也變大；在采空區(qū)中心區(qū)域，上覆巖層在垮塌過程中破裂較為劇烈，在上覆巖層的自重壓實作用下，巖塊尺寸變小，孔隙度也越小。因此，有研究表明采空區(qū)孔隙度分布與邊界距離呈現(xiàn)“O”型分布形狀[24]。

滲透率是在壓差作用下，介質(zhì)中允許流體通過的能力，是用于描述流體在介質(zhì)內(nèi)傳導(dǎo)液體能力的參數(shù)。其中，滲透率越大，流體通過多孔介質(zhì)的速度越快。采空區(qū)孔隙率ε 與碎脹系數(shù)Kp如式（7），采空區(qū)滲透率α 根據(jù)Blake-Kozeny 公式計算得到[25]：

式中：Dp為垮落帶煤巖塊平均直徑，m。

2 大面積采空區(qū)數(shù)值模型

2.1 工程實例

鮑店煤礦位于山東省濟寧市，隸屬兗州煤業(yè)股份有限公司，鮑店煤礦煤樣的基礎(chǔ)參數(shù)為：水分1.910%，灰分12.790%，揮發(fā)分13.040%，固定碳52.760%，碳元素72.660%，氫元素4.274%，氮元素1.276%，臨界溫度70 °C。該地區(qū)煤樣為氣煤，易發(fā)生自燃。5316 工作面鄰近多個采空區(qū)，煤層埋藏深、通風系統(tǒng)和地質(zhì)條件復(fù)雜、采空區(qū)范圍大、內(nèi)部遺煤多、漏風通道多。5316 工作面推采過程中鄰近采空區(qū)漏風流場極為復(fù)雜，存在煤自燃隱患。以U 型通風5316 工作面采空區(qū)和鄰近五采區(qū)的5309、5310、5311 和5312 串聯(lián)采空區(qū)為物理模型進行3D 數(shù)值模擬研究。大面積采空區(qū)模型如圖1。

圖1 大面積采空區(qū)模型Fig. 1 Large area goaf model

圖1（a）是工作面的二維示意圖，工作面回采率約90%，兩側(cè)留有5 m 煤柱。5309、5310、5311 和5312 采空區(qū)已封閉，5316 為生產(chǎn)面。CFD 模型如圖1（b），大面積采空區(qū)采用六面體網(wǎng)格，為提高計算精度，對邊界層網(wǎng)格進行加密處理。

2.2 模型驗證

5309 和5310 工作面采空區(qū)預(yù)留鉆孔，通過氣相色譜監(jiān)測采空區(qū)內(nèi)部氣體體積分數(shù)，反映采空區(qū)氣體滲流情況，大面積采空區(qū)氧氣體積分數(shù)的實測值與模擬值的對比驗證如圖2。

圖2 大面積采空區(qū)氧氣體積分數(shù)的實測值與模擬值的對比驗證（風速v=1.5 m/s，高度切片Z=1 m）Fig. 2 Verification of the measured and simulated values of oxygen concentration in large area goaf（ventilation flux v=1.5 m/s, extracted slices Z=1 m）

數(shù)值模擬結(jié)果與現(xiàn)場實測O2體積分數(shù)變化規(guī)律基本一致，由于煤柱漏風，隨著埋深的增加，采空區(qū)內(nèi)O2體積分數(shù)迅速下降。風流從進風口到回風口的過程，5309 和5310 工作面O2體積分數(shù)也有所降低。對于5316 工作面采空區(qū)，當進風側(cè)采空區(qū)距工作面72 m 時，從散熱區(qū)到氧化升溫區(qū)，O2體積分數(shù)降低到18%，距工作面150 m 時，從氧化區(qū)到窒息區(qū)，O2體積分數(shù)降低至8%，氧化帶寬為78 m；回風側(cè)距工作面82 m 時，O2體積分數(shù)降低至18%，從散熱區(qū)到氧化升溫區(qū)，距工作面109 m 時，O2體積分數(shù)降低8%，從氧化區(qū)到窒息區(qū)，氧化帶寬度為27 m。

3 結(jié)果與討論

3.1 煤自燃環(huán)境下大面積采空區(qū)溫度場演化規(guī)律

大面積采空區(qū)溫度隨時間的演化云圖如圖3，大面積采空區(qū)溫度隨高度變化云圖如圖4，大面積采空區(qū)溫度分布曲線如圖5。

圖3 大面積采空區(qū)溫度隨時間的演化云圖（風速v=1.5 m/s，高度切片z=1 m）Fig.3 Cloud images of temperature evolution with time in large area goaf（ventilation flux v=1.5 m/s,extracted slices z=1 m）

圖4 大面積采空區(qū)溫度隨高度變化云圖（風速v=1.5 m/s，第30 d）Fig.4 Cloud images of temperature variation with height in large area goaf（ventilation flux v=1.5 m/s,30th day）

圖5 大面積采空區(qū)溫度分布曲線（風速v=1.5 m/s，第30 d，切片y=320 m）Fig.5 Temperature distribution curves of large area goaf（ventilation flux v=1.5 m/s, 30th day,extracted slices y=320 m）

由圖3 可知：隨著時間的推移，采空區(qū)高溫區(qū)域向四周擴散，并隨風流運動方向遷移，通過溫度的等值線可以看出溫度向外擴散和蔓延的幅度較小，與環(huán)境溫度（300 K）形成較大的溫度梯度；由于風流流經(jīng)隅角，在隅角處積聚，風流對于熱量耗散作用較小，高溫點逐漸趨向隅角處，形成高溫點向隅角偏移現(xiàn)象。

由圖4 和圖5 可知：受傳導(dǎo)和對流作用影響，溫度在采空區(qū)豎直方向上進行傳遞，高溫區(qū)域面積沒有明顯改變；在切片位置為z=1、5、10 m 處的最高溫度分別為346、339、336 K，溫度上升梯度分別為1.75 K/m 和0.6 K/m；在離開遺煤區(qū)域后，距火源距離越遠，溫度上升梯度越小，這是由于采空區(qū)屬于多孔介質(zhì)空間，內(nèi)部的巖石和空氣的導(dǎo)熱能力較差，導(dǎo)致熱量不斷積聚，從而形成高溫點并能夠持續(xù)的維持其影響范圍。

3.2 煤自燃環(huán)境下大面積采空區(qū)氣體運移規(guī)律

壓力是促使大面積采空區(qū)內(nèi)部氣體運移的直接動力，掌握采空區(qū)氧氣體積分數(shù)分布和流動規(guī)律是防治自然發(fā)火的關(guān)鍵技術(shù)基礎(chǔ)。大面積采空區(qū)不同高度下靜壓分布云圖如圖6，大面積采空區(qū)靜壓分布曲線如圖7，大面積采空區(qū)氧氣濃度分布云圖如圖8，大面積采空區(qū)氧氣濃度分布曲線如圖9，不同風速條件下大面積采空區(qū)氧氣云圖如圖10。

圖6 大面積采空區(qū)不同高度下靜壓分布云圖（風速v=1.5 m/s，第30 d）Fig.6 Cloud images of static pressure distribution at different heights in large area goaf（ventilation flux v=1.5 m/s, 30th day）

圖8 大面積采空區(qū)氧氣體積分數(shù)分布云圖（風速v=1.5 m/s，第30 d）Fig.8 Cloud images of oxygen concentration distribution in large area goaf（ventilation flux v=1.5 m/s, 30th day）

圖9 大面積采空區(qū)氧氣體積分數(shù)分布曲線（風速v=1.5 m/s，第30 d）Fig.9 Oxygen concentration distribution curves of large area goaf（ventilation flux v=1.5 m/s, 30th day）

圖10 不同風速條件下大面積采空區(qū)氧氣云圖（切片高度z=10 m）Fig.10 Oxygen cloud images of large area goaf under different wind speeds（extracted slices z=10 m）

由圖6 和圖7 可知：①無煤自燃的情況下，通風作用會在進風側(cè)形成較大壓力，在回風側(cè)形成負壓。靠近進風側(cè)采空區(qū)深部壓力隨著高度增加而增大，而靠近回風側(cè)采空區(qū)壓力隨高度增大而急劇減小，這是由于進風側(cè)風流隨漏風通道進入采空區(qū)，向采空區(qū)上部運移，導(dǎo)致采空區(qū)上部壓力增大，最終風流匯聚在回風側(cè)流出，致使壓力驟降；②煤自燃產(chǎn)生高溫會使氣體膨脹，密度減小，從而在采空區(qū)底部形成負壓區(qū)域，相應(yīng)的受浮力作用影響，氣體受熱膨脹向上運動，在采空區(qū)頂部積聚形成高壓區(qū)域。

由圖8 和圖9 可知：①氧氣進入采空區(qū)，受孔隙率分布影響，在向采空區(qū)深部擴散過程中呈現(xiàn)出不同的濃度分布，且隨著采空區(qū)豎直高度越高氧氣體積分數(shù)越低，這是由于氧氣密度大于空氣，受重力影響，導(dǎo)致氧氣在采空區(qū)底部積聚；②在煤自燃作用下，可以看到溫度產(chǎn)生的火風壓作用使局部風流方向發(fā)生改變，使得采空區(qū)煤自燃區(qū)域不斷卷吸周邊氧氣，產(chǎn)生強對流作用涌入采空區(qū)上部區(qū)域并在隅角處積聚。

加強通風會導(dǎo)致采空區(qū)漏風強度增大，進而影響采空區(qū)氣體的運移和積聚。

由圖10 可知：隨著進風風速的增加，進風側(cè)的壓力明顯增大，采空區(qū)漏風更加嚴重；因此，進風風速越大，氧化升溫帶越寬，并且明顯從工作面向采空區(qū)深部移動；隨著進風風速從0.5 m/s 增加到2.0 m/s，5309 采空區(qū)進風側(cè)氧化升溫帶從距工作面55 m移動到65 m；由于氧化升溫帶向深部轉(zhuǎn)移，采空區(qū)煤自燃的風險會更加嚴重。

4 結(jié) 語

建立了壓力場、溫度場、氣體濃度場的多場耦合模型，并將重力作用納入模擬煤自燃環(huán)境下大面積采空區(qū)模型，研究煤自燃環(huán)境下大面積采空區(qū)氣體運移規(guī)律。同時，討論了不同風速條件下O2濃度場的變化規(guī)律。

1）大面積采空區(qū)的漏風和煤自燃高溫引起的熱浮力效應(yīng)是影響采空區(qū)氣體運移主要因素。采空區(qū)發(fā)生煤自燃時，由于內(nèi)部巖石和空氣導(dǎo)熱系數(shù)較低，高度方向上溫度梯度從1.75 K/m 下降至0.6 K/m，熱量在遺煤區(qū)域積聚，引起周圍氣體的漏風和運移。

2）煤自燃高溫產(chǎn)生的火風壓作用使局部風流方向發(fā)生改變，使得采空區(qū)煤自燃區(qū)域不斷卷吸周邊氧氣，產(chǎn)生強對流作用涌入采空區(qū)上部區(qū)域并在隅角處積聚。

3）隨著漏風源風速的增大，大面積采空區(qū)進風側(cè)氧化升溫帶從距工作面55 m 移動到65 m，增大了大面積采空區(qū)煤自燃的風險，同時加快了采空區(qū)深部有害氣體向工作面的運移。