杜 斌，趙 磊

（1.中煤科工集團沈陽研究院有限公司，遼寧 撫順 113122；2.煤礦安全技術(shù)國家重點實驗室，遼寧 撫順 113122；3.平頂山天安煤業(yè)股份有限公司，河南 平頂山 467099）

采空區(qū)遺煤自然發(fā)火，威脅井下安全生產(chǎn)，還可能引發(fā)采空區(qū)瓦斯爆炸，其條件之一是采空區(qū)漏風持續(xù)地為遺煤供氧。因此國內(nèi)外許多學者通過現(xiàn)場實測、數(shù)值模擬等方法對采空區(qū)內(nèi)氣體流動的速度場、壓力場、氧濃度場開展了研究。

采空區(qū)內(nèi)氣體流動規(guī)律在宏觀尺度服從N-S方程，氣體流動驅(qū)動力主要是工作面上下隅角的氣壓差。然而，煤與氧的物理吸附、化學吸附、化學反應過程主要在煤基質(zhì)的表面（氣-固界面）進行，包括煤塊與采空區(qū)氣體的宏觀氣-固界面和煤塊中孔隙氣體與基質(zhì)的介觀氣-固界面。因此除了需要對采空區(qū)宏觀氣體流動規(guī)律進行研究，還需要對煤孔隙中的介觀氣體流動規(guī)律進行研究。然而由于研究尺度不同，煤孔隙中的氣體流動將不再像宏觀氣體流動那樣服從N-S 方程，而是服從達西定律。梁運濤等[1-2]將發(fā)生自然發(fā)火的遺煤按照從小到大的層次劃分為5 個尺度：反應粒子、顆粒內(nèi)孔、單煤顆粒、顆粒聚團、堆積媒體，提出使用關(guān)聯(lián)型多尺度方法，對自然發(fā)火過程進行跨尺度研究。但是目前尚未有學者對2 種尺度下多孔介質(zhì)內(nèi)的氣體流動進行關(guān)聯(lián)型分析。另外，煤孔隙中的殘余水分占據(jù)了氣體流動通道，甚至存在“水鎖”效應，阻礙了氣體流動，同時，氣體流動產(chǎn)生的蒸發(fā)效應，加快了殘余水分的蒸發(fā)，減小了氣體流動通道內(nèi)阻礙物質(zhì)的體積，增加了有效滲透率，利于氣體流動，因此需要考慮殘余水分與氣體流動間的耦合影響，即氣液兩相流之間的耦合影響；同時水分蒸發(fā)吸熱使煤體溫度降低，又反作用于氣體流動和水分的蒸發(fā)速率，使氣體分子的熱運動和蒸發(fā)速率降低。因此，對介觀尺度下煤中熱-流-濕多場耦合機制進行研究。

1 采空區(qū)多尺度氣體流動規(guī)律

在礦井采場空間范圍內(nèi)共存在5 個層次的研究尺度，其中邊界尺度為反應粒子、單煤顆粒、堆積煤場，介于邊界尺度之間的為顆粒內(nèi)孔（介尺度I）和顆粒聚團（介尺度II）。從宏觀尺度入手，按照尺度逐漸減小的順序，對采空區(qū)松散煤體及之間的空隙、遺煤基質(zhì)及孔隙裂隙內(nèi)的氣體流動情況展開研究。

1.1 采空區(qū)宏觀尺度氣體流動規(guī)律

選取長度300 m、寬度180 m 的采空區(qū)堆積煤體范圍作為研究對象，其中0～40 m 為自然堆積區(qū)、40～100 m 為載荷影響區(qū)、100～300 m 為壓實穩(wěn)定區(qū)[3]。進風巷和回風巷寬度為5 m。將采空區(qū)內(nèi)巖石垮落堆積狀態(tài)及其空隙視為多孔介質(zhì)，相關(guān)模擬參數(shù)為：①空氣摩爾質(zhì)量：29 g/mol；②采空區(qū)空隙率：0.29；③動力黏度：2.01×10-5Pa·s；④氣體擴散系數(shù)：2.88×10-2；⑤采空區(qū)區(qū)域：300 m×180 m；⑥自然堆積區(qū)滲透率：5.38×10-6m2；⑦載荷影響區(qū)滲透率：2.6×10-6m2；⑧壓實穩(wěn)定區(qū)滲透率：1.3×10-6m2；⑨CO2生成速率：1.026×10-3mol/（m3·s）；⑩CO 生成速率：3.4×10-4mol/（m3·s）。

