秦劍云 ，宋雙林 ，王永敬 ，李小超 ，黨 龍 ，朱 鵬 ，李世豪

（1.庫車縣榆樹嶺煤礦有限責(zé)任公司，新疆 庫車 842099；2.中煤科工集團(tuán)沈陽研究院有限公司，遼寧 撫順 113122；3.煤礦安全技術(shù)國家重點實驗室，遼寧 撫順 113122）

煤炭作為我國主體能源，其安全供給直接關(guān)系到國家能源安全[1]。煤礦火災(zāi)是威脅煤礦安全生產(chǎn)的主要原因之一，已嚴(yán)重威脅到煤礦企業(yè)的安全生產(chǎn)，其中煤炭自燃是引發(fā)煤礦火災(zāi)的最主要原因[2-3]，采動影響下的垮落巖體和遺煤在地下形成的多孔松散采空區(qū)是瓦斯與煤自燃致災(zāi)的重災(zāi)區(qū)，井下采空區(qū)遺煤自燃占井下火災(zāi)事故的60%[4]。

針對采空區(qū)的自然發(fā)火問題，國內(nèi)外學(xué)者開展了一系列的研究，包括防滅火措施、自燃規(guī)律的數(shù)值模擬及實驗分析等。在數(shù)值模擬方面，前期的研究以研究風(fēng)流在采場的穩(wěn)態(tài)分布規(guī)律為主[5-7]，根據(jù)求解的漏風(fēng)速度對自燃區(qū)域進(jìn)行大致區(qū)分。目前的研究多是針對采空區(qū)自然發(fā)火的多場耦合模擬研究[8-9]，將流場、溫度場、濃度場、應(yīng)力場等進(jìn)行耦合，從礦井尺度上研究煤自燃高溫區(qū)域的發(fā)生及發(fā)展的規(guī)律[10-12]。由于采空區(qū)的物理尺寸和物理特性都會隨著工作面的推進(jìn)而改變，一些學(xué)者通過引入移動坐標(biāo)系來處理物理邊界的變化[8,13]，或者利用CFD 軟件的動網(wǎng)格技術(shù)[14]，得到孔隙度在時間和空間上的非均質(zhì)分布，體現(xiàn)高溫區(qū)域在小范圍推進(jìn)距離下的動態(tài)移動規(guī)律，模擬結(jié)果更符合實際。文獻(xiàn)[15-17]通過COMSOL Multiphysics?與MATLAB 聯(lián)用，在生成新的采空區(qū)網(wǎng)格后，調(diào)用上一時間步的場信息作為當(dāng)前時間步的初值，實現(xiàn)模擬的物理量在時間和空間上的連續(xù)性，較大程度呈現(xiàn)了采空區(qū)煤自燃升溫過程，但這些研究未將工作面與采空區(qū)的流動、傳熱及傳質(zhì)問題聯(lián)合求解，從而也缺少圍巖溫度、通風(fēng)溫度等因素對采空區(qū)煤自燃升溫過程影響的研究。為此，考慮采空區(qū)動態(tài)演化過程，以巷道出入口作為邊界，將工作面與采空區(qū)的物理場聯(lián)合求解，模擬工作面動態(tài)推進(jìn)下采空區(qū)煤巖裂隙場中煤自燃災(zāi)害的發(fā)生、發(fā)展及演化規(guī)律；通過對采空區(qū)遺煤氧化升溫情況的預(yù)測預(yù)報，可確定采空區(qū)煤自燃的危險性和高溫?zé)嵩吹奈恢茫M(jìn)而為采空區(qū)這一隱蔽火區(qū)煤自燃的治理提供及時有效的手段。

1 模型構(gòu)建

U 型采場可分為幾個區(qū)域，即煤體固相區(qū)、采空多孔介質(zhì)區(qū)和采煤作業(yè)區(qū)。煤體固相區(qū)包括面間煤柱、邊界煤柱、待回采實體煤、頂板巖層和底板巖層；采煤作業(yè)區(qū)包括進(jìn)風(fēng)巷、工作面和回風(fēng)巷。采空區(qū)二維示意圖如圖1。

