微尺度通道內(nèi)稀薄氣體高階努森數(shù)滲透率修正模型

盧 銀 彬

西安石油大學機械工程學院

0 引言

頁巖氣儲集在低孔隙度、特低滲透率暗色泥頁巖或高碳泥頁巖層系中[1-2]，這些層系孔隙的尺寸主要屬于微納米級。已有的研究表明，當氣體流經(jīng)常規(guī)通道時，由于流動通道的特征尺寸遠大于氣體分子運動的平均自由程，可視氣體為連續(xù)介質(zhì)，采用經(jīng)典Hagen-Poiseuille理論公式計算氣體流量[3]。但是，當氣體在微納米級孔隙通道中流動時，由于通道尺寸接近甚至小于氣體分子運動的平均自由程，此時氣體流動會產(chǎn)生稀薄效應，宏觀上表現(xiàn)為相同壓降梯度下獲得的氣體流量大于Hagen-Poiseuille理論公式計算的流量，連續(xù)性假設(shè)不再適用，傳統(tǒng)Navier-Stokes方程失效。劉杰等[4]認為在頁巖納米級孔隙中氣體流態(tài)為滑脫流或過渡流，不存在連續(xù)流，且孔隙越小，壓力越低，滑脫流越弱，氣體稀薄效應越強。定量研究頁巖氣流動的稀薄效應，對頁巖氣井產(chǎn)能的精確評價有重要意義[5]。目前已有較多一階和二階滑移模型用于描述氣體稀薄效應，這些模型通過采用格子Boltzmann方法（LBM）、直接模擬蒙特卡洛法（Direct Simulation Monte Carlo，簡稱DSMC法）或者實驗等方法求得模型中的滑移系數(shù)，但所建立的模型均缺乏普適性。為此，筆者首先采用R26矩方法對平板微通道中氣體的流動進行數(shù)值模擬，并與DSMC法、R13矩方法的模擬結(jié)果進行對比，然后基于R26矩方法的模擬結(jié)果建立平板微通道與圓管微通道的高階努森數(shù)氣體滲透率修正模型，計算不同努森數(shù)對應的氣體滲透率修正系數(shù)，并與Tang模型預測結(jié)果、實驗數(shù)據(jù)及線性Boltzmann方程解進行對比。所建立的模型預測精度高且具有普適性。

1 含一階和二階滑移系數(shù)的氣體滲透率修正模型

1879年，Maxwell[6]通過理論分析研究，首次提出氣體流動時在壁面存在滑移效應。而后，Knudsen[7]于1909年提出采用努森數(shù)（Kn）來描述氣體稀薄效應，Kn的計算式為：

式中λ表示分子平均自由程，m；R表示流動通道特征尺寸，m，對于平板通道為平板間距，對于圓管通道為圓管直徑，m；μ表示氣體動力黏度，Pa·s；p表示氣體平均壓力，Pa；Rg表示氣體常數(shù)，J/（mol·K）；T表示溫度，K；M表示氣體摩爾質(zhì)量，g/mol。

Kn越大，氣體稀薄效應越顯著。根據(jù)Kn的取值范圍劃分氣體流動區(qū)間：當Kn＜10-3時，為無滑移區(qū)，即連續(xù)區(qū)；當10-3≤Kn≤10-1時，為滑移區(qū)；當10-1＜Kn≤10時，為過渡區(qū)；當Kn＞10時，為自由分子區(qū)。

1999年，Beskok和Karniadakis[8]推導出稀薄氣體在毛細管中的流量計算式，即

式中Q表示氣體體積流量，m3/s； 表示稀薄氣體系數(shù)，無量綱，在滑移區(qū)取值為-1；r表示圓管半徑，m；Δp表示流體在毛細管進出口處的壓降，Pa；Lt表示毛細管總長度，m。

