王登科，李文睿，魏建平,4，于 充，張凱旋

(1.河南理工大學(xué) 河南省瓦斯地質(zhì)與瓦斯治理重點實驗室-省部共建國家重點實驗室培育基地，河南 焦作 454000; 2.中國礦業(yè)大學(xué) 深部巖土力學(xué)與地下工程國家重點實驗室，江蘇 徐州 221116; 3.河南理工大學(xué) 安全科學(xué)與工程學(xué)院，河南 焦作 454000; 4.煤炭安全生產(chǎn)河南省協(xié)同創(chuàng)新中心，河南 焦作 454000)

多孔介質(zhì)中流體的質(zhì)量傳遞、滲透率、流動阻力和擴散性能均與介質(zhì)內(nèi)部的微觀結(jié)構(gòu)密切相關(guān)。煤作為一種具有從微孔到大孔范圍的多孔介質(zhì)材料，其內(nèi)部流體的傳熱傳質(zhì)一直以來都是研究的熱點和難點。由于煤體內(nèi)部存在有機質(zhì)和無機質(zhì)分布的不均勻性，流動通道呈現(xiàn)出復(fù)雜多變的結(jié)構(gòu)特征，多采用理想水動力管束模型和宏觀流動方程研究氣體傳輸機理，難以對真實孔隙流動通道內(nèi)的傳輸機理進行準(zhǔn)確定義，隨著科研技術(shù)的不斷進步，這一限制條件正在被打破，已有專家學(xué)者采用先進的SEM[1]、顯微CT[2]、3D打印[3-4]和原子力顯微鏡[5]等技術(shù)手段對煤體、頁巖和砂巖內(nèi)部孔裂隙結(jié)構(gòu)進行了定性和定量描述。此外，已有學(xué)者證明多孔介質(zhì)煤巖材料內(nèi)部的孔隙結(jié)構(gòu)可以采用分形幾何理論進行表征，輔以計盒法可求得其分形維數(shù)[6-7]。XU和YU[8]的研究結(jié)果表明孔隙結(jié)構(gòu)不僅是分形的，而且分形維數(shù)還可表示為結(jié)構(gòu)參數(shù)的函數(shù);并在此基礎(chǔ)上建立了多孔介質(zhì)孔隙毛細管束分形模型，用于多孔介質(zhì)的分形維數(shù)與滲透率、擴散系數(shù)之間關(guān)系的研究。此外，基于隨機分形Menger海綿體模型[9]，王恩元和何學(xué)秋還建立了煤巖孔隙結(jié)構(gòu)的隨機分形模型[10]。李子文等[11-12]研究發(fā)現(xiàn)，只有當(dāng)孔隙壓力低于10 MPa時以壓汞法測量出的孔隙才具有隨機分形模型的特征，此時對應(yīng)的孔隙半徑約為60 nm。鄭仟[13]、員美娟[14]、李玉丹[15]、王洪磊[16]等通過理論和實驗研究表明，毛細管束分形模型對孔隙結(jié)構(gòu)的表征更加準(zhǔn)確且適用的孔隙尺寸范圍更廣，由此建立的多孔介質(zhì)內(nèi)部流體滲透率模型也更加準(zhǔn)確[17]。國外學(xué)者通過類比孔隙結(jié)構(gòu)分形，發(fā)現(xiàn)孔隙內(nèi)部表面粗糙元或者接觸點也滿足一定的分形標(biāo)度規(guī)律[18-19]。楊姍姍[17,20]建立了錐體的粗糙元分形模型，得到了摩擦因子和泊肅葉數(shù)的分析解模型，并討論了孔隙結(jié)構(gòu)參數(shù)對泊肅葉數(shù)的影響，這使得微觀孔隙結(jié)構(gòu)的分形研究得到了進一步的發(fā)展和延伸。

筆者以粗糙微通道中孔隙形態(tài)和內(nèi)壁粗糙元分形分析為基礎(chǔ)，推導(dǎo)了粗糙微納米孔隙中的瓦斯氣體傳輸方程，將瓦斯氣體等效視滲透率模型與實驗結(jié)果進行對比，從孔隙結(jié)構(gòu)參數(shù)的精確求解和分子流動機理兩方面印證了筆者所建瓦斯氣體傳輸方程的合理性和準(zhǔn)確性。

1 煤樣孔隙結(jié)構(gòu)實驗結(jié)果

