王昌進(jìn)，張 賽，徐靜磊

（昆明理工大學(xué)機(jī)電工程學(xué)院，昆明 650500）

0 引言

多孔介質(zhì)中的氣體擴(kuò)散問題是能源開采、燃料電池、化學(xué)工程等領(lǐng)域研究的重點(diǎn)，但是氣體在多孔介質(zhì)中的擴(kuò)散包含了很多復(fù)雜的傳質(zhì)機(jī)理，如何建立多孔介質(zhì)中氣體有效擴(kuò)散系數(shù)模型成為研究熱點(diǎn)［1］。

傳統(tǒng)的物質(zhì)擴(kuò)散定律一般由菲克定律來描述，該定律描述了物質(zhì)擴(kuò)散通量與濃度梯度之間的定量關(guān)系。Cao 等［2］通過建立多孔介質(zhì)三維模型，證明在穩(wěn)態(tài)情況下，菲克定律能夠在多孔介質(zhì)氣體擴(kuò)散中適用，然而該模型使用的是100×100×100 的立方網(wǎng)絡(luò)模型，對于多孔介質(zhì)的不規(guī)則性和復(fù)雜性方面的描述不夠準(zhǔn)確。氣體在多孔介質(zhì)中復(fù)雜的輸運(yùn)機(jī)理與多孔介質(zhì)孔隙大小和結(jié)構(gòu)差異有很大關(guān)系，根據(jù)多孔介質(zhì)孔隙大小和結(jié)構(gòu)的不同，多孔介質(zhì)中的氣體流動(dòng)與擴(kuò)散常常劃分為4 個(gè)流動(dòng)區(qū)域：連續(xù)介質(zhì)區(qū)、滑移區(qū)、過渡區(qū)和自由分子區(qū)［3］。Beskok等［4］通過大量實(shí)驗(yàn)給出了一維圓管中氣體體積流量的表達(dá)式，該表達(dá)式引入克努森數(shù)kn和滲透率修正因子f（kn），可以適用于所有流動(dòng)區(qū)域。符東宇等［5］運(yùn)用格子Boltzmann 方法研究了頁巖氣藏的滲流規(guī)律，發(fā)現(xiàn)克努森數(shù)kn是影響頁巖氣滲流的主要參數(shù)。Roy 等［6］使用Ar，N2，O2和Al2O3作為實(shí)驗(yàn)材料進(jìn)行了擴(kuò)散實(shí)驗(yàn)，引入了克努森擴(kuò)散系數(shù)，建立了微納米孔道中的氣體擴(kuò)散方程，但是該方程只能適用于自由分子區(qū)。張烈輝等［7］考慮了多重運(yùn)移機(jī)制，建立了頁巖氣藏的綜合滲流模型。李亞雄等［8］通過研究頁巖氣的滑脫和擴(kuò)散機(jī)理，提出了“壁聯(lián)擴(kuò)散”的新概念來表征克努森擴(kuò)散和表面擴(kuò)散的總效應(yīng)，但該方法的數(shù)學(xué)模型與論證過程還有待完善。在一些科學(xué)研究和實(shí)際情況中，例如川南地區(qū)頁巖孔隙結(jié)構(gòu)和川西坳陷致密砂巖孔隙結(jié)構(gòu)，它們的孔道直徑很小，達(dá)到了微納米級別，此時(shí)，孔隙結(jié)構(gòu)顯示出分形的特性。陳居凱等［9］，朱漢卿等［10］和鄧浩陽等［11］均使用分形理論研究了頁巖與致密砂巖等致密多孔介質(zhì)的孔隙結(jié)構(gòu)。目前，也有很多文獻(xiàn)將分形理論運(yùn)用到氣體滲透和擴(kuò)散中［12-14］。Zheng等［15］基于分形理論和毛細(xì)管束模型，得到了氣體的有效擴(kuò)散系數(shù)模型。該模型能與已有實(shí)驗(yàn)數(shù)據(jù)相吻合，并且能適用于不同的多孔介質(zhì)。Shi 等［16］在考慮了2 種運(yùn)輸機(jī)制的情況下，運(yùn)用分形理論建立了H2在燃料電池氣體擴(kuò)散層中的有效擴(kuò)散系數(shù)模型，并且分析了分形理論相關(guān)參數(shù)（孔隙度、分形維數(shù)等）的影響。Zhang［17］考慮了3 種多孔介質(zhì)中的氣體流動(dòng)機(jī)理，得到了滲透率和氣體有效擴(kuò)散分形模型，但忽略了氣體的滑脫效應(yīng)。Woignier 等［18］研究了分形多孔介質(zhì)中氣體在滑移區(qū)的擴(kuò)散機(jī)理并建立了滑移模型。張賽等［19］和牟新竹等［20］在研究多尺度分形多孔介質(zhì)氣體擴(kuò)散時(shí)加入了孔隙連通性參數(shù)，得到了新的數(shù)學(xué)模型。以上研究大多只考慮了1 種或2 種氣體擴(kuò)散機(jī)理，不能完整地描述氣體在多孔介質(zhì)中的擴(kuò)散問題。