根據(jù)實際情況，在COMSOL Multiphysics 數(shù)值模擬軟件中設定采空區(qū)入口風速為0.9 m/s，分析采空區(qū)內(nèi)氣體流速分布情況。采空區(qū)氣體流場示意圖如圖1。

圖1 采空區(qū)氣體流場示意圖Fig.1 Distribution of gas flow field in goaf

從圖1 可以發(fā)現(xiàn)：在流線密集處，采空區(qū)漏風風速較大；在流線稀疏處，采空區(qū)漏風風速較??；流線密集處主要集中在采空區(qū)的進風側(cè)、回風側(cè)以及采空區(qū)的淺部。為了更詳細具體地分析采空區(qū)內(nèi)氣體流場的分布情況，進而根據(jù)采空區(qū)漏風風速劃分采空區(qū)自燃“三帶”的分布范圍，在采空區(qū)中部沿工作面走向方向設置了1 條二維截線，借此分析采空區(qū)漏風風速在二維截線上各點的分布情況；二維截線的起點坐標為（0，80），終點坐標為（260，80）；二維截線的位置如圖1 中紅色線段。

根據(jù)國內(nèi)外學者對采空區(qū)漏風的研究，采空區(qū)“三帶”的范圍可以根據(jù)采空區(qū)漏風流速劃分[5]。采空區(qū)漏風風速沿二維截線的分布情況如圖2。

圖2 沿采空區(qū)二維截線（y=80 m）風流速度分布曲線Fig.2 Curve of gas velocity in goaf along two-dimensional section line（y＝80 m）

從圖2 中可以看出：在0～40 m 之間風速較大，在0～40 m 范圍內(nèi)的遺煤由于散熱速率大于產(chǎn)熱速率，煤體氧化產(chǎn)生的熱量大部分都被漏風攜帶走，不具備良好的蓄熱條件，因此將這一區(qū)間劃分為“散熱帶”；而在135 m 之后，漏風風速較小，煤體氧化產(chǎn)生的熱量大部分留存在松散的煤體內(nèi)，散熱速率小于產(chǎn)熱速率，具有良好的蓄熱條件，但卻由于漏風風速偏小，不能及時地為煤提供足量的氧氣，造成煤的氧化反應進程受到抑制，因此將這一區(qū)間劃分為“窒息帶”；在40～135 m 之間的采空區(qū)遺煤由于漏風風速適中，既能提供較充足氧氣供給煤的氧化反應進程，又具有較好的儲熱條件，因此將這一區(qū)間劃分為“氧化升溫帶”。

1.2 采空區(qū)介觀尺度氣體流動規(guī)律

采空區(qū)某處遺煤及附近區(qū)域幾何結(jié)構(gòu)示意圖如圖3。

圖3 采空區(qū)某處遺煤及附近區(qū)域幾何結(jié)構(gòu)示意圖Fig.3 Schematic diagram of geometric structure of residual coal and nearby area in goaf

為了進一步研究介觀尺度下顆粒聚團和顆粒內(nèi)孔內(nèi)部氣體流動狀態(tài)，選取采空區(qū)氧化升溫帶內(nèi)某處（x=40 m，y=80 m）遺煤及附近區(qū)域作為研究對象。研究區(qū)域的長度為80 cm，高度為50 cm，假設其中的煤塊為橢圓形，煤塊長軸為24 cm，短軸為8 cm，與采空區(qū)底板的夾角為40°。實驗測得孔隙率為9.5%，滲透率為1×10-16m2[6]。根據(jù)關(guān)聯(lián)型多尺度研究方法，煤塊作為顆粒聚團尺度，其入口風速由上一級尺度所決定，顆粒聚團處的風速為0.04 m/s。

模擬后得到顆粒聚團及附近區(qū)域氣體流動的速度場和壓力場分布情況如圖4 和圖5。