利用河南理工大學(xué)的JSM-6390LV掃描電鏡(圖1)對煤樣的孔隙進行了掃描，結(jié)果如圖2所示。

圖1 JSM-6390LV掃描電鏡

根據(jù)電鏡掃描得到的煤體孔隙內(nèi)部具體形貌細節(jié)(圖2)，可以發(fā)現(xiàn)煤的孔隙及其內(nèi)部粗糙結(jié)構(gòu)與納米管束分形模型的基本假設(shè)相吻合[21]。因此，借鑒前人的研究成果[22]，可將孔隙通道簡化為圓管，內(nèi)表面粗糙元以簡化后的半球體代替，如圖3所示，并可按照分形標(biāo)度定律通過結(jié)構(gòu)參數(shù)分形對孔隙和粗糙元進行分形表征。

圖2 煤中粗糙孔隙微通道形貌結(jié)構(gòu)

圖3 煤中孔隙微通道結(jié)構(gòu)簡化模型

2 粗糙微通道的分形表征

2.1 孔隙分布的分形表征

YU[23]的研究表明，多孔介質(zhì)中孔隙孔徑分布與孔隙結(jié)構(gòu)滿足分形特征，孔隙數(shù)量分布函數(shù)表達式為

N(≥dp)=(dm/dp)Df

(1)

式中，dp為多孔介質(zhì)納米孔的直徑，m;N為直徑≥dp的孔隙數(shù)量;dm為多孔介質(zhì)孔隙直徑最大值，m;Df為孔隙直徑的分形維數(shù)，在二維和三維空間內(nèi)的取值分別為0

在式(1)中可把孔隙數(shù)量表達式看作一個連續(xù)的可微分函數(shù)，則直徑為dp～(dp+ddp)的孔隙數(shù)量為

(2)

式中，-dN>0 表明孔隙數(shù)量隨著孔隙半徑減小而增加的關(guān)系。

對于具有自相似特征的孔隙結(jié)構(gòu)，其直徑大小[23]可用式(3)計算孔隙分形維數(shù)

(3)

式中，χ為歐幾里德維數(shù)，在二維平面中χ=2，在三維平面中χ=3;φ為多孔介質(zhì)的孔隙度，%;σ為孔隙最小直徑與最大直徑之比。

2.2 孔隙內(nèi)部粗糙元分形分析

對于粗糙元底面直徑大小的分形維數(shù)表達式可類比于式(3)，但粗糙元的分形維數(shù)D的取值為0

(4)

式中，ε為孔隙粗糙體微元底面直徑的最小值與最大值之比，ε=dmin/dmax;φ′為在一個分形集合中，所有半球體的底面總面積占整個分形集合單元總面積的百分數(shù)。

同理，直徑為dr～d+ddr的粗糙元數(shù)量為

(5)

其中，dr為粗糙元底部直徑，m;dmax為粗糙元底部直徑最大值，m;N′為直徑≥dr的粗糙元數(shù)量;-dN’>0 表明粗糙元數(shù)量隨著孔隙半徑減小而增加的關(guān)系;D為粗糙元的分形維數(shù)。

一個分形集合的單元中所有半球體的體積可以通過從最小粗糙元到最大粗糙元的積分求和得到:

(6)

一個分形集合單元中所有半球體的底面積總和為

(7)

這里，考慮到孔隙直徑屬于二維平面，所以χ=2。將式(5)代入式(7)中，可得到

(8)

因此，由式(8)可得出一個孔隙結(jié)構(gòu)單元中的總面積為

(9)

所以，在孔隙內(nèi)部半球體粗糙元的有效平均高度為

(10)

3 粗糙孔隙內(nèi)氣體傳輸方程

3.1 滑移流動方程

微孔隙通道中的瓦斯氣體傳輸不同于傳統(tǒng)的水動力學(xué)流體傳輸規(guī)律。尤其當(dāng)微通道孔隙的孔徑與瓦斯氣體分子自由程處于同一數(shù)量級時，瓦斯氣體分子與孔隙壁面之間的碰撞愈發(fā)劇烈，氣體稀薄效應(yīng)也更加顯著，其流動速度明顯受到孔隙結(jié)構(gòu)影響[25-26]。因此，學(xué)者們提出了以Maxwell分子碰撞散射模型去揭示稀薄氣體分子在微通道孔隙中的真實傳輸情況。其基本思想為:分子的碰撞運動由兩種反射組合而成，假設(shè)來流分子中發(fā)生漫反射的部分為t，則剩余部分發(fā)生鏡面反射。

