盧 銀 彬

西安石油大學(xué)機(jī)械工程學(xué)院

0 引言

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

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

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

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

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

1999年，Beskok和Karniadakis[8]推導(dǎo)出稀薄氣體在毛細(xì)管中的流量計(jì)算式，即

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

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

式中Ka表示表觀滲透率，m2；K∞表示巖石絕對(duì)滲透率，m2；Q∞表示無稀薄效應(yīng)時(shí)氣體體積流量，m3/s。

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

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

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

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

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

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

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

2 矩方法

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

2.1 矩方程

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

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

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

2.2 數(shù)值求解方法

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

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

2.3 邊界條件

模擬計(jì)算氣體流動(dòng)時(shí)需要限定氣體流動(dòng)區(qū)域，為構(gòu)造R26矩方程的邊界條件，利用分子分布函數(shù)在五階截?cái)嗟慕票磉_(dá)式[19]，結(jié)合麥克斯韋動(dòng)理論邊界條件[26]，獲得描述矩變量[20,25]的動(dòng)理論模型。

3 平板微通道中氣體稀薄效應(yīng)

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

3.1 R26矩方法準(zhǔn)確性驗(yàn)證

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

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

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

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

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

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

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

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

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

泰勒展開式（9），得

可見，當(dāng)泰勒展開后的高階努森數(shù)氣體滲透率修正模型截?cái)嘀罧n時(shí)，即為一階努森數(shù)氣體滲透率修正模型，當(dāng)截?cái)嘀罧n2時(shí)，即為二階努森數(shù)氣體滲透率修正模型。泰勒展開后項(xiàng)數(shù)越多，模型預(yù)測(cè)結(jié)果越準(zhǔn)確。由此可見，本文提出高階努森數(shù)氣體滲透率修正模型具有兩大優(yōu)點(diǎn)：①預(yù)測(cè)精度高；②具有普適性，可以直接使用。

3.3 模型驗(yàn)證

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

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

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

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

4 圓管微通道中氣體稀薄效應(yīng)

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

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

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

泰勒展開式（11），有

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

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

5 結(jié)論

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

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