賈男

（1.中煤科工集團沈陽研究院有限公司，遼寧 撫順 113122；2.煤礦安全技術國家重點實驗室，遼寧 撫順 113122）

0 引言

高效開發(fā)非常規(guī)天然氣是增強我國能源保障、落實綠色發(fā)展戰(zhàn)略、助力“雙碳”目標的重要途徑[1]。煤系氣是指蘊藏于煤系中的非常規(guī)天然氣，其中煤層瓦斯是主要組成部分[2]。煤層瓦斯高效抽采是煤礦安全生產(chǎn)的重要前提[3]，然而目前缺乏對煤層瓦斯運移機制的充分認識。

目前，學者們在研究煤體瓦斯微觀運移特性時，多以等效孔隙網(wǎng)絡模型為切入點，對微觀孔隙結構進行滲流過程研究[4]。劉學鋒[5]綜合考慮浸潤性、流體動力學黏度等因素，基于雙相格子Boltzmann 模型對灰?guī)r數(shù)字巖心油水驅(qū)替過程進行模擬。崔冠哲等[6]基于砂巖孔隙空間結構網(wǎng)絡模型，使用格子波爾茲曼法對各孔徑范圍的孔隙結構滲流規(guī)律進行研究，發(fā)現(xiàn)微孔數(shù)量對煤巖體連通性的影響較大。王剛等[7]基于CT 三維重建模型重構了6 種煤樣的孔隙網(wǎng)絡結構模型，探究非達西滲流中速度、壓力、非達西系數(shù)等參數(shù)對滲流的影響。白若男[8]模擬研究了微觀尺度下煤巖體孔隙結構單相、多向低壓水滲流過程。雖然諸多學者探討了煤微觀孔隙瓦斯運移特性[9]，但忽略了瓦斯解吸-擴散過程。

本文以原生結構焦煤孔隙系統(tǒng)為研究對象，綜合運用壓汞測試、納米級工業(yè)CT 掃描和數(shù)值仿真，精準重構煤微觀孔隙空間結構，并對其進行定量表征，研究瓦斯解吸-擴散-滲流的演化過程，初步探討煤孔隙空間結構對瓦斯運移的影響，旨在為煤層瓦斯運移理論提供新的見解。

1 樣品及實驗

1.1 煤樣處理

本文實驗樣品取自華晉焦煤有限責任公司沙曲一號煤礦四采區(qū)5 號煤和6 號煤，均屬焦煤。其中，5 號煤賦存于山西組下部，煤層平均厚度為3.30 m，屬穩(wěn)定大部可采煤層，瓦斯含量為11.16 m3/t，頂板為泥巖、砂質(zhì)泥巖及極少量中-細粒砂巖，底板為粉砂巖、砂質(zhì)泥巖、泥巖。6 號煤賦存于太原組上部，上距5 號煤平均距離為16.53 m，煤層平均厚度為0.55 m，屬不穩(wěn)定的局部可采煤層，瓦斯含量為10.05 m3/t，頂板為石灰?guī)r，底板為泥巖或粉砂巖。原始煤樣采集自工作面新鮮暴露的煤壁，每個煤層分別采集3 塊10 cm×10 cm×10 cm 塊狀樣品，保鮮膜包裹后運至實驗室。塊狀煤樣經(jīng)破碎后打磨加工成5 mm×5 mm×10 mm 的長方體，用于先后完成X 射線CT（X-ray μCT）掃描實驗和壓汞實驗，實驗樣品如圖1 所示。

圖1 實驗樣品Fig.1 Experimental samples

1.2 實驗方法

X-ray μCT 掃描實驗用于獲取并重構煤體微觀孔隙結構，進而開展瓦斯?jié)B流-擴散數(shù)值模擬。實驗采用Phoenix v|tome|xs 型多功能高分辨率CT 掃描儀，該設備配備有高功率納米焦點X 射線管，空間對比能力和圖像對比分辨率較強，如圖2（a）所示。實驗過程中選用的工作電壓為130 kV，工作電流為70 mA，分辨率為0.5 μm，曝光時間為500 ms，采用局部掃描方式。煤樣經(jīng)過X-ray μCT 掃描后，進行壓汞測試。煤樣經(jīng)過二次干燥和脫氣后，完成壓汞測試，從而獲取煤樣孔隙率，用于煤微觀孔隙精準重構。實驗儀器為AutoPore IV 9500 型壓汞儀，如圖2（b）所示。分析孔徑范圍為5 nm～800 μm，最大進汞壓力為228 MPa。實驗前樣品經(jīng)過12 h 烘干，并在膨脹儀中抽真空6 h。