假設(shè)w0是氣體相對于壁面的速度，則在流動方向存在法向速度梯度，dw/dz，法向方向速度為

(11)

另外，根據(jù)Maxwell散射模型可得孔壁表面分子速度w0為

(12)

聯(lián)立式(11)和(12)可得

(13)

需要說明的是，在微孔隙中，除了分子與壁面的碰撞以外，分子與分子之間還發(fā)生著不同程度的碰撞。熱力學(xué)理論認為，單位時間內(nèi)分子碰撞的次數(shù)等于分子與分子之間碰撞次數(shù)與分子與壁面碰撞次數(shù)之和[27]，因此可得出

(14)

式中，T1為單位時間內(nèi)分子與分子碰撞消耗的時間;T2為單位時間內(nèi)分子與壁面碰撞消耗的時間;Ttotal為分子碰撞的某一單位時間。

假設(shè)分子發(fā)生碰撞的初始熱運動速度均為v，則各形式碰撞次數(shù)可表示為

(15)

式中,λ為氣體分子平均自由程;r為圓管半徑。

將式(15)代入(14)中可得

(16)

由此，可以得到分子與壁面的碰撞次數(shù)占總碰撞次數(shù)的比例系數(shù)c為

(17)

其中，Kn為克努森數(shù)，等于分子平均自由程與孔隙半徑之比，無量綱，Kn越大表示氣體傳輸過程中分子碰撞越劇烈。

在式(11)中分子碰撞后在法相方向上的位移可表示為

(18)

進一步，將式(18)代入式(13)可得粗糙微通道壁面分子的速度為

(19)

如圖4所示，在圓柱徑向x～x+dx，軸向r～r+dr方向流體受到壁面的切應(yīng)力為

(20)

式中，μ為氣體動力黏度，Pa·s。

圖4 圓管內(nèi)微元體傳輸結(jié)構(gòu)

(21)

根據(jù)式(21)進行一次積分后可得到

(22)

通過對式(22)進行再次積分進一步得到孔隙流體速度

(23)

當(dāng)r=0時，壁面分子流速為定值，則a=0;當(dāng)滿足真實孔隙半徑大小邊界條件，即式(19)中r=R時，通過聯(lián)立式(23)，管壁面的速度滿足:

(24)

由此可得常數(shù)b為

(25)

再將式(25)代入式(23)進而可得到微通道內(nèi)某一質(zhì)點的流體速度表達式

(26)

當(dāng)某一孔隙半徑為定值r0時，通過單個微通道的體積流量可以表示為

(27)

根據(jù)孔隙通道的真實情況，在此對氣體傳輸通道的有效半徑進行細化可得

(28)

式中，re為粗糙孔隙微通道內(nèi)的有效半徑，m。

因此，微通道內(nèi)具有滑移流效應(yīng)的平均流速可表示為

(29)

相應(yīng)地，考慮到真實氣體效應(yīng)影響，孔隙微通道內(nèi)的質(zhì)量流量可以表示為

(30)

式中，J為粗糙孔隙微通道內(nèi)的氣體質(zhì)量通量，kg/(m2·s);Q為氣體體積流量，m3/s;ρ為游離瓦斯氣體真實密度，kg/m3。

同時考慮到多孔介質(zhì)的孔隙率影響，式(30)氣體質(zhì)量通量方程可轉(zhuǎn)化為等量視滲透率為

(31)

式中,Ks為連續(xù)性滑流氣體傳輸?shù)攘恳暆B透率，m2;φ為多孔介質(zhì)孔隙率，%。

3.2 非連續(xù)流氣體擴散傳輸方程

受煤體內(nèi)孔徑連通性的影響，孔隙壓力處于一個動態(tài)松弛變化的過程，導(dǎo)致孔隙內(nèi)的氣體擴散傳輸呈現(xiàn)出由瞬態(tài)激發(fā)到緩慢衰減的一個漫長過程[28]。若不考慮煤的理化性質(zhì)對氣體吸附性能差異性的影響[29]，單純研究孔隙內(nèi)可流動氣體的傳輸特性就變得直觀具體，僅就某一孔隙作為氣體流動通道載體來說，穩(wěn)定壓力梯度驅(qū)動下的氣體擴散傳輸可通過理論公式進行定量研究。

