微尺度通道內稀薄氣體高階努森數滲透率修正模型

盧 銀 彬

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

0 引言

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

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

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

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

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

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

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

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

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

對于無稀薄效應的氣體或者液體，其滲透率與流體物性無關，僅與通道結構相關，與巖石絕對滲透率相等。

根據公式（3）并結合式（2）得到與流體物性相關的氣體滲透率修正系數的計算式，即

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

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

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

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

表1 滑移系數C1、C2統計表

2 矩方法

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

2.1 矩方程

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

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

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

2.2 數值求解方法

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

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

2.3 邊界條件

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

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

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

3.1 R26矩方法準確性驗證

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

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

3.2 高階努森數氣體滲透率修正模型

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

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

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

式中a、b、c分別表示修正系數。

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

圖2 平板微通道中Kn與關系曲線圖

泰勒展開式（9），得

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

3.3 模型驗證

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

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

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

圖4 本文模型、Tang模型和實驗數據結果對比圖（平板微通道）

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

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

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

圖5 圓管微通道中Kn與關系曲線圖

泰勒展開式（11），有

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

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

5 結論

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

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