圖2 實驗儀器Fig.2 Experimental instruments

1.3 CT 數(shù)據(jù)處理

基于X-ray μCT 掃描獲得的原始CT 數(shù)據(jù)存在像素噪聲，影響后續(xù)孔隙重構的準確性，有必要進行降噪處理[10]。常見的消除噪聲方法有高斯濾波、均值濾波和中值濾波[11]。高斯濾波主要用于處理符合正態(tài)分布的噪聲，具有一定局限性[12]。均值濾波能擴大目標像素點數(shù)，以8 個像素點的灰度平均值作為目標點的灰度值，誤差較大[13]。中值濾波可在有效去除噪點的同時保護圖像邊緣，數(shù)據(jù)重建更為真實[14]。因此，本文采用中值濾波對煤樣的原始CT 圖像進行預處理，處理結果如圖3 所示?？煽闯觯幚砗蟮腃T 圖像更清晰、噪點更少。

圖3 降噪前后CT 圖像對比Fig.3 CT image comparison before and after noise reduction

CT 圖像中，煤孔隙的識別依賴于閾值分割。常用的閾值分割法主要包括視覺判別、經(jīng)驗選取和模型反演[15]。本文采用文獻[16]提出的Bi-PTI 模型識別煤孔隙結構，Bi-PTI 模型對缺陷位置拓撲特征具有較好的識別能力，可真實還原孔裂隙拓撲體系的孔喉豐度和連通特性，對沙曲焦煤的適應性良好。Bi-PTI 模型數(shù)學公式為

式中：Φ（Gm）為煤樣孔隙率；Gm為孔隙的CT 灰度閾值，HU；ω為閾值權重；γ1，γ2分別為孔裂隙和基質(zhì)骨架的閾值均值，HU；ζ1，ζ2分別為孔裂隙和基質(zhì)骨架的灰度敏感系數(shù)。

利用圖形處理軟件提取CT 數(shù)據(jù)中不同孔隙率對應的灰度閾值，然后用Bi-PTI 模型進行擬合。將壓汞孔隙率代入模型并反演孔隙灰度閾值，擬合結果如圖4 所示，擬合參數(shù)見表1。

表1 Bi-PTI 模型擬合參數(shù)Table 1 Fitting parameters of Bi-PTI model

圖4 孔隙灰度閾值的Bi-PTI 擬合結果Fig.4 Bi-PTI fitting resules of pore gray threshold

結合圖4 和表1 可知，煤體CT 數(shù)據(jù)孔隙率與灰度閾值的Bi-PTI 擬合相關系數(shù)均超過0.99，擬合效果良好。5 號煤的灰度閾值為10 475 HU，6 號煤的灰度閾值略大于5 號煤，為12 923 HU。

2 數(shù)值模擬

2.1 數(shù)學模型

瓦斯在煤微觀孔隙中的連續(xù)滲流過程符合低雷諾數(shù)和低速流的條件[17]，因此滲流過程可采用Navier-Stokes（N-S）方程描述：

式中：ρ為瓦斯密度，kg/m3；v為流體速度，m/s；t為瓦斯擴散時間，s；?為拉普拉斯算子；P為瓦斯壓力，MPa；μ為動力黏度，Pa·s；F為體積力，N。

由質(zhì)量守恒定律可知，瓦斯在任意時刻流入和流出代表性體積（Representative Elementary Volume，REV）單元體截面的質(zhì)量差與單元體內(nèi)部的瓦斯質(zhì)量變化量相等[18]，故瓦斯在REV 單元體內(nèi)流動的連續(xù)性方程可描述為