圖1 采空區(qū)二維示意圖Fig.1 Two-dimensional schematic diagram of the mining area

1.1 基本假設(shè)

1）因采空區(qū)長度和寬度遠(yuǎn)大于煤層開挖引起的上覆煤巖層垮落高度，將采空區(qū)視為高度上平均的二維非均質(zhì)多孔介質(zhì)。

2）假設(shè)氣體為理想氣體，并忽略水蒸氣對煤自熱過程的影響。

1.2 采空區(qū)孔隙率及滲透率的動態(tài)演化

隨著工作面向前推進(jìn)，自由垮落帶的矸石、遺煤在重力和礦壓作用下會被逐漸壓實，垮落碎脹系數(shù)在壓實過程中不斷減小，其數(shù)值滿足式（1）[18]：

式中：Kp(x,y)為 采空區(qū)垮落碎脹系數(shù)；Kp,max為初始垮落碎脹系數(shù)；Kp,min為垮落壓實后的碎脹系數(shù)； αx、 αy分別為采空區(qū)碎脹系數(shù)傾向和走向方向的衰減率，1/m；dx、dy分別為采空區(qū)任意點(x,y)到采空區(qū)某一固定圍巖邊界和工作面邊界的距離，m； ξ1為控制“O”型圈模型分布形態(tài)的調(diào)整系數(shù)， ξ1<1。

如果工作面推進(jìn)速度為u，則隨著時間推進(jìn)，采空區(qū)內(nèi)任意一點到采空區(qū)上下邊界和工作面的距離分別可表示為：

式中：u為工作面推進(jìn)速度，m/s；t為工作面回采時間，d；L為工作面長度，m；x0為工作面初始位置，m。

垮落采空區(qū)的高度H和孔隙率 εb的 表達(dá)式分別為：

式中：M為工作面采高，m； εb為采空區(qū)孔隙率。

根據(jù)Ergun 的實驗研究[19]，采空區(qū)滲透率K可以根據(jù)孔隙率、體積平均粒徑計算得到，計算公式為：

式中：dp為采空區(qū)多孔介質(zhì)的平均顆粒直徑，m。

1.3 采空區(qū)流動及傳熱方程

采空區(qū)內(nèi)滲透率具有明顯的非均質(zhì)特征，固壁處的滲透率和采空區(qū)深處的孔隙率相差1 個數(shù)量級，研究表明在靠近工作面的采空區(qū)內(nèi)具有較大的空氣流速，空氣處于湍流狀態(tài)，而采空區(qū)內(nèi)部的空氣流速較低，空氣處于層流或者蠕動流狀態(tài)，因此使用Darcy-Brinkman 方程描述采空區(qū)內(nèi)氣體的流動[20]。對于采空區(qū)內(nèi)空氣與煤巖體之間的換熱過程，采用局部熱平衡假設(shè)，即假定局部流體溫度與煤巖體溫度一致。采空區(qū)內(nèi)的流動、傳熱及傳質(zhì)方程[21-23]為：

氧氣消耗量與氧氣濃度和煤自燃參數(shù)有關(guān)，煤氧反應(yīng)熱與氧氣消耗量和反應(yīng)熱有關(guān)，分別計算如下：

式中：RO2為氧氣消耗量，mol；Qa為煤氧僅應(yīng)熱，J；為氧氣濃度，mol/m3；A為指前因子，s-1；E為活化能，J/mol；R為氣體常數(shù)，J/(mol·K)；n為反應(yīng)級數(shù)； ΔQ為煤氧化學(xué)反應(yīng)過程中生熱量，J/mol；sV為形狀因子。

1.4 工作面及巷道控制方程