對于非連續(xù)流的分子擴散傳輸，主要以分子間碰撞和分子與壁面碰撞兩種形式存在。在兩種分子擴散同時存在情況下，以濃度梯度為驅(qū)動力的擴散傳輸通量方程[25]可表示為

(32)

同理，氣體擴散質(zhì)量通量傳輸方程式(31)轉(zhuǎn)化為等量視滲透率為

(33)

式中，Kd為分子擴散傳輸通量的等量視滲透率，m2。

3.3 微納米孔隙氣體視滲透率模型

將滑流方程和擴散方程的等效視滲透率相加即可得到納米孔氣體傳輸統(tǒng)一視滲透率數(shù)學(xué)模型:

(34)

式中，Kapp為考慮粗糙孔隙分形特征的瓦斯氣體傳輸統(tǒng)一視滲透率，m2。

對于未經(jīng)修正的多孔介質(zhì)固有滲透率，其表達式為

(35)

(36)

對于A的取值，AGRAWAL和PRABHU[30]系統(tǒng)性的研究了參數(shù)A中t的取值，建議取值為0.9較為合適，即A=11/9;此后，DARABI[31]等經(jīng)過數(shù)值模擬、計算分析與實驗對比驗證了取值的合理性，因此在下文模型計算過程中則直接取該值代入計算。

圖5 粗糙形貌的3D掃描成像儀

3.4 微納米孔隙相對粗糙度的測試與求解

利用河南理工大學(xué)的Optical Tensiometer設(shè)備(圖5)對取自煤體孔隙內(nèi)壁的煤片進行微米級3D形貌結(jié)構(gòu)掃描分析，利用自帶的OneAttension軟件進行詳盡的定量化表征，二維成圖效果如圖6所示，三維成圖效果如圖7所示。

圖6 煤體孔隙壁面粗糙形貌結(jié)構(gòu)

圖7 孔隙壁面粗糙體三維分布特征

圖6中左圖為粗糙表面二維結(jié)構(gòu)圖，右圖為對應(yīng)的圖片二值化處理結(jié)果。在二維平面中，圖中白色斑點部分為光滑無粗糙區(qū)域，暗部為粗糙結(jié)構(gòu)分布區(qū)域，通過MATLAB編程技術(shù)求得粗糙元投影底面總面積占整個平面比例φ′為77.08%。由于在二維孔隙直徑內(nèi)部粗糙元分形維數(shù)D=2表示粗糙元鋪滿整個平面，所以此處根據(jù)77.08%比例關(guān)系求得對應(yīng)D=1.541 6。

圖7中以粗糙體峰值0高度作為孔隙壁面分割基準(zhǔn)面，以上部分為上凸粗糙體，以下為下凸粗糙體。由于第2節(jié)中將粗糙體假想為半球微元體，所以粗糙體微元底面直徑最小值與最大值之比ε就轉(zhuǎn)化為粗糙元高度之比。在實際巖石孔隙中，存在孔隙內(nèi)壁粗糙元直徑大小與孔隙半徑相當(dāng)，此時粗糙元最大直徑dmax即為孔隙半徑。因此，在φ′=77.08%，D=1.541 6，dmax=r的賦值條件下，可求出相對粗糙度值，求解公式為

(37)

為便于計算這里將粗糙元最低輪廓(圖7中淺色部分所構(gòu)成的宏觀輪廓)偏離基準(zhǔn)面的算術(shù)平均值作為粗糙體理想化最小高度值(該值由設(shè)備自帶的OneAttension軟件處理得到)，相對粗糙度計算結(jié)果見表1。

表1 壁面粗糙結(jié)構(gòu)的參數(shù)信息

Table 1 Parameters about the structure roughness of inner surface in pore

表面輪廓偏離分割基準(zhǔn)面的算術(shù)平均差/μm最大峰值/μm粗糙元高度之比ε相對粗糙度γ4.17446.2630.0090.352

本文所提出的相對粗糙度計算公式(式(37))，相較于DARABI等[31]提出的孔隙粗糙度經(jīng)驗性表達式(分子直徑與孔隙局部直徑比值的指數(shù)式修正系數(shù)，見表2中Darabi模型)，以真實孔隙簡化后的微結(jié)構(gòu)分形幾何解析表達式，更能直觀反映出煤體中微納米孔隙的結(jié)構(gòu)參數(shù)，能更好地反映粗糙元底面直徑、排布方式、內(nèi)壁面積覆蓋度等因素對氣體傳輸通道有效半徑的影響。