圖4 采空區(qū)某顆粒聚團附近氣體流動速度場分布示意圖Fig.4 Distribution diagram of gas velocity field near a particle cluster in goaf

圖5 采空區(qū)某顆粒聚團附近氣體壓力場分布示意圖Fig.5 Distribution diagram of gas pressure field near a particle cluster in goaf

1.2.1 介觀尺度垂直方向氣體流動規(guī)律

為了進一步研究顆粒內(nèi)孔尺度下，孔隙內(nèi)部的氣體在垂直方向的流動情況，在圖3 中設置1 條垂直于底板且縱貫顆粒聚團的二維截線，函數(shù)表達式為x=40 cm（6 cm＜y＜14.5 cm），顆粒聚團內(nèi)部垂直方向二維截線位置圖如圖6，垂直方向二維截線如圖6中黃色線段。

圖6 顆粒聚團內(nèi)部垂直方向二維截線位置圖Fig.6 Location diagram of two-dimensional vertical section line in a particle cluster

顆粒內(nèi)孔氣體流動速度沿垂直二維截線的分布曲線如圖7。從圖7 可以發(fā)現(xiàn)：越靠近顆粒聚團的表面，顆粒內(nèi)孔中氣體流動的速度越快。

圖7 顆粒內(nèi)孔氣體流動速度沿垂直二維截線的分布曲線Fig.7 Distribution curve of gas velocity in pores along vertical two-dimensional section line

1.2.2 介觀尺度水平方向氣體流動規(guī)律

為了進一步研究顆粒內(nèi)孔尺度下孔隙內(nèi)部的氣體在水平方向的流動情況，在圖3 中設置1 條平行于底板且貫穿顆粒聚團的二維截線，函數(shù)表達式為y=7 cm（32 cm＜x＜42 cm），顆粒聚團內(nèi)部水平方向二維截線位置示意圖如圖8，水平方向二維截線如圖8 中黃色線段。

圖8 顆粒聚團內(nèi)部水平方向二維截線位置示意圖Fig.8 Location diagram of two-dimensional horizontal section line in a particle cluster

顆粒內(nèi)孔氣體流動速度沿水平二維截線的分布曲線如圖9。從圖9 可以發(fā)現(xiàn)：越靠近顆粒聚團的迎風側(cè)，顆粒內(nèi)孔的氣體流動速度越快。通過對介尺度層面顆粒聚團內(nèi)部孔隙氣體流動速度進行分析可以發(fā)現(xiàn)：顆粒內(nèi)孔氣體流動速度的單位大約在納米級別，這和顆粒內(nèi)孔的孔徑是在一個尺度層面的，證明了模擬結(jié)果正確性。另外可以發(fā)現(xiàn)：越靠近顆粒聚團表面，越靠近顆粒聚團的迎風側(cè)，顆?？紫吨械臍怏w流動速度越大。

圖9 顆粒內(nèi)孔氣體流動速度沿水平二維截線的分布曲線Fig.9 Distribution curve of gas velocity in pores along horizontal two-dimensional section line

2 采空區(qū)氣體流場對遺煤濕分的熱蒸發(fā)作用分析

從圖7 和圖9 中可以看出：顆粒內(nèi)孔中氣體流動的速度大部分在10～20 pm/s 之間，且越靠近氣-固交界面，氣體流速越快。煤作為一種孔隙-裂隙雙孔隙率多孔介質(zhì)，孔隙-裂隙不僅是氣體發(fā)生流動的主要結(jié)構(gòu)，水分也主要在這些孔隙結(jié)構(gòu)中進行儲存和輸運。滕騰[6]認為煤體裂隙中的殘余水分在表面張力的作用下會鑲嵌在迂曲基質(zhì)塊表面的溝塹和凹槽里，或者黏附在基質(zhì)表面，也即裂隙表面；殘余水分的存在不僅降低了基質(zhì)表面的裸露面積，同時還減小了裂隙的有效開度；并且認為裂隙網(wǎng)絡中的非連續(xù)殘余水分會隨著氣體流動發(fā)生蒸發(fā)遷移。賀瓊瓊[7]通過采用NMR 方法分析褐煤中水分的賦存情況，發(fā)現(xiàn)勝利褐煤中的水分主要集中在0.003～0.3 μm 孔隙范圍內(nèi)，昭通褐煤中的水分主要集中在0.01～1 μm 孔隙范圍內(nèi)，褐煤中1 個有效親水點位形成的水簇尺寸大約由6～7 個水分子組成；然后使用N2吸附差值法研究了介孔中水分的分布情況，發(fā)現(xiàn)大部分水分集中在26～100 nm 的孔隙中；大孔范圍和裂隙并非儲水結(jié)構(gòu)，而是僅作為脫水過程中的傳質(zhì)通路。介觀尺度下的氣體也主要在這些傳質(zhì)通路中流動。介觀尺度下煤中水分的儲水及輸運模型如圖10。