當氣體在微納米尺度通道中流動時存在稀薄效應，采用含努森數(shù)的函數(shù)對氣體在通道內(nèi)的流量（或滲透率）進行修正，即在相同壓力梯度下氣體實際流量與采用Hagen-Poiseuille理論公式計算流量的比值（或氣體表觀滲透率與巖石絕對滲透率的比值）可用含努森數(shù)的函數(shù)f（Kn）表示，即

式中Ka表示表觀滲透率，m2；K∞表示巖石絕對滲透率，m2；Q∞表示無稀薄效應時氣體體積流量，m3/s。

對于無稀薄效應的氣體或者液體，其滲透率與流體物性無關(guān)，僅與通道結(jié)構(gòu)相關(guān)，與巖石絕對滲透率相等。

根據(jù)公式（3）并結(jié)合式（2）得到與流體物性相關(guān)的氣體滲透率修正系數(shù)的計算式，即

2005年，Tang等[9]將二階滑移邊界條件應用于Navier-Stokes方程中，獲得稀薄氣體在平板通道中的滲透率修正模型，即

式中C1和C2分別表示一階、二階滑移系數(shù)。

同時，Tang等[9]將該修正模型推導至圓管通道，獲得圓管中氣體滲透率修正模型，即

表1中列出了幾組常見的滑移系數(shù)C1和C2。

表1 滑移系數(shù)C1、C2統(tǒng)計表

2 矩方法

矩方程方法（簡稱矩方法）是一種介于傳統(tǒng)流體力學計算方法和粒子算法之間的方法。1949年，Grad[17]首次推導出包含13個矩變量的Grad十三階矩方法（簡稱G13矩方法），其后Struchtrup[18]發(fā)展了正則化十三階矩方法（簡稱R13矩方法），R26矩方法[19]是基于R13矩方法獲得的。2009年，Gu和Emerson[20]采用R26矩方法對非平衡氣體流動傳熱問題進行數(shù)值模擬，指出該方法能夠準確捕捉到氣體在滑移和早期過渡區(qū)的稀薄行為，且計算量遠小于粒子方法（如DSMC法），是一種值得采用的研究方法。

2.1 矩方程

當氣體在通道中流動而速度、溫度等參數(shù)不隨時間發(fā)生改變時，可視氣體流動已達到平衡狀態(tài)。當氣體流動處于平衡狀態(tài)或者接近于平衡狀態(tài)時，對氣體在流動相空間分子分布函數(shù)（g）的5個矩變量（密度、溫度以及三維空間x、y、z方向上的速度）建立質(zhì)量、動量和能量守恒方程。當氣體流動偏離平衡狀態(tài)時，稀薄效應凸顯，需要在矩方程中引入更高階的矩變量來描述氣體流動。矩方法中，Grad[19]通過引入麥克斯韋分布函數(shù)（gM）[21]，采用Hermite多項式將分子分布函數(shù)展開，即

式中g(shù)表示分子分布函數(shù)；gM表示麥克斯韋分布函數(shù)；N表示展開后的總項數(shù)；an表示第n項矩的線性組合系數(shù)；Hn表示第n項Hermite多項式。

當N趨于∞，采用Hermite多項式可準確重構(gòu)分子分布函數(shù)。而實際計算時，無法求解無窮項，必須根據(jù)實際問題，確定需要展開的項數(shù)后截斷g函數(shù)。在Grad的理論中，截斷后的g函數(shù)記作gG,其中包含的矩被稱為“Grad矩集合”，即Grad Moment Manifold（簡稱GMM）[19]，這些矩的控制方程可通過Boltzmann方程推導獲得，簡稱為矩方程。采用四階截斷分布函數(shù)可以獲得26個矩變量構(gòu)成的方程組，通過增加矩變量數(shù)量可以提高氣體在滑移和早期過渡區(qū)流動行為的描述精度。

2.2 數(shù)值求解方法