3.5 視滲透率模型驗證性實驗

對于孔隙瓦斯氣體傳輸方程，采用等效的滲透率方程驗證是最為直接有效的方法[32]。因此，在未經(jīng)修正的多孔介質(zhì)固有滲透率模型基礎(chǔ)之上，學(xué)者們又提出了多種較為精確的視滲透率模型，具體見表2。

本文孔隙滲透率驗證實驗步驟為:首先以納米薄膜材料為孔隙載體，通過SEM掃描技術(shù)對納米孔隙薄膜孔徑分布進行表征(圖8);然后利用MATLAB對圖片進行二值化處理求得孔隙率(參數(shù)見表3);并利用河南理工大學(xué)PMI微滲透測試設(shè)備(圖9)進行滲透率實驗。最后，將本文視滲透率模型和幾種典型視滲透率模型的計算結(jié)果與實驗結(jié)果進行對比，對比結(jié)果如圖10所示。

表2 幾種典型孔隙流體的視滲透率模型

Table 2 Several typical apparent permeability models of gas flow in pore

模型簡稱視滲透率模型數(shù)學(xué)表達式闡述Darabi模型[31]Kapp=é?êê1+4AKn+64φKn3τπδr()Df-2ù?úúKDacy,δr()Df-2為考慮粗糙度影響的修正系數(shù),Df=2.5,孔隙迂曲度τ=2.2,甲烷分子直徑δ=0.38nmCivan模型[32]Kapp=[1+β(Kn)Kn]1+4Kn1+Kn()KDacy,β(Kn)=6415π1+0.199Kn0.365()Beskok模型[33]Kapp=[1+α(Kn)Kn]1+4Kn1+Kn()KDacy,α(Kn)=12815π2arctan(4Kn0.4)

圖8 納米孔隙薄膜SEM圖和對應(yīng)的二值化圖

表3 納米薄膜物理特性參數(shù)

圖9 納米薄膜微滲透性測試設(shè)備

圖10 理論模型與實驗數(shù)據(jù)比較結(jié)果

不論是微米級的孔隙通道還是納米級孔隙通道，其內(nèi)部形貌結(jié)構(gòu)在一定程度上都較為類似。PFUND等[34]通過微米級管道的流體分析，得出相對粗糙度隨著孔隙尺度的減小而增大，孔隙壁面粗糙將對流體傳輸造成顯著影響。圖10中的模型對比結(jié)果有效驗證了這一觀點，相較于其他理論模型不考慮粗糙度影響或者采用粗糙度經(jīng)驗性公式的視滲透率模型，本文模型更加精確可靠。同時從圖10中也可以看出，本文所建立的視滲透率傳輸模型與實驗結(jié)果吻合效果良好，充分證明了用分形方法表征孔隙微觀結(jié)構(gòu)進而研究粗糙微納米氣體傳輸?shù)暮侠硇院蜏?zhǔn)確性。

4 結(jié) 論

(1)以分形表征得到的微納米孔隙的粗糙元有效平均高度，精細地反映了微納米孔隙表面的粗糙結(jié)構(gòu)機制，并且每個結(jié)構(gòu)參數(shù)的物理意義明確，能夠揭示影響有效孔隙半徑的微觀結(jié)構(gòu)機制，這對于以分形拓撲理論構(gòu)建自仿射粗糙微通道進行煤層瓦斯氣體流動模擬研究具有重要意義。

(2)相較于直接在氣體傳輸方程中采用粗糙度經(jīng)驗公式修正微納米孔隙瓦斯氣體傳輸通量的方法，利用粗糙形貌3D掃描成像設(shè)備，對孔隙粗糙表面微結(jié)構(gòu)進行定量表征，所求得的相對粗糙度值更加科學(xué)可靠。

(3)所建立的基于分形表征的粗糙孔隙瓦斯氣體傳輸方程綜合考慮了孔隙粗糙度、氣體速度滑移和分子擴散對瓦斯氣體流動的重要影響，計算結(jié)果更為準(zhǔn)確合理，有效揭示了瓦斯氣體微觀傳輸機制。