煤礦井下采場工作面和巷道內(nèi)區(qū)域氣體一般都處于湍流狀態(tài)，因此需要使用湍流模型進(jìn)行模擬。選用L-VEL 湍流模型[24]，該模型基于普朗特混合長度理論，根據(jù)局部流速和與最近壁面的距離來計算湍流黏度，穩(wěn)定且計算效率高，適合用于計算內(nèi)部流動。流動、傳熱及傳質(zhì)方程如下：

式中： μT為有效湍流黏度，kg·m2/s；I為壓力在空向上的分布；ρ為流體密度，kg/m3；Cp為空氣比熱容，J/(kg·K)；λf為空氣導(dǎo)熱系數(shù)，W/(m·K)；Di為組分i的擴(kuò)散系數(shù)，m2/s。

1.5 圍巖（煤）傳熱控制方程

圍巖（煤）固相傳熱以熱傳導(dǎo)為主，可用下式表示：

式中： (ρC)s為 圍巖（煤）熱容，J/(m3·K)； λs為導(dǎo)熱系數(shù)，J/(m·K)。

1.6 多場耦合關(guān)系

建立的模型具有3 個方面的耦合關(guān)系。首先，在采空區(qū)內(nèi)部，采空區(qū)能量方程和氣體組分方程相互耦合，并且均與采空區(qū)流動方程所求解的速度場相關(guān)聯(lián)。其次，工作面和巷道氣體傳熱方程和氣體組分方程均與其內(nèi)部速度場相關(guān)聯(lián)。最后，溫度、風(fēng)速、壓力、氣體體積分?jǐn)?shù)在采空區(qū)和工作面區(qū)域的界面上數(shù)值相等，即具有連續(xù)性。

1.7 工作面推進(jìn)過程中計算網(wǎng)格的處理

隨著工作面向前推進(jìn)，采空區(qū)面積不斷增大，計算區(qū)域物理尺寸不斷變化，因此計算網(wǎng)格也應(yīng)隨之調(diào)整。利用COMSOL Multiphysics?軟件的變形幾何模型，通過設(shè)定計算邊界的移動速度，實現(xiàn)網(wǎng)格的自由變形，從而實現(xiàn)對采空區(qū)動態(tài)增長過程的模擬。

1.8 初始條件及邊界條件

通過模型耦合關(guān)系的分析，明確了采空區(qū)、面間煤柱、邊界煤柱、底板巖層、頂板巖層、待回采實體煤、進(jìn)風(fēng)巷、回風(fēng)巷等地點的初始條件和邊界條件，如下：

進(jìn)風(fēng)巷及回風(fēng)巷面間煤柱固壁、工作面固壁、始采線邊界煤柱固壁和終采線邊界煤柱固壁滿足無滑移邊界條件。

2 模型驗證

利用在開灤集團(tuán)崔家寨煤礦E12604 工作面獲得的實測數(shù)據(jù)與建立的數(shù)學(xué)模型計算結(jié)果進(jìn)行對比分析。崔家寨煤礦E12604 工作面推進(jìn)距離570 m，傾向長度135 m，煤層厚度3.0～3.2 m，屬易自燃煤層，最短自然發(fā)火期90 d，采煤方法為綜采傾斜長壁、一次采全高，采高為2.8～3.0 m?，F(xiàn)場數(shù)據(jù)監(jiān)測點距離E12604 工作面開切眼120 m，距回風(fēng)巷底板2.5 m。數(shù)值模擬參數(shù)為：①入口風(fēng)速vin=1.01 m/s；②空 氣 入 口 溫 度Tin=295.85 K；③活化能E=21.6 kJ/mol；④進(jìn)風(fēng)流O2濃度C0=9.375 mol/m3；⑤指前因子A=0.85 s-1；⑥氧化反應(yīng)熱400.11 kJ/mol；⑦圍巖初始溫度295.85 K；⑧煤 巖 密 度1 300 kg/m3；⑨煤 巖 比 熱 容1 003 J/（kg·K）；⑩煤巖導(dǎo)熱系數(shù)λs=0.2 W/（m·K）；?初始碎脹系數(shù)Kp,max=1.5；?壓實碎脹系數(shù)Kp,min=1.15；?碎脹系數(shù)衰減率αx=0.268 m-1；?碎脹系數(shù)衰減率αy=0.0368 m-1；?工作面采高3 m；?工作面推進(jìn)速度7.5 m/d；?煤矸石直徑dp=0.04 m；?工作面長度238 m。崔家寨煤礦E12604 工作面采空區(qū)氣體溫度和O2體積分?jǐn)?shù)的模擬與實測結(jié)果如圖2。