依據(jù)前述假設條件，瓦斯密度ρ為常數(shù)，故可將式（3）簡化為

受煤巖體孔隙結構非均質(zhì)性影響，瓦斯在孔隙結構內(nèi)部流動過程中，任意一點的速度大小和方向均不相同，隨著空間位置和時間的變化而變化[19]，因此可將式（4）轉(zhuǎn)換為

聯(lián)立式（3）、式（5）、式（6）可得到N-S 方程的一般形式：

煤基質(zhì)內(nèi)瓦斯的解吸-擴散過程遵循含氣體源項的菲克第二定律：

式中：C為物質(zhì)的量濃度，與氣體所處的孔隙空間坐標位置（x,y,z）及擴散時間t有關，mol/L；D為氣體擴散系數(shù)，m2/s；S為氣體源項。

對于瓦斯氣體解吸-擴散耦合方程，式（7）中的氣體源項可用瓦斯?jié)舛入S時間的變化情況來描述[20]：

式中Cad為煤基質(zhì)內(nèi)瓦斯物質(zhì)的量濃度，mol/L。

式中：nad為煤巖基質(zhì)內(nèi)瓦斯氣體量，mol；Ve為數(shù)值模擬計算中網(wǎng)格單位的體積，m3；s為網(wǎng)格單位所包含孔隙表面積，nm2；M為單位孔隙表面積吸附的氣體量，mol；Vm為標準狀況下的氣體摩爾體積，22.4 L/mol。

單位孔隙表面積吸附的氣體量M為[12]

式中：Vad為煤基質(zhì)吸附的瓦斯體積，m3；Nsolid為固體體素總數(shù)量；Vvoxel為單位體素的體積，m3；ρtrue為煤巖密度，kg/m3；f為計算單元內(nèi)部孔隙總表面積，m2。

假設瓦斯只吸附在孔隙內(nèi)表面，當煤體內(nèi)部任意位置體素點的灰度g（x,y,z）滿足下列方程中的一個時，即可識別出孔隙空間內(nèi)表面[15]?？紫恫糠煮w素的灰度為1，其余部分的灰度為0。

煤基質(zhì)吸附的瓦斯體積Vad可用朗格繆爾擴展型方程表示：

式中：a為甲烷的朗格繆爾體積常數(shù)，m3/kg；b為甲烷的朗格繆爾壓力常數(shù)，MPa-1。

聯(lián)立式（7）—式（17）可得煤基質(zhì)內(nèi)部瓦斯的解吸-擴散耦合方程：

由于瓦斯具有可壓縮性，其狀態(tài)方程為

式中：n為氣體物質(zhì)的量，mol；Z為壓縮因子，與瓦斯溫度和壓力有關；R為氣體常數(shù)，J/（K·mol）；T為溫度，K。

在考慮瓦斯氣體壓縮效應的情況下，煤基質(zhì)內(nèi)瓦斯氣體解吸-擴散耦合方程為

2.2 模型及參數(shù)

基于重構的煤微觀孔隙結構，借助COMSOL 仿真軟件開展煤體瓦斯微觀解吸-擴散-滲流數(shù)值計算。模型選取尺寸為60×60×60 像素的立方體單元。煤孔隙幾何模型重建后需要優(yōu)化調(diào)試，主要包括消除重合線條、倒角、缺口等，經(jīng)優(yōu)化后的四面體網(wǎng)格可直接用于后續(xù)數(shù)值模擬計算中[13]。

數(shù)值計算所需的參數(shù)見表2，其中甲烷氣體相關參數(shù)均按照標準狀況下的氣體屬性設置，密度為0.714 kg/m3，動力黏度為1.80×10-5Pa·s，溫度為293 K。對于瓦斯解吸-擴散數(shù)值模擬，初始條件下，各煤體孔隙幾何模型中孔隙內(nèi)表面飽和瓦斯氣體的初始壓力為1×10-2Pa，且孔隙內(nèi)部壓力設置為0。瓦斯邊界條件如圖5 所示。對于瓦斯?jié)B流數(shù)值模擬，瓦斯?jié)B流過程為單向滲流。其中，在5 號煤中沿著x方向進行，其中孔隙外表面x+、x-分別設置為滲流進出口，對應瓦斯初始壓力分別為1.1，0.1 MPa；在6 號煤中沿著z方向進行，其中孔隙外表面z-、z+分別設置為滲流進出口，對應瓦斯初始壓力分別為1.1，0.1 MPa。將其余4 個孔隙外表面設置為自由滑移壁面，孔隙內(nèi)表面設置為無滑移壁面。