圖10 基質(zhì)內(nèi)孔與裂隙結(jié)構(gòu)的“樹杈”型儲水及輸運模型[2]Fig.10“Tree branch”type water storage and transport model of pore and fracture structure in matrix

當發(fā)生氧化蓄熱升溫現(xiàn)象時，顆粒內(nèi)孔中的水分子會在熱蒸發(fā)效應下發(fā)生相變成為氣相水分子，之后通過擴散作用混合進入干燥氣體中，并沿著輸運通路進行運移。

3 熱蒸發(fā)過程物理模型

采空區(qū)氣體流場與遺煤之間的熱濕傳遞過程主要涉及3 種物理過程：流體流動、傳熱、液體的相變以及與空氣之間的輸運。在自由流動域主要發(fā)生的是空氣的流動、水分與空氣之間的輸運等物理效應；在多孔介質(zhì)域主要發(fā)生的是水分的流動輸運、相變吸熱等物理效應。

3.1 液 相

在多孔介質(zhì)域存在氣相和液相2 種不同的流動相，且滿足以下約束：

式中：Sg為氣相的飽和度，對應于多孔介質(zhì)中的濕空氣；Sl為液相的飽和度，對應多孔介質(zhì)中的水分。

多孔介質(zhì)孔隙中氣相為不飽和狀態(tài)流動，所以使用Brinkman 方程計算多孔介質(zhì)中濕空氣的流場ug和壓力分布pg。與濕空氣速度相比，液相速度較小，因此根據(jù)氣相壓力梯度定義的Darcy 定律[8]可用于計算水流速度ul：

式中：κl為液相在煤體中的有效滲透率，m2；μl為液相的動力黏度系數(shù)，（N·s）/m2；▽pg為氣相壓力分布梯度，Pa/m。

多孔介質(zhì)中的液相為不飽和流動，其有效滲透率與煤體絕對滲透率、相對滲透率、飽和度有關(guān)，因此式（2）可以進一步表示為：

式中：κ 為煤體絕對滲透率，m2；κrl為液相的相對滲透率；μw為水的動力黏度系數(shù)，（N·s）/m2；ε 為孔隙率。

液相的飽和度與水的濃度呈線性關(guān)系：

式中：cw為液相在煤體中的濃度，mol/m3；Mw為水的摩爾質(zhì)量，g/mol；ρw為水的密度，kg/m3。

相對滲透率通常根據(jù)經(jīng)驗或者實驗來確定，其形式在很大程度上取決于多孔材料屬性和液體本身。根據(jù)參考文獻［9］，將相對滲透率始終定義為正：

式中：κrg為氣相的相對滲透率；eps 為浮點相對精度。

式中：Sli為殘留在多孔介質(zhì)中的液相的不可約飽和度。

3.2 水蒸氣相

多孔介質(zhì)中水蒸氣的傳遞涉及2 個過程：第1個傳遞機制是根據(jù)Brinkman 方程求出的流體速度場，將其值賦給蒸汽相；第2 個傳遞機制是水蒸氣及其氣相干燥空氣的二元擴散產(chǎn)生的速度分量。這2個組分之間的有效擴散系數(shù)Deff通常可以采用Millington & Quirk 方程進行描述：

式中：Dva為蒸汽-空氣擴散率，取2.6×10-5m2/s。

將這2 種擴散效應得出的速度場應用于水蒸氣傳遞方程，可以得到：