圖2 模擬與實測數(shù)據(jù)對比Fig.2 Comparison of simulated and measured data

從圖2 可以看出：在工作面回采初期采空區(qū)溫度模擬結(jié)果與實測數(shù)據(jù)比較接近，隨著推進(jìn)距離的增加，測點位置的模擬溫度值逐漸高于現(xiàn)場實測值，這可能是由于實際采空區(qū)靠近壁面附近存在鄰近采空區(qū)的漏風(fēng)源，導(dǎo)致實測溫度偏低；O2體積分?jǐn)?shù)的模擬值與實測值的變化規(guī)律一致，工作面回采初期O2體積分?jǐn)?shù)較高，隨著推進(jìn)距離增加，采空區(qū)內(nèi)部O2逐漸被煤氧反應(yīng)所消耗，體積分?jǐn)?shù)不斷降低?？傮w上，可以認(rèn)為建立的模型能夠較好地描述采空區(qū)煤自熱過程隨工作面推進(jìn)的動態(tài)變化。

3 工作面推進(jìn)條件下多物理場的動態(tài)演化

工作面推進(jìn)過程中采空區(qū)空隙率的變化如圖3。

圖3 工作面推進(jìn)過程中采空區(qū)空隙度的動態(tài)變化Fig.3 Dynamic change of porosity in the mining area during the advance of the working face

從圖3 可以看出：采空區(qū)空隙率呈現(xiàn)非均勻分布特征；總體上看，距離工作面和兩巷煤壁較近的區(qū)域，由于液壓支架及巷道煤柱的支撐作用，空隙率較大；距離工作面和兩巷煤壁較遠(yuǎn)的區(qū)域，由于礦山壓力作用明顯，壓實程度較高，空隙率較?。徊煽諈^(qū)空隙率分布規(guī)律具有典型的“O”形圈特征；隨著工作面不斷向前推進(jìn)，采空區(qū)范圍不斷擴(kuò)大，采空區(qū)走向長度增加，“O”形圈尺寸不斷增大，但采空區(qū)空隙率分布規(guī)律始終保持不變；工作面回采時間分別為10 d 和60 d 時，采空區(qū)的空隙率最小值分別為0.18 和0.14，而在工作面液壓支架附近空隙率始終保持在0.33 左右，這表明隨著工作面不斷推進(jìn)，采空區(qū)內(nèi)逐漸壓實。

工作面推進(jìn)過程中采空區(qū)溫度場的動態(tài)變化如圖4。

圖4 工作面推進(jìn)過程中溫度場的動態(tài)變化Fig.4 Dynamic change of temperature field during the advance of working face

從圖4 可以看出：隨工作面推進(jìn)距離的增加，高溫區(qū)域的形狀有較大變化；開采初期，采空區(qū)走向長度相對于傾向長度較短，采空區(qū)整體范圍內(nèi)氧氣體積分?jǐn)?shù)較高，因此升溫區(qū)域及等溫線分布明顯地向回風(fēng)側(cè)延展；在動態(tài)采空區(qū)煤自燃的過程中，高溫區(qū)域隨工作面推進(jìn)向前移動，且與工作面的相對距離保持不變；隨著采空區(qū)走向長度的增加和氧化升溫時間延續(xù)，采空區(qū)固相煤體溫度上升且高溫區(qū)域逐漸擴(kuò)大。然而，由于采空區(qū)固相煤體與底板之間的換熱效應(yīng)，高溫煤體在進(jìn)入采空區(qū)深部窒息區(qū)后有一定的降溫過程。因此，采空區(qū)高溫區(qū)域的傾向?qū)挾妊亻_切眼方向逐漸減小，并且由于熱量的傳遞需要一段時間，存在“熱慣性”現(xiàn)象，而采空區(qū)深部風(fēng)流速度較小，熱量不宜散失，因此在采空區(qū)的進(jìn)風(fēng)側(cè)溫度場形成拖尾現(xiàn)象。