基于Beskok 等［4］給出的表達(dá)式、迂曲毛細(xì)管束模型、分形理論和菲克定律，建立氣體在多孔介質(zhì)中擴(kuò)散的分形模型；該模型包含面積分形維數(shù)、迂曲分形維數(shù)、孔隙度、最大最小孔隙直徑、氣體密度、黏度、摩爾質(zhì)量等參數(shù)，然后分析這些參數(shù)與氣體有效擴(kuò)散系數(shù)之間的函數(shù)關(guān)系與變化趨勢，以期能夠正確地描述氣體在多孔介質(zhì)中的擴(kuò)散規(guī)律，并應(yīng)用于油氣田的開發(fā)。

1 多孔介質(zhì)的分形特征

1.1 孔隙數(shù)目與橫截面積描述

在具有分形特性的多孔介質(zhì)中，孔隙直徑大于等于λ的孔隙數(shù)量的總和與孔隙直徑大小滿足下面的分形標(biāo)度關(guān)系［21］

式中：N為孔隙長度大于孔隙直徑的多孔介質(zhì)的孔隙總數(shù)目；L為毛細(xì)管長度，m；λmax與λ分別為最大孔隙直徑和孔隙直徑，m；Df為無因式多孔介質(zhì)的面積分形維數(shù)。

在真實(shí)多孔介質(zhì)中，孔隙數(shù)目十分巨大，根據(jù)統(tǒng)計(jì)理論和方法，可以將式（1）看成是連續(xù)和可微的函數(shù)，對式（1）進(jìn)行微分可以得到

式（2）得到了在(λ,λ+dλ) 的孔隙數(shù)目，表明了孔隙數(shù)目隨著孔隙直徑的增大而減少，-dN＞0。由于截面上的孔隙可以看作是直徑不同的圓，因此通過式（2），用微元法可得到分形集中截面上總的孔隙面積AP為

式中：AP為截面上的總的孔隙面積，m2；λmin為最小孔隙直徑，m。

則總的截面積At可以計(jì)算得到

式中：At為總的截面積，m2；ε為面孔隙度。

一般情況下，體孔隙度和面孔隙度不相等，兩者之間相差一個(gè)迂曲度因子

式中：Φ為體孔隙度；為迂曲度因子［1］。在計(jì)算多孔介質(zhì)相關(guān)參數(shù)時(shí)通常使用的是體孔隙度，通過式（5）可以實(shí)現(xiàn)面孔隙度與體孔隙度的轉(zhuǎn)化。

上述式（1）—（5）描述了多孔介質(zhì)內(nèi)部結(jié)構(gòu)和孔隙分布情況，是描述多孔介質(zhì)內(nèi)部輸運(yùn)物理現(xiàn)象的重要公式。

1.2 孔隙迂曲特征描述

假設(shè)多孔介質(zhì)是由一束橫截面積不同的迂曲毛細(xì)管道組成，在該模型中，單根迂曲毛細(xì)管道的長度L(λ) 與管道直徑λ的關(guān)系滿足分形冪規(guī)律，表示為