式中：uv為濕空氣的速度場；ug為水蒸氣的速度場；Ma為干燥空氣的摩爾質(zhì)量，g/mol；Mma為濕空氣的摩爾質(zhì)量，g/mol；ρma為濕空氣的密度，kg/m3。

將濕空氣視為理想氣體，濕空氣的密度和壓強之間具有以下關(guān)系：

式中：ρg為氣相的密度，kg/m3；pg為氣相的壓強，Pa；R 為理想氣體常數(shù)，8.314 J/（mol·K）；Tg為氣相的溫度，K；Mg為氣相的摩爾質(zhì)量，g/mol。

3.3 相 變

為了計算蒸發(fā)到空氣中的水量，并考慮到液態(tài)水的減少和濕空氣中水蒸氣濃度的增加，定義了以下計算式：

式中：mevap為蒸發(fā)量，mol/（m3·s）；K 為蒸發(fā)率，1/s；aw為蒸汽系數(shù)；csat為飽和情況下的蒸汽濃度，mol/m3；c 為當前蒸汽濃度，mol/m3。

將水蒸氣視為理想氣體，可得其飽和水蒸氣濃度[10]為：

式中：csat為飽和水蒸氣濃度，mol/m3；psat（T）為飽和水蒸氣分壓，Pa；T 為溫度，K。

不同溫度時所對應的飽和水蒸氣分壓近似為[7]：

對飽和蒸氣分壓的數(shù)學表達式進行數(shù)學分析可以發(fā)現(xiàn)其變化趨勢主要取決于指數(shù)部分。將指數(shù)部分做如下數(shù)學變換：

當溫度在［273.15 K，373.15 K］之間變化時，飽和蒸氣分壓psat（T）的變化區(qū)間為［610.7 Pa，100 kPa］。

4 采空區(qū)氣體流場對遺煤的熱蒸發(fā)效應數(shù)值模擬

為了研究采空區(qū)內(nèi)氣體流場和濕度場之間的相互關(guān)系，使用多物理場模擬軟件，首先對采空區(qū)氣體流動狀態(tài)開展穩(wěn)態(tài)研究，然后在穩(wěn)態(tài)的基礎上，研究熱濕傳遞過程的瞬態(tài)變化。采空區(qū)內(nèi)氣體流動穩(wěn)態(tài)結(jié)果如圖4 和圖5，分別表示采空區(qū)內(nèi)氣體流動的速度和壓力分布情況。瞬態(tài)模擬的時間跨度為［0 h，24 h］，研究在此時間區(qū)間內(nèi)受氣體流動和蒸發(fā)吸熱影響下，自由流動域和多孔介質(zhì)域內(nèi)的水蒸氣濃度和濕度的變化情況。

4.1 采空區(qū)氣體流場對遺煤中水分的影響

4.1.1 多孔介質(zhì)域中水分含量的變化情況

多孔介質(zhì)域水分含量隨時間的變化情況如圖11。多孔介質(zhì)域中迎風側(cè)某點水分含量的變化曲線如圖12。多孔介質(zhì)域內(nèi)含水量隨時間變化情況如圖13。

圖11 多孔介質(zhì)域水分含量隨時間的變化情況Fig.11 Variation of water content with time in porous media domain

圖12 多孔介質(zhì)域中迎風側(cè)某點水分含量的變化曲線Fig.12 Variation curve of moisture content on the windward side in porous media domain

從圖11 可以發(fā)現(xiàn)：煤中水分的蒸發(fā)首先從煤塊外表面開始，之后逐步影響煤塊內(nèi)部；而且迎風側(cè)比背風側(cè)水分含量小。

從圖12 中可以看到：觀測點處的水分含量變化趨勢為“下降-上升-下降”，且2 個下降階段的曲線斜率大致相同，所以認為這2 個階段中水分含量下降主要是由于熱蒸發(fā)效應導致，中間階段水分含量上升主要是受多孔介質(zhì)內(nèi)流場的影響，將迎風側(cè)的水分向觀測點處堆積遷移所導致的。