停采后第10 d 的溫度場分布如圖5。

圖5 停止開采后第10 d 溫度場分布Fig.5 Temperature field distribution on the 10th day after mining stopped

從圖5 可以看出：停采后，采空區(qū)內(nèi)部由于缺少氧氣供應(yīng)，煤氧反應(yīng)減弱，溫度逐漸降低；由于工作面附近氧氣供應(yīng)充足，高溫區(qū)域逐漸向工作面靠近，高溫中心距離工作面約100 m，最高溫度64.2 ℃。因此，對于工作面推進(jìn)過程中形成的高溫區(qū)域，如果沒有形成災(zāi)變，那么在采空區(qū)的堆積壓實和圍巖散熱作用下溫度會逐漸降低，但開采停止后高溫區(qū)域會向工作面遷移。

工作面推進(jìn)過程中采空區(qū)氧氣濃度場的動態(tài)變化如圖6。

圖6 工作面推進(jìn)過程中采空區(qū)氧氣濃度場的動態(tài)變化Fig.6 Dynamic changes of the oxygen concentration field in the mining area during the advance of the working face

從圖6 可以看出：隨工作面推進(jìn)距離的增加，采空區(qū)中部及深部區(qū)域固相煤體逐漸被壓實，漏風(fēng)阻力增大，采空區(qū)氧氣濃度呈現(xiàn)不對稱性分布；進(jìn)風(fēng)側(cè)的氧氣濃度分布區(qū)域較回風(fēng)側(cè)更深入，氧氣濃度的等值線分布與工作面的相對距離固定，采空區(qū)氧氣濃度場分布趨于穩(wěn)定。

停止開采后采空區(qū)第10 d 的氧氣濃度分布如圖7。

圖7 停止開采后第10 d 氧氣濃度場分布Fig.7 Oxygen concentration field distribution on the 10th day after the mining stopped

從圖7 可以看出：停止開采后，氧氣分布帶形狀變化不大，仍然呈不對稱分布；進(jìn)風(fēng)側(cè)寬度大于回風(fēng)側(cè)寬度，但總的寬度有所減小，這是由于采空區(qū)內(nèi)部的氧氣逐漸耗盡，而采空區(qū)內(nèi)部壓實，風(fēng)流向采空區(qū)內(nèi)部輸送的氧氣量非常有限；工作面封閉后，進(jìn)風(fēng)側(cè)較回風(fēng)側(cè)積存的氧氣多，煤氧化產(chǎn)生的熱量不能及時帶到回風(fēng)側(cè)，導(dǎo)致進(jìn)風(fēng)側(cè)溫度逐漸升高。

4 敏感性分析

當(dāng)圍巖溫度為20 ℃時，不同通風(fēng)溫度下采空區(qū)內(nèi)最高溫度如圖8。

圖8 不同通風(fēng)溫度下采空區(qū)內(nèi)最高溫度Fig.8 Maximum temperature in the extraction zone at different ventilation temperature conditions