式中：L0為毛細(xì)管流動(dòng)方向的特征長度，m。根據(jù)分形理論，L0可以通過L0=計(jì)算得到。須要注意的是L0不是多孔介質(zhì)的實(shí)際邊長［21］。Dt是迂曲分形維數(shù)，Dt越高表示毛細(xì)管道越彎曲。在二維平面內(nèi)，當(dāng)Dt=1 時(shí)，表示毛細(xì)管道是直管道，當(dāng)Dt=2 時(shí)，表示毛細(xì)管道非常彎曲，已經(jīng)填充了整個(gè)二維平面。根據(jù)式（6），迂曲度可以定義為

2 多孔介質(zhì)中氣體擴(kuò)散的分形模型

氣體在多孔介質(zhì)中的擴(kuò)散問題影響因素眾多，其中孔隙結(jié)構(gòu)為主要的影響因素。當(dāng)孔隙直徑不同時(shí)，氣體擴(kuò)散的規(guī)律也不同。Beskok 等［4］通過實(shí)驗(yàn)給出了描述氣體在一維圓管中氣體體積流量的模型，該模型根據(jù)Hagen-Poiseuille 方程修正得到，其表達(dá)式為

式中：kn為克努森數(shù)，kn=l/λ，定義為氣體分子的平均自由程l與孔隙直徑λ的比值，l的單位為m；Δp為毛細(xì)管中的氣體壓力差，Pa；μ為氣體黏度，Pa·s。f(kn)為滲透率修正因子，此因子只與kn有關(guān)，其表達(dá)式為

根據(jù)kn的不同，氣體在致密多孔介質(zhì)中可以分為連續(xù)介質(zhì)區(qū)、滑移區(qū)、過渡區(qū)和自由分子區(qū)等4 個(gè)流動(dòng)區(qū)域，但是，將式（9）和式（10）代入式（8）后得到的公式具有很強(qiáng)的非線性特征，不易求解，應(yīng)用價(jià)值不大［22］。為了得到廣泛適用的滲透率修正因子，使用matlab 多項(xiàng)式曲線擬合工具，對式（9）進(jìn)行分段擬合。

當(dāng)kn≤0.001和0.001＜kn≤0.100時(shí)，氣體流動(dòng)處于連續(xù)介質(zhì)區(qū)和滑移區(qū)。運(yùn)用泰勒公式將式（10）化簡后再代入式（9），然后忽略滲透率修正因子的二階和高階修正，保留一級修正。此時(shí)，f(kn)=1+4kn，據(jù)文獻(xiàn)［23］報(bào)道，該公式正是Klinkengerg 于1941 年給出的考慮滑脫效應(yīng)的氣測滲透率的數(shù)學(xué)表達(dá)式，這里的4kn即為滑移因子，但是忽略了二階和高級修正會導(dǎo)致簡化模型值小于Beskok模型值，而通過曲線擬合的方式得到的擬合公式f(kn)=1+4.8 kn其簡化模型值與Beskok 模型值比較吻合。相對誤差為0.060%～0.272%，度量擬合優(yōu)度的統(tǒng)計(jì)量r2可達(dá)0.999 4。

當(dāng)0.1 ＜kn≤10.0 和kn＞10 時(shí)，氣體處于過渡區(qū)和自由分子區(qū)。過渡區(qū)的擬合公式為f(kn)=0.8+6.2kn，相對誤差為1.3%～8.5%，統(tǒng)計(jì)量r2為0.992。自由分子區(qū)的擬合公式為f(kn)=-4.1 +6.65kn，相對誤差為0.11%～1.93%，統(tǒng)計(jì)量r2為0.999。

綜上，通過對滲透率修正因子曲線擬合而簡化得到的模型如下

通過對式（8）的擬合得到分段函數(shù)與原函數(shù)的對比（圖1），可以看出簡化后的模型與原模型擬合誤差較小，精確度較高。通過對滲透率修正因子的簡化，不僅為后續(xù)理論計(jì)算和求解提供了基礎(chǔ)，同時(shí)也有利于Beskok 模型的推廣。

圖1 Beskok 模型與簡化模型的對比Fig.1 Comparison between Beskok model and simplified model