為了在傳統(tǒng)流體力學方程基礎(chǔ)上求解矩方程，需要根據(jù)矩方程中非梯度輸運項與梯度輸運項的特點修改原始矩變量，從而使矩方程能描述對流擴散過程，經(jīng)改寫后的方程組可用于模擬氣體低速流動下的非平衡行為。

采用有限體積法數(shù)值求解改寫后的矩方程，應用中心差分格式離散擴散項和源項、CUBISTA格式離散矩方程中的對流項[22]，同時利用交錯網(wǎng)格上的SIMPLE算法耦合速度與壓力[23-24]，使用多塊網(wǎng)格技術(shù)生成復雜區(qū)域中的計算網(wǎng)格。如需了解更多細節(jié)，可參閱本文參考文獻[20,25]，此不贅述。

2.3 邊界條件

模擬計算氣體流動時需要限定氣體流動區(qū)域，為構(gòu)造R26矩方程的邊界條件，利用分子分布函數(shù)在五階截斷的近似表達式[19]，結(jié)合麥克斯韋動理論邊界條件[26]，獲得描述矩變量[20,25]的動理論模型。

3 平板微通道中氣體稀薄效應

通過R26矩方法模擬氣體在平板微通道中的流動行為。兩平行平板長度（L）均為100 μm，板間間距（h）為1 μm。模擬氣體為氮氣，溫度（T）為300 K。同時，將h作為特征尺寸，采用Kn描述氣體的非平衡狀態(tài)。對于平板微通道，K∞為h2的1/12[26]，Ka可由Darcy定律計算獲得。

3.1 R26矩方法準確性驗證

DSMC法從微觀角度出發(fā)，考慮分子之間的作用力以及分子與壁面之間的碰撞作用，通過模擬大量分子的運動行為，統(tǒng)計獲得氣體的宏觀物性以及運動等參數(shù)[27-28]。通常認為DSMC法的模擬結(jié)果非常準確，本文參考文獻[29]采用DSMC法，模擬Kn為0.038 6、0.178 5和0.537 1下氣體在平板微通道中的流動。為驗證R26矩方法的準確性，將模擬結(jié)果與本文參考文獻[29-30]的計算數(shù)據(jù)進行對比（圖1），圖1中縱坐標為，其中u表示流體流速，umax表示兩平板中心處的氣體最大流速；橫坐標為 ，Y表示距離平板中心線的距離?？梢钥闯?，R26矩方法的模擬結(jié)果與DSMC計算數(shù)據(jù)能夠很好地吻合，而R13矩方法的模擬結(jié)果[30]在Kn較大（取值為0.537 1）時，與DSMC法的模擬結(jié)果存在較大偏差，準確性較R26矩方法偏低，這是因為R26矩方法能夠較準確地捕捉到緊靠固體壁面克努森層的氣體滑移特性[31]?？梢姡捎肦26矩方法模擬研究氣體稀薄效應具有較高的計算精度。

圖1 R26矩方法、R13矩方法與DSMC法模擬結(jié)果對比圖

3.2 高階努森數(shù)氣體滲透率修正模型

目前大部分氣體表觀滲透率修正模型多為一階和二階修正，且需要選擇合適的滑移系數(shù)代入模型，由于提出的系數(shù)眾多，增加了選擇難度。2007年，Zhu等[32]首次提出高階努森數(shù)氣體滲透率修正模型，但未給出式中常數(shù)A的具體數(shù)值及參數(shù)α的函數(shù)式，模型計算式為：

結(jié)合分析Tang等[9]提出氣體滲透率修正模型，如式（5）、（6）所示，筆者認為該模型中的參數(shù)A應取值為6（針對平板）和8（針對圓管），但仍需對指數(shù)進行確定。

筆者提出平板微通道中高階努森數(shù)氣體滲透率修正模型表達式為：

式中a、b、c分別表示修正系數(shù)。

為求取式（9）中的3個參數(shù)a、b和c，采用R26矩方法模擬Kn介于0.01～1.00時氣體在平板微通道中的流動，計算得到氣體表觀滲透率，結(jié)合式（3）求得氣體滲透率修正系數(shù)，將數(shù)據(jù)點擬合，得到式（9）中a為3.94、b為0.333和c為-3.94（圖2）。

