煤孔隙結(jié)構(gòu)對(duì)瓦斯解吸-擴(kuò)散-滲流過(guò)程的影響

賈男

（1.中煤科工集團(tuán)沈陽(yáng)研究院有限公司，遼寧 撫順 113122；2.煤礦安全技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 撫順 113122）

0 引言

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

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

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

1 樣品及實(shí)驗(yàn)

1.1 煤樣處理

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

圖1 實(shí)驗(yàn)樣品Fig.1 Experimental samples

1.2 實(shí)驗(yàn)方法

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

圖2 實(shí)驗(yàn)儀器Fig.2 Experimental instruments

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

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

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

CT 圖像中，煤孔隙的識(shí)別依賴于閾值分割。常用的閾值分割法主要包括視覺判別、經(jīng)驗(yàn)選取和模型反演[15]。本文采用文獻(xiàn)[16]提出的Bi-PTI 模型識(shí)別煤孔隙結(jié)構(gòu)，Bi-PTI 模型對(duì)缺陷位置拓?fù)涮卣骶哂休^好的識(shí)別能力，可真實(shí)還原孔裂隙拓?fù)潴w系的孔喉豐度和連通特性，對(duì)沙曲焦煤的適應(yīng)性良好。Bi-PTI 模型數(shù)學(xué)公式為

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

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

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

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

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

2 數(shù)值模擬

2.1 數(shù)學(xué)模型

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

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

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

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

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

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

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

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

對(duì)于瓦斯氣體解吸-擴(kuò)散耦合方程，式（7）中的氣體源項(xiàng)可用瓦斯?jié)舛入S時(shí)間的變化情況來(lái)描述[20]：

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

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

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

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

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

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

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

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

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

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

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

2.2 模型及參數(shù)

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

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

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

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

3 結(jié)果及討論

3.1 煤微觀孔隙結(jié)構(gòu)定量表征

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

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

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

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

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

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

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

3.2 瓦斯解吸-擴(kuò)散過(guò)程

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

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

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

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

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

3.3 瓦斯?jié)B流過(guò)程

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

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

4 結(jié)論

1）從微觀角度看，瓦斯在解吸-擴(kuò)散時(shí)由孔隙中心位置向邊緣位置進(jìn)行。中心位置瓦斯壓力相對(duì)較高，而邊緣位置瓦斯壓力相對(duì)較低。大尺寸連通孔隙結(jié)構(gòu)能夠?yàn)闅怏w分子提供不同方向上的多元化運(yùn)動(dòng)空間，削弱尺寸效應(yīng)對(duì)氣體分子擴(kuò)散廣度的影響，進(jìn)而促進(jìn)瓦斯解吸-擴(kuò)散速率。

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