表2 數(shù)值模擬參數(shù)Table 2 Numerical simulation parameters

圖5 瓦斯邊界條件Fig.5 Gas boundary condition

3 結果及討論

3.1 煤微觀孔隙結構定量表征

借助Thermo Sientific Avizo 軟件對數(shù)字巖心進行重構，重構體尺寸均為500×500×1 000 像素。CT 數(shù)據(jù)體的REV 單元是探索孔裂隙精細結構特征的一種有效手段，通常用于地球物理領域的數(shù)字煤巖結構研究[20]。在重構的數(shù)字巖心結構中，選取感興趣的4 個目標點，選取的REV 單元尺寸均設置為200×200×200 像素，位置如圖6 所示。

基于各煤樣最佳孔隙灰度閾值，通過Avizo 構建REV 表征單元內(nèi)部的微觀孔隙空間結構，如圖7所示。分析可知，5 號煤含有大量孔隙，A2 內(nèi)部有多條連通的裂隙，A4 裂隙含量最少，均存在較多孔隙團。6 號煤含有大量孔隙，4 個REV 表征單元內(nèi)部的裂隙表現(xiàn)出不同形態(tài)，B1 內(nèi)部2 條呈V 字形的連通裂隙貫穿整個表征單元，B4 下方有多條呈片狀橫向平行分布的裂隙。

圖7 REV 表征單元孔隙空間結構Fig.7 Representative elementary volume(REV) characterisation of unit pore space structure

各煤樣REV 單元微觀孔隙等效直徑的分布規(guī)律如圖8 所示。分析可知，5 號煤4 個REV 單元的孔隙等效直徑分布范圍較為雜亂，4 μm 的孔隙均占比最高，其在A4 中的占比達26%。A2 的等效直徑分布范圍最窄，為1～8 μm，等效直徑為4 μm 的孔隙占比約為21%。與5 號煤相比，6 號煤微觀孔隙等效直徑分布范圍、優(yōu)勢孔隙分布均非常相似，表明5 號煤的非均質(zhì)性較強，其中等效直徑為3.5 μm 的孔隙占比約為23%，相對應的孔隙累計占比為57%。可見，5 號煤和6 號煤的孔隙等效直徑分布規(guī)律大致相同。

圖8 孔隙等效直徑分布規(guī)律Fig.8 Distribution law of pore equivalent diameter

煤微觀孔隙體積分布規(guī)律如圖9 所示。分析可知，5 號煤的孔隙體積主要分布在1～800 μm3，平均值為419 μm3，A1 中的優(yōu)勢孔最大體積為530 μm3，相對應孔隙體積在1～530 μm3的孔隙累計占比達63%。6 號煤的孔隙體積分布規(guī)律與5 號煤較為相似，不同之處在于6 號煤中孔隙體積為260 μm3的孔隙占比較高，表明6 號煤體積的小孔隙占主導地位?？梢姡? 種煤樣的孔隙等效體積分布范圍不存在明顯差異。

圖9 孔隙體積分布規(guī)律Fig.9 Distribution law of pore volume

3.2 瓦斯解吸-擴散過程

瓦斯在煤微觀孔隙REV 單元中的壓力分布情況如圖10 所示?？煽闯觯S著瓦斯解吸-擴散時間的增加，瓦斯壓力自REV 單元中心至邊緣逐漸降低，同時在不同時間及位置上瓦斯壓力分布規(guī)律差異明顯。造成瓦斯壓力分布差異性的原因在于各REV單元中孔隙與喉道的半徑、長度、形狀、連通性能不同。

圖10 瓦斯壓力空間分布Fig.10 Space distribution of gas pressure