所以一維圓管中氣體體積流量公式也可以用一個(gè)分段函數(shù)來描述

由于多孔介質(zhì)具有不規(guī)則性和復(fù)雜性，在同一多孔介質(zhì)中孔隙大小差別很大，這使得多種流動(dòng)區(qū)域在同一多孔介質(zhì)中同時(shí)存在，根據(jù)氣體流動(dòng)區(qū)域的劃分，計(jì)算氣體通過橫截面的流量公式如下［17］

式中：λk=0.1l為滑移區(qū)和過渡區(qū)的臨界直徑；λp=10l為過渡區(qū)和自由分子區(qū)的臨界直徑。

聯(lián)合式（2）、式（5）、式（11）和式（12）可得

根據(jù)菲克定律，氣體在多孔介質(zhì)中擴(kuò)散的總的摩爾流率為［15］

式中：QF為氣體摩爾流率，mol/s；At為多孔介質(zhì)總的截面積，m2；Deff為有效擴(kuò)散系數(shù)，m2/s；ΔC為多孔介質(zhì)孔道兩端的濃度差，ΔC=Δp/(RT)，其中，R為氣體常數(shù)8.314J/(mol·K) ；T為溫度，K。計(jì)算時(shí)須要將氣體摩爾流率轉(zhuǎn)化為氣體體積流量。

聯(lián)合式（5）、式（6）、式（13）和式（14），可以求出

式（16）是基于滲透率修正因子、迂曲毛細(xì)管束模型、分形理論和菲克定律求出來的多孔介質(zhì)氣體擴(kuò)散模型。此模型的每一個(gè)參數(shù)都具有物理意義，不包含經(jīng)驗(yàn)常數(shù)。式（16）表明，氣體有效擴(kuò)散系數(shù)與多孔介質(zhì)面積分形維數(shù)、迂曲分形維數(shù)、孔隙度、溫度、最大與最小孔隙直徑、氣體分子平均自由程、氣體摩爾質(zhì)量、氣體密度和黏度等均有關(guān)，并且描述了在不同流動(dòng)區(qū)域共同作用下多孔介質(zhì)中氣體擴(kuò)散的詳細(xì)機(jī)理。

通過式（16）可得知這些參數(shù)對氣體有效擴(kuò)散系數(shù)的影響，并應(yīng)用于工程實(shí)踐中的相關(guān)參數(shù)分析。

3 分析和討論

本文建立的氣體有效擴(kuò)散系數(shù)包含了2 類參數(shù)。一類是關(guān)于多孔介質(zhì)孔隙數(shù)目和結(jié)構(gòu)的參數(shù)，如多孔介質(zhì)面積分形維數(shù)、迂曲分形維數(shù)、特征長度、孔隙度等。另一類是關(guān)于氣體的性質(zhì)參數(shù)，如密度、摩爾質(zhì)量、黏度、氣體分子平均自由程等。

多孔介質(zhì)的結(jié)構(gòu)參數(shù)通過分形理論的描述和計(jì)算得出，面積分形維數(shù)由以下公式來確定

式中：DE為歐幾里德維數(shù)，在二維空間中，DE=2；在三維空間中，DE=3。該文是在二維平面上建立的模型，所以DE=2。Feng 等［24］歸納了孔隙度和面積分形維數(shù)的實(shí)驗(yàn)數(shù)據(jù)，結(jié)果顯示當(dāng)λmin/λmax≈10-3時(shí)，理論模型值與實(shí)驗(yàn)結(jié)果比較吻合。

對于迂曲分形維數(shù)，尹帥等［25］分析了常規(guī)及非常規(guī)儲層巖石的分形特征。運(yùn)用壓汞法和吸附法研究得到了具有分形特性的多孔介質(zhì)的迂曲分形維數(shù)，其中砂巖多孔介質(zhì)的迂曲分形維數(shù)為1.207～1.704，頁巖多孔介質(zhì)的迂曲分形維數(shù)為1.668～2.187。