從圖8 可以看出：首先在不同的通風(fēng)溫度下，隨著時間推移，采空區(qū)內(nèi)最高溫度先迅速增加，后緩慢下降，這是因為在采空區(qū)推進(jìn)的前期，由于采空區(qū)走向長度較短，風(fēng)流阻力小，風(fēng)速較大，采空區(qū)內(nèi)氧氣供應(yīng)條件較好，煤氧反應(yīng)劇烈，煤氧反應(yīng)產(chǎn)熱與散熱相比占主導(dǎo)地位；隨著采空區(qū)不斷向前推進(jìn)，采空區(qū)走向長度逐漸增加，采空區(qū)內(nèi)風(fēng)流阻力逐漸變大，風(fēng)速減小，采空區(qū)內(nèi)氧氣供應(yīng)條件變差，煤氧反應(yīng)相對較弱，煤氧反應(yīng)熱小于煤巖散熱，因此溫度略有下降；其次，通風(fēng)溫度越高，采空區(qū)內(nèi)最高溫度越大，但最高溫度的增幅要小于通風(fēng)溫度的增幅，且工作面推進(jìn)的前期最高溫度非常接近；通風(fēng)溫度為10、30 ℃時，第30 d 的最高溫度分別為55.2、59.4 ℃，通風(fēng)溫度提升了10 ℃，但是最高溫度的增幅小于5℃，這是由于空氣在進(jìn)入采空區(qū)前會與巷道壁面進(jìn)行換熱，進(jìn)入采空區(qū)后溫度降低，并且由于煤矸石的熱容量遠(yuǎn)大于空氣的熱容量。

不同圍巖溫度下采空區(qū)內(nèi)最高溫度的演化如圖9。

圖9 不同圍巖溫度下采空區(qū)內(nèi)最高溫度演化Fig.9 Evolution of maximum temperature in mining area under different surrounding rock temperature conditions

從圖9 可以看出：在開采的初期，采空區(qū)內(nèi)最高溫度主要受圍巖溫度的影響，兩者較為接近。隨著開采向前推進(jìn)，不同圍巖溫度下采空區(qū)內(nèi)最高溫度逐漸接近，這是由于遺煤氧化產(chǎn)生的熱量對采空區(qū)溫度的影響逐漸占主導(dǎo)地位；另外，采空區(qū)內(nèi)最高溫度先迅速上升后緩慢下降的趨勢與圖8 類似。

不同推進(jìn)速度下采空區(qū)內(nèi)最高溫度的演化如圖10。

圖10 不同推進(jìn)速度下采空區(qū)內(nèi)最高溫度Fig.10 Maximum temperature in the extraction zone at different advance speed conditions

從圖10 可以看出：工作面推進(jìn)速度越快，采空區(qū)內(nèi)最高溫度越低。例如推進(jìn)速度為2 m/d 時，采空區(qū)最高溫升接近194 ℃，而推進(jìn)速度為8 m/d 時，采空區(qū)最高溫不超過47 ℃；這是由于工作面推進(jìn)速度越快，采空區(qū)遺煤越快進(jìn)入窒息帶，煤氧反應(yīng)受到抑制，釋放的熱量相對較少，采空區(qū)環(huán)境溫度上升較慢。

5 結(jié) 語

1）由于“熱慣性”的存在，采空區(qū)的進(jìn)風(fēng)側(cè)溫度場形成高溫區(qū)域拖尾現(xiàn)象。工作面推進(jìn)過程中形成的高溫區(qū)域在采空區(qū)的堆積壓實和圍巖散熱作用下溫度會逐漸降低，開采停止后高溫區(qū)域會向工作面遷移。

2）相同的圍巖溫度下，不同通風(fēng)溫度對采空區(qū)最高溫度的影響，在開采的初期幾乎無影響，而在開采的后期，通風(fēng)溫度越高，采空區(qū)內(nèi)最高溫度越高。

3）工作面推進(jìn)過程中，采空區(qū)內(nèi)最高溫度先迅速增加，然后緩慢下降。當(dāng)開采停止后，采空區(qū)內(nèi)繼續(xù)升溫但高溫區(qū)域向工作面遷移。

4）工作面推進(jìn)速度越快，采空區(qū)遺煤越快進(jìn)入窒息帶，采空區(qū)內(nèi)最高溫度越低。