動態推進U 型采場煤自燃多場耦合數值研究

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

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

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

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

1 模型構建

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

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

1.1 基本假設

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

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

1.2 采空區孔隙率及滲透率的動態演化

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

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

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

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

垮落采空區的高度H和孔隙率 εb的 表達式分別為：

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

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

式中：dp為采空區多孔介質的平均顆粒直徑，m。

1.3 采空區流動及傳熱方程

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

氧氣消耗量與氧氣濃度和煤自燃參數有關，煤氧反應熱與氧氣消耗量和反應熱有關，分別計算如下：

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

1.4 工作面及巷道控制方程

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

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

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

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

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

1.6 多場耦合關系

建立的模型具有3 個方面的耦合關系。首先，在采空區內部，采空區能量方程和氣體組分方程相互耦合，并且均與采空區流動方程所求解的速度場相關聯。其次，工作面和巷道氣體傳熱方程和氣體組分方程均與其內部速度場相關聯。最后，溫度、風速、壓力、氣體體積分數在采空區和工作面區域的界面上數值相等，即具有連續性。

1.7 工作面推進過程中計算網格的處理

隨著工作面向前推進，采空區面積不斷增大，計算區域物理尺寸不斷變化，因此計算網格也應隨之調整。利用COMSOL Multiphysics?軟件的變形幾何模型，通過設定計算邊界的移動速度，實現網格的自由變形，從而實現對采空區動態增長過程的模擬。

1.8 初始條件及邊界條件

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

進風巷及回風巷面間煤柱固壁、工作面固壁、始采線邊界煤柱固壁和終采線邊界煤柱固壁滿足無滑移邊界條件。

2 模型驗證

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

圖2 模擬與實測數據對比Fig.2 Comparison of simulated and measured data

從圖2 可以看出：在工作面回采初期采空區溫度模擬結果與實測數據比較接近，隨著推進距離的增加，測點位置的模擬溫度值逐漸高于現場實測值，這可能是由于實際采空區靠近壁面附近存在鄰近采空區的漏風源，導致實測溫度偏低；O2體積分數的模擬值與實測值的變化規律一致，工作面回采初期O2體積分數較高，隨著推進距離增加，采空區內部O2逐漸被煤氧反應所消耗，體積分數不斷降低。總體上，可以認為建立的模型能夠較好地描述采空區煤自熱過程隨工作面推進的動態變化。

3 工作面推進條件下多物理場的動態演化

工作面推進過程中采空區空隙率的變化如圖3。

圖3 工作面推進過程中采空區空隙度的動態變化Fig.3 Dynamic change of porosity in the mining area during the advance of the working face

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

工作面推進過程中采空區溫度場的動態變化如圖4。

圖4 工作面推進過程中溫度場的動態變化Fig.4 Dynamic change of temperature field during the advance of working face

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

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

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

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

工作面推進過程中采空區氧氣濃度場的動態變化如圖6。

圖6 工作面推進過程中采空區氧氣濃度場的動態變化Fig.6 Dynamic changes of the oxygen concentration field in the mining area during the advance of the working face

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

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

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

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

4 敏感性分析

當圍巖溫度為20 ℃時，不同通風溫度下采空區內最高溫度如圖8。

圖8 不同通風溫度下采空區內最高溫度Fig.8 Maximum temperature in the extraction zone at different ventilation temperature conditions

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

不同圍巖溫度下采空區內最高溫度的演化如圖9。

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

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

不同推進速度下采空區內最高溫度的演化如圖10。

圖10 不同推進速度下采空區內最高溫度Fig.10 Maximum temperature in the extraction zone at different advance speed conditions

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

5 結 語

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

2）相同的圍巖溫度下，不同通風溫度對采空區最高溫度的影響，在開采的初期幾乎無影響，而在開采的后期，通風溫度越高，采空區內最高溫度越高。

3）工作面推進過程中，采空區內最高溫度先迅速增加，然后緩慢下降。當開采停止后，采空區內繼續升溫但高溫區域向工作面遷移。

4）工作面推進速度越快，采空區遺煤越快進入窒息帶，采空區內最高溫度越低。