然后將已有實(shí)驗(yàn)數(shù)據(jù)與本文提出的有效擴(kuò)散系數(shù)的理論模型進(jìn)行對比。Currie［26］將H2作為擴(kuò)散氣體，使用了多種不同的多孔介質(zhì)來進(jìn)行擴(kuò)散實(shí)驗(yàn)。須要注意的是，實(shí)驗(yàn)中一般使用體孔隙度，而理論模型中是用面孔隙度。須要借助關(guān)系式Φ=，把面孔隙度轉(zhuǎn)化為體孔隙度。其中=1+0.63ln(1/Φ)，該公式由實(shí)驗(yàn)數(shù)據(jù)擬合得到［1］。

表1 為建立的模型的部分參數(shù)：

表1 分形模型中的部分參數(shù)與數(shù)值Table 1 Some parameters and values in the fractal model

將以上數(shù)據(jù)和公式（17）代入式（16）可以得到圖2 中的函數(shù)關(guān)系，圖2 顯示氣體有效擴(kuò)散系數(shù)的預(yù)測值隨著孔隙度的增加而增加，因?yàn)殡S著孔隙度的增加，氣體擴(kuò)散的區(qū)域和路徑會增多，有效擴(kuò)散系數(shù)相應(yīng)增加。當(dāng)孔隙度超過0.65 時(shí)，模型預(yù)測值會低于實(shí)際值，因?yàn)樵撃Ｐ椭械淖畲罂紫吨睆绞褂玫氖枪潭ㄖ?。通過對比模型預(yù)測值與實(shí)驗(yàn)數(shù)據(jù)，可知該模型能夠準(zhǔn)確地計(jì)算氣體在分形多孔介質(zhì)中的有效擴(kuò)散系數(shù)。本次建立的理論模型可以用于多孔介質(zhì)微納米級孔隙情況下的單相氣體擴(kuò)散計(jì)算，可以為非常規(guī)天然氣資源開采提供關(guān)鍵參數(shù)，也可以為氣液兩相或者多相滲透與擴(kuò)散提供參考。

圖2 分形模型的氣體有效擴(kuò)散系數(shù)與文獻(xiàn)［26］實(shí)驗(yàn)數(shù)據(jù)對比Fig.2 Comparison of effective gas diffusion coefficient of fractal model with experimental data in reference

圖3 為孔道迂曲分形維數(shù)分別取1.5，1.6，1.7時(shí)，氣體有效擴(kuò)散系數(shù)隨著多孔介質(zhì)面積分形維數(shù)的變化趨勢。從圖3 可以得出，面積分形維數(shù)越大，氣體有效擴(kuò)散系數(shù)也越大。因?yàn)槊娣e分形維數(shù)增加時(shí)，孔隙數(shù)目也會增加，孔隙度也會增大，導(dǎo)致擴(kuò)散相對更加容易。當(dāng)面積分形維數(shù)一定時(shí)，迂曲分形維數(shù)增加，氣體有效擴(kuò)散系數(shù)會減小。因?yàn)橛厍中尉S數(shù)越大，表明孔道毛細(xì)管越彎曲，擴(kuò)散阻力越大，導(dǎo)致擴(kuò)散性能減小。

圖3 面積分形維數(shù)和迂曲分形維數(shù)對氣體有效擴(kuò)散系數(shù)的影響Fig.3 Effect of area fractal dimension and tortuous fractal dimension on effective gas diffusion coefficient

圖4 為考慮滑移因子的模型與不考慮滑移因子模型的對比。據(jù)文獻(xiàn)［23］報(bào)道，Klinkengerg 于1941 年發(fā)現(xiàn)了滑移現(xiàn)象，具體表現(xiàn)為氣體的表觀滲透率明顯大于絕對滲透率。從圖4 可以看出當(dāng)考慮滑移因子時(shí)，氣體的有效擴(kuò)散系數(shù)大于未考慮滑移因子的氣體有效擴(kuò)散系數(shù)，這與Klinkengerg 描述的現(xiàn)象一致。由圖4 分析可知，在孔隙直徑較大的情況下，滑移因子的影響越來越小，當(dāng)孔隙直徑較小時(shí)，滑移因子對有效擴(kuò)散系數(shù)的影響較大，這與滑脫效應(yīng)一般發(fā)生在低滲透率和低孔隙度相一致。從圖4 還可以看出，考慮了滑移因子的分形模型與不同多孔介質(zhì)的實(shí)驗(yàn)數(shù)據(jù)更加吻合，這表明考慮了滑移因子的分形模型能夠適用于更多種類的多孔介質(zhì)。