為進一步分析瓦斯解吸-擴散過程中瓦斯壓力在不同微觀孔隙系統(tǒng)中的分布情況，對切片中瓦斯壓力場進行定量分析。分別截取各REV 單元體在x軸5，10，15，20，25，30 μm 位置的y-z平面，觀察不同位置切片壓力場隨時間變化情況，如圖11 所示。

圖11 REV 單元y-z 截面孔隙壓力分布Fig.11 Pore pressure distribution in y-z section of REV unit

分析可知，在同一REV 切片中，瓦斯在解吸-擴散時由壓力大的中心位置向壓力較小的邊緣位置擴散。同一時間，在不同切片中，中心位置（10，15，20，25 μm）切片瓦斯壓力相對較高，而邊緣位置（5，30 μm）切片瓦斯壓力相對較低?？梢?，瓦斯氣體在孔隙系統(tǒng)中由內(nèi)而外地逐步完成解吸-擴散。此外，尺寸較大的連通孔隙能夠為氣體分子提供不同方向上的多元化運動空間，進而提高瓦斯解吸-擴散速率。但連通孔隙的小尺寸結構限制了氣體分子空間擴散的廣度，極大阻礙瓦斯的解吸-擴散進程。

3.3 瓦斯?jié)B流過程

在煤三維孔隙空間結構中，由于孔隙壓力分布、孔隙通道連通性的差異性，造成瓦斯氣體在煤巖體孔隙空間中的滲流速度大小和方向發(fā)生明顯變化。為研究瓦斯?jié)B流速度變化規(guī)律，保持瓦斯初始壓力和擴散系數(shù)不變，分別截取各REV 單元不同時間節(jié)點（1，10，20，40 s）的速度場流線分布情況，如圖12所示，圖中不同顏色的箭頭指示線條表示不同時間節(jié)點瓦斯在煤巖體中的滲流過程，顏色由深藍色逐漸過渡到深紅色表示瓦斯流動速度不斷增大。分析發(fā)現(xiàn)，5 號煤中瓦斯?jié)B流速度流線的分布范圍、密度較6 號煤更大，表明5 號煤的連通孔隙數(shù)量更多，瓦斯?jié)B流通道更加豐富。5 號煤和6 號煤滲流過程中瓦斯在各煤體三維孔隙空間內(nèi)的滲流速度大小和方向表現(xiàn)出顯著的差異性。在諸如5 號煤的強非均質(zhì)連通孔隙結構中，瓦斯?jié)B流分散而高效，能通過廣泛溝通煤基質(zhì)完成氣體由擴散到滲流的轉(zhuǎn)變，提升瓦斯傳質(zhì)效率。然而，在諸如6 號煤的弱非均質(zhì)連通孔隙結構中，氣體滲流路徑單一、流線集中，滲流傳質(zhì)阻力較大，氣體分子由擴散到滲流的轉(zhuǎn)變效率低，不利于瓦斯高效運移。

圖12 瓦斯?jié)B流速度流線分布Fig.12 Velocity streamline distribution of gas seepage

4 結論

1）從微觀角度看，瓦斯在解吸-擴散時由孔隙中心位置向邊緣位置進行。中心位置瓦斯壓力相對較高，而邊緣位置瓦斯壓力相對較低。大尺寸連通孔隙結構能夠為氣體分子提供不同方向上的多元化運動空間，削弱尺寸效應對氣體分子擴散廣度的影響，進而促進瓦斯解吸-擴散速率。

2）瓦斯在微觀孔隙空間結構內(nèi)的滲流速度大小和方向表現(xiàn)出明顯的差異性，這主要是由于不同煤體微觀孔隙系統(tǒng)存在高度非均質(zhì)性。強非均質(zhì)連通孔隙結構中的滲流傳質(zhì)阻力較小，能通過廣泛溝通煤基質(zhì)完成氣體由擴散到滲流的高效轉(zhuǎn)變，提升瓦斯傳質(zhì)效率。研究結果能從微觀角度豐富煤體瓦斯運移理論，為瓦斯抽采工程實踐提供一定理論基礎。