圖2 平板微通道中Kn與關(guān)系曲線圖

泰勒展開式（9），得

可見，當泰勒展開后的高階努森數(shù)氣體滲透率修正模型截斷至Kn時，即為一階努森數(shù)氣體滲透率修正模型，當截斷至Kn2時，即為二階努森數(shù)氣體滲透率修正模型。泰勒展開后項數(shù)越多，模型預測結(jié)果越準確。由此可見，本文提出高階努森數(shù)氣體滲透率修正模型具有兩大優(yōu)點：①預測精度高；②具有普適性，可以直接使用。

3.3 模型驗證

采用本文參考文獻[33-34]中的實驗數(shù)據(jù)以及線性Boltzmann方程解[35]驗證模型的有效性，如圖3所示，線性Boltzmann方程解具有較高的計算精度，本文模型能夠較好地預測實驗結(jié)果，并且與線性Boltzmann方程解非常吻合。將表1中的滑移系數(shù)[10-16]代入Tang模型[9]，如式（5）所示，可以看出，Tang模型[9]采用不同滑移系數(shù)進行預測，多組預測結(jié)果與實驗結(jié)果存在較大偏差，雖然采用Kim和Pitsch[16]提出的參數(shù)預測的結(jié)果與實驗結(jié)果較吻合，但仍與線性Boltzmann方程解[35]存在一定差距。

圖3 基于本文模型、考慮不同滑移系數(shù)的Tang模型、實驗數(shù)據(jù)和線性Boltzmann方程解的Kn與對比圖（平板微通道）

為進一步驗證模型的準確性以及確認Kim和Pitsch[16]提出系數(shù)的準確性，再次將本文模型預測結(jié)果與本文參考文獻[36]的實驗數(shù)據(jù)進行對比，可以看出，本文模型預測結(jié)果與實驗結(jié)果總體較吻合，但是Kim和Pitsch提出的系數(shù)應用于Tang模型[9]在努森數(shù)較大時存在較大偏差（圖4）。而且，有學者指出應針對不同氣體選擇不同的系數(shù)C1與C2[37]，這更增加了選擇難度。若對滑移系數(shù)感興趣，可參閱本文參考文獻[38]。

圖4 本文模型、Tang模型和實驗數(shù)據(jù)結(jié)果對比圖（平板微通道）

4 圓管微通道中氣體稀薄效應

下面模擬N2在直徑為1 μm的毛細管中的流動情況，氣體努森數(shù)介于0.01～1.00。

根據(jù)模擬結(jié)果擬合獲得高階努森數(shù)氣體滲透率修正模型（圖5），具體表達式為：

圖5 圓管微通道中Kn與關(guān)系曲線圖

泰勒展開式（11），有

將本文模型預測結(jié)果和線性Boltzmann方程解[39]進行對比（圖6），可以看出，圓管中本文模型與線性Boltzmann方程解能夠很好地吻合，證實了本文模型的準確性。另外，不同滑移系數(shù)[10-16]應用到Tang模型后預測結(jié)果差異大，且大多數(shù)預測結(jié)果與線性Boltzmann方程解存在較大偏差。

圖6 本文模型、Tang模型和線性Boltzmann方程解結(jié)果對比圖（圓管微通道）

5 結(jié)論

1）采用R26矩方法描述氣體稀薄效應，其模擬結(jié)果與DSMC計算數(shù)據(jù)吻合情況良好，且計算精度高于R13矩方法。

2）針對平板微通道和圓管微通道分別建立高階努森數(shù)氣體滲透率修正模型，前者預測結(jié)果與實驗結(jié)果、線性Boltzmann方程解吻合情況好，后者預測結(jié)果也與線性Boltzmann方程解吻合情況好，所建立的模型預測精度高且具有普適性。