圖4 滑移因子對氣體有效擴(kuò)散系數(shù)的影響Fig.4 Effect of slip factor on effective gas diffusion coefficient

圖5 為在不同最大孔隙直徑下，有效擴(kuò)散系數(shù)隨孔隙度變化的情況。

圖5 最大孔隙直徑對氣體有效擴(kuò)散系數(shù)的影響Fig.5 Effect of maximum pore diameter on effective gas diffusion coefficient

從圖5 可以得出，當(dāng)最大孔隙直徑減小時(shí)，氣體有效擴(kuò)散系數(shù)也會減小。因?yàn)楦鶕?jù)分形理論，分形多孔介質(zhì)的面積分形維數(shù)、迂曲分形維數(shù)、最小最大孔徑比、擴(kuò)散截面面積等都與最大孔隙直徑有關(guān)，所以最大孔隙直徑是氣體有效擴(kuò)散系數(shù)的主要影響因素。

圖6 為使用不同的擴(kuò)散氣體時(shí)，有效擴(kuò)散系數(shù)隨孔隙度的變化。由圖6 可知，不同的氣體在多孔介質(zhì)中的有效擴(kuò)散系數(shù)不同，當(dāng)擴(kuò)散的氣體的摩爾質(zhì)量，或者氣體密度增加時(shí)，有效擴(kuò)散系數(shù)會減小。因?yàn)闅怏w質(zhì)量越大，擴(kuò)散速率就會越小，有效擴(kuò)散系數(shù)也會越小，這與氣體擴(kuò)散定律描述的現(xiàn)象一致。

圖6 不同氣體對氣體有效擴(kuò)散系數(shù)的影響Fig.6 Effect of different kinds of gases on effective gas diffusion coefficient

4 結(jié)論

（1）基于多孔介質(zhì)的分形特征、迂曲毛細(xì)管束模型、菲克定律、滲透率修正因子建立了多孔介質(zhì)中氣體擴(kuò)散的分形模型。該模型包含了多孔介質(zhì)面積分形維數(shù)、迂曲分形維數(shù)、孔隙度、最大孔隙直徑、氣體密度和摩爾質(zhì)量等參數(shù)。通過參數(shù)分析可以發(fā)現(xiàn)：多孔介質(zhì)孔隙度、面積分形維數(shù)、最大孔隙直徑、氣體的物理性質(zhì)等對有效擴(kuò)散系數(shù)均有顯著影響。

（2）對比分形多孔介質(zhì)氣體擴(kuò)散模型與Currie關(guān)于氣體擴(kuò)散率的數(shù)據(jù)發(fā)現(xiàn)：體孔隙度為0.10～0.65 時(shí)，模型與多種不同的多孔介質(zhì)實(shí)驗(yàn)數(shù)據(jù)均具有較好的一致性，說明該模型可以用于計(jì)算具有分形特性的多孔介質(zhì)的氣體有效擴(kuò)散系數(shù)，對非常規(guī)天然氣在實(shí)際多孔介質(zhì)中的有效擴(kuò)散系數(shù)的計(jì)算具有一定的參考意義。同時(shí)也能為后續(xù)研究多孔介質(zhì)氣液兩相或多相滲透與擴(kuò)散提供參考。

（3）考慮滑移因子的氣體有效擴(kuò)散系數(shù)大于未考慮滑移因子的氣體有效擴(kuò)散系數(shù)，在較小的孔隙直徑下，兩者的差距比較明顯，隨著孔隙直徑的增大，兩者的差異逐漸減小，直至相等?？紤]了滑移因子的氣體有效擴(kuò)散系數(shù)與Currie 關(guān)于氣體擴(kuò)散率的實(shí)驗(yàn)數(shù)據(jù)更加吻合。