從圖13 可以看出：通過對多孔介質(zhì)域的水分含量分別在0 h 和24 h 進行積分運算，可以得到在［0 h，24 h］的時間區(qū)間內(nèi)，水分含量由于蒸發(fā)效應減少了0.046 mol/m。這與劉劍等[11]之前使用熱重分析得到的結(jié)果是相同的。即在煤常溫氧化的早期階段，煤重量的減少主要是由于煤中水分的蒸發(fā)所導致的。

4.1.2 自由流動域水蒸氣濃度的變化情況

采空區(qū)內(nèi)水蒸氣濃度隨時間變化情況如圖14。

圖14 采空區(qū)內(nèi)水蒸氣濃度隨時間變化情況Fig.14 Variation of water vapor concentration in goaf with time

從圖14 可以看出：在0.5 h 時自由流動域中水蒸氣濃度最高處位于煤塊背風側(cè)風流尾部，在6、12、18 h 時自由流動域中水蒸氣濃度最高處位于煤塊的外表面，說明蒸發(fā)現(xiàn)象早期主要發(fā)生在煤塊表層，導致了自由流動域氣-固交界面處水蒸氣濃度較低；后期水分才逐漸從煤塊內(nèi)部遷移至外表面處蒸發(fā)，導致了外表面的水蒸氣濃度較高，說明在早期的表層蒸發(fā)和后期的水分遷移之間存在1 個啟動過程。

4.2 熱蒸發(fā)效應對采空區(qū)遺煤溫度的影響

在多孔介質(zhì)域，由于多孔介質(zhì)域內(nèi)的水分蒸發(fā)吸熱，會使多孔介質(zhì)域內(nèi)的溫度分布發(fā)生改變，多孔介質(zhì)域中某點的蒸發(fā)吸熱量Q 與蒸發(fā)量有關(guān)：

式中；Hevap為蒸發(fā)潛熱，J/mol；mevap為蒸發(fā)量，mol/（m3·s）。而多孔介質(zhì)域內(nèi)的溫度服從固體熱傳導規(guī)律，符合傅立葉傳熱方程。為了計算多孔介質(zhì)域內(nèi)的溫度分布情況，需要將蒸發(fā)熱作為源項插入到傅立葉傳熱方程中。需要注意的是，采空區(qū)氧化升溫帶內(nèi)的遺煤在蓄熱條件良好的情況下會將氧化蓄熱積存在煤體內(nèi)部，導致煤體溫度緩慢升高，但這一過程不在研究范圍內(nèi)，此處僅考慮由于蒸發(fā)吸熱導致的煤體溫度變化情況。

多孔介質(zhì)域蒸發(fā)吸熱對煤體溫度變化的影響如圖15。

圖15 多孔介質(zhì)域蒸發(fā)吸熱對煤體溫度變化的影響Fig.15 Effect of evaporation and heat absorption in porous media on temperature of coal

模擬設定的環(huán)境溫度為293.15 K，在2 h 時，煤塊內(nèi)部溫度高于環(huán)境溫度，在煤塊的上部和下部溫度較高。在6～12 h 時，煤塊內(nèi)部溫度依然高于環(huán)境溫度，但是從表層開始溫度逐漸降低，煤塊內(nèi)部高于環(huán)境溫度的區(qū)域越來越小。在18、24 h 時，煤塊內(nèi)部溫度由于蒸發(fā)吸熱，會出現(xiàn)溫度低于環(huán)境溫度的情況。

5 結(jié) 語

1）利用關(guān)聯(lián)型多尺度研究方法，基于采空區(qū)氣體流動在宏觀尺度服從的N-S 方程和在介觀尺度服從的達西定律，使用多物理場模擬軟件研究了采空區(qū)流場速度和壓力的跨尺度分布情況。

2）基于多孔介質(zhì)模型總結(jié)了采空區(qū)遺煤在熱蒸發(fā)效應下所涉及的氣-液兩相非飽和流動、相變、蒸發(fā)吸熱等物理過程，構(gòu)建了關(guān)于流體場、濕度場、溫度場的多物理場耦合模型。并使用多物理場模擬軟件研究了采空區(qū)多孔介質(zhì)域孔隙內(nèi)水分含量的變化情況以及自由流動域內(nèi)水蒸氣濃度的變化情況。

3）初步探討了多孔介質(zhì)域由于蒸發(fā)吸熱導致的煤體溫度分布變化情況。