多孔介質內甲烷碳同位素遷移擴散的分子動力學模擬研究

劉 暢，張慧，趙耀中，陳威，2

(1.蘇州大學能源學院，江蘇 蘇州 215006；2.中國航天科技集團有限公司氫能工程技術研發中心，北京 100074；3.北京航天試驗技術研究所，北京 100074)

隨著常規油氣資源的日趨緊缺和出于能源清潔的考慮，頁巖氣這一清潔能源已成為多國學者研究的熱點之一[1]。頁巖氣作為一種非常規能源，由于其具有儲量豐富，燃燒時污染程度低的優良性能[2]，在能源領域發揮著越來越重要的作用。近年來，水平鉆井和壓裂技術的發展提高了頁巖氣的開采能力，促進了美國“頁巖氣革命”的成功。學者們對致密頁巖內的氣體流動展開了相關研究，但是目前對于頁巖氣在納米級孔尤其是在含孔喉納米孔中的吸附、擴散和同位素分餾特征仍不十分清晰[3]。

穩定碳同位素氣體的組成、分餾可以為各成熟階段的天然氣提供有價值的地球化學信息，如氣體的成因、示蹤運移、聚集和散失等[4-7]，也被作為探索頁巖內部氣體分子微觀尺度非線性輸運過程的新方法。甲烷氣體產出具有明顯的多尺度流動的特征，包含了滲流、吸附/脫附和努森擴散等大量分子層次的非線性流動，其流體的滲透率和有效擴散系數等參數不再是常數[8]。在納米級孔內不但存在烷烴氣體分子的相對運動，而且存在甲烷同位素氣體(12CH4和13CH4)的分餾(即在流動、吸附和熱解等物理或者化學過程中，輕、重同位素分子的豐度在物相之間重新分配，12CH4和13CH4相對比例發生浮動)。由于實驗技術和測量手段的限制，目前還無法直接觀測熱采工況下致密頁巖內部甲烷的賦存狀態和微觀流動現象。因此，模擬成為研究同位素微觀狀態特征和擴散形式的有效手段。ZHANG等人使用分子動力學模擬的方法研究了頁巖氣在孔喉處的吸附特性，并且指出了頁巖氣藏中孔喉的下限[9]；NADER等人利用晶格玻爾茲曼模型研究了喉-孔比、Kn數對頁巖氣流動特性和滑移速率的影響[10]；一般認為甲烷碳同位素在受限空間中的分子擴散效應和物理吸附/脫附是造成其分餾的兩個主要原因。XIA等人通過建立多孔流動模型研究吸附/脫附過程對同位素分餾造成的影響，LIU等人通過建立數學模型比較了烷烴碳同位素的分餾程度和吸附氣比例的關系[11]。研究頁巖氣甲烷碳同位素的吸附和擴散過程能夠強化對其流動特性的理解，有助獲得更多頁巖氣的流動和儲存信息；探究頁巖氣甲烷碳同位素在頁巖氣藏中的分餾程度對追蹤頁巖氣的生成和演化、評價頁巖氣的相對吸附能力具有重要的意義。

因此，本文借助分子動力學的方法研究了甲烷同位素氣體在頁巖納米級孔中(包含孔喉)的吸/脫附和擴散現象，獲得不同工況下12CH4和13CH4擴散能力和分餾特征，旨在分析不同孔結構(包括孔喉大小)對甲烷同位素分餾的影響規律。

1 分子動力學理論基礎與模型

1.1 分子動力學理論

本文應用分子動力學模擬軟件(LAMMPS)研究甲烷碳同位素在有無孔喉存在的納米級有機孔中的擴散速率和分餾特征。模擬甲烷碳同位素在壓差的驅動下沿著碳納米管流動。粒子間相互作用采用截斷及移位的Lennard-Jones(LJ)勢能函數[12]，即：

(1)

公式(1)中，E為分子勢能，r為分子間距，rc為截斷半徑，本文中的rc=10 ?，σ是相互作用的勢能正好為零的兩體距離，ε是勢能阱的深度。相同粒子間的相互作用勢參數σ和ε[12]見表1，不同粒子間的σ和ε參數由Lorentz-Berthelot混合規則確定[13]：

表1 相同粒子間的相互作用勢參數

(2)

1.2 同位素分餾模型

分餾模型采用兩種不同結構的碳納米管，用來模擬甲烷碳同位素的分餾。如圖1(a)是含有孔喉存在的碳納米管，其長度為20.4 nm，直徑為4 nm，孔喉的長度為0.4 nm，直徑為1 nm，孔喉左右兩側碳納米管長度相同，均為10 nm。為了對比孔喉的影響，構建相同尺寸的不含有孔喉的碳納米管用來研究甲烷同位素的分餾，即：碳納米管的長度為20.4 nm，直徑為4 nm(圖1(b))。

圖1 有/無孔喉存在的碳納米管模型

為了比較甲烷同位素在不同孔內的分餾特征，12CH4和13CH4同位素在兩種管內的擴散過程被分別模擬。首先使用巨正則蒙特卡羅方法，模擬了頁巖在地下貯存353 K左右的溫度和20 MPa的壓力條件，使得模擬結果更加接近實際情況。在溫度為353 K、20 MPa的壓強下在含有孔喉的碳納米管中生成一定數目的12CH4分子。模擬初始階段，先通過Nose-Hoover控溫器使管內氣體分子在NVT系綜下維持溫度為353 K，經過時長為2 ns的弛豫時間后粒子溫度、壓力等熱力學量到達穩定，且在一定時間內不發生波動，則系統達到熱平衡，然后給每個粒子在流動方向施加外力，達到模擬孔隙內甲烷氣體壓差作用下流動的目的；13CH4粒子在同樣的條件下生成并在相同壓差條件下流動，分別統計并計算13CH4和12CH4的遷移擴散系數進行分析。本文考慮了非理想氣體效應，利用Peng-Robinson狀態方程求出逸度系數來表征實際氣體偏離理想氣體的程度[14]。在不含孔喉的情況下，以同樣的方法生成粒子并附加相同的條件，分別統計并計算13CH4和12CH4的遷移擴散系數，與含有孔喉的情況進行對比分析，旨在探究孔喉的存在對甲烷碳同位素分子分餾特性的影響。

2 氣體流動特性和擴散系數的計算

2.1 氣體流動特性

與傳統宏觀的連續流體流動狀態不同，由于高努森數效應的存在，達西方程已經不適用于描述甲烷同位素氣體在納米孔中的流動。頁巖氣的流動狀態可由努森數來判定。努森數(Kn)的定義為[15]：

(3)

(4)

(5)

公式(3)中，λ為平均分子自由程，H為流動通道的特征參數。公式(4)中，R為氣體常數，T為絕對溫度，Mw為氣體分子的摩爾質量，μ為粘度系數，P為氣體壓強。公式(5)中，KB為玻爾茲曼常數，dm為分子直徑。當Kn<10-3時，流動屬于連續流；當10-310時，流動形式為分子自由擴散[16]。本文模擬的Kn數的范圍為0.06～0.18，屬于滑脫流和過渡流的范圍，與傳統無滑移的連續流不同，需要考慮壁面的滑移效應。

2.2 擴散系數

按流動方向(z軸)統計區域均分成長度為4 ?的小區域，每個小區域通過的摩爾流率計算公式如下：

(6)

公式(6)中，A為碳納米管的截面積，L為小區域的長度(此處為4 ?)，Ni為小區域內甲烷分子的個數。當系統達到穩定之后，由于邊緣效應的存在，進口區域和出口區域的摩爾流率波動很大，但中間區域摩爾流率較穩定，因此選取穩定區域的摩爾流率平均值作為頁巖有機納米孔的平均摩爾流率，從而可求擴散系數[17]：

(7)

公式(7)中，D為12CH4的擴散系數，D′為13CH4的擴散系數。其中R為氣體常數，T為溫度，Lpore為碳納米管的長度，F和F′分別為12CH4和13CH4分子的摩爾流率，ε為孔隙度，A為碳納米管的截面積，ΔP是通過碳納米管的壓差，壓差與外力的關系由下式給出：

(8)

公式(8)中，n為外力驅動的粒子數目，Fex為外加驅動力。由于本文中的驅動力以壓差的形式表示，孔喉兩端施加的10 MPa壓力差為粒子擴散的驅動力。因此，12CH4和13CH4的遷移擴散系數的公式可以表示如下：

(9)

3 結果和討論

3.1 甲烷同位素分子在孔內的徑向密度分布

為了直觀地反映納米級孔內甲烷碳同位素分子的吸附特性和密度分布，圖2給出了甲烷同位素分子在無孔喉存在的納米通道的遷移擴散現象，圖3給出了甲烷同位素分子在無孔喉模型中徑向分布函數收斂后的相關的數據(壓差為10 MPa)。甲烷同位素分子在靠近壁面處形成第一和第二吸附層，局部密度遠高于孔中間分子游離區的密度。由于壁面與甲烷碳同位素分子之間存在較強范德華力的作用，粒子易在吸附位點取得平衡并形成吸附層。并且由于甲烷碳同位素分子數量多，在壁面處形成競爭吸附，由此形成第一吸附層和第二吸附層，靠近壁面的第一吸附層甲烷碳同位素分子的密度遠大于第二吸附層甲烷碳同位素分子的密度。而當甲烷碳同位素遠離壁面時，壁面對分子的吸附作用減弱，而分子間的作用力占主導。

圖2 無孔喉模型中CH4在10 MPa壓差下沿著Z軸方向的流動

圖3 無孔喉模型中12CH4和13CH4在10 MPa壓差下的徑向密度分布

圖4為甲烷同位素在含有孔喉的納米孔中的流動模型。由圖4可知，壁面處13CH4分子的密度更高，這是由于12CH4和13CH4在相同的溫度條件下，質量大的13CH4粒子處于更深的勢阱，更容易被壁面所束縛。

圖4 甲烷同位素在含有孔喉的納米孔中的流動模型

孔喉的存在通常會引起甲烷同位素密度分布在其前后產生差異。甲烷同位素在10 MPa的壓差下沿著Z軸正方向遷移擴散。圖5為甲烷通過孔喉前后的徑向密度分布。從圖5中可以看出，流過孔喉之前甲烷同位素的密度分布明顯高于流過孔喉之后的密度分布，這是由于孔喉對甲烷碳同位素的通過具有一定的選擇性，只有帶有足夠能量的分子才能越過孔喉的能量壁壘，進入到孔喉的另一側[18]。而且由于孔喉處直徑較小，一部分粒子會與孔喉附近壁面相撞而改變運動路徑，從而降低了自身通過孔喉的概率。

圖5 甲烷通過孔喉前后的徑向密度分布

3.2 甲烷同位素分子在孔內的速度分布

在本文模擬中，甲烷碳同位素分子(12CH4和13CH4)在納米通道中的流動形式處于滑脫流和過渡流兩種流態(0.06

此外，對比圖6(a)、圖6(b)可知，甲烷碳同位素分子在含孔喉有機孔中流動的徑向分布速度小于無孔喉有機孔中的速度。這是由于氣體分子會與流動方向上的壁面發生碰撞，粒子存在一定的概率不能通過孔喉。而且由于孔喉部分直徑小，粒子可以自由移動的區域減小，由此阻滯了甲烷碳同位素分子的遷移擴散過程，從而降低了徑向分布速度。

圖6 12CH4和13CH4分子在有無孔喉存在的碳納米管內流動的徑向速度分布

3.3 同位素遷移擴散系數和分餾程度

頁巖氣隨著生產過程的進行，孔內壓力逐漸降低，由于甲烷同位素分子脫附能力的差異導致其分餾程度發生變化[19]。通過改變納米孔中對甲烷的加力大小，可以模擬甲烷同位素在不同壓差下的遷移擴散。圖7(a)、7(b)分別表示在有孔喉和無孔喉頁巖有機孔中的甲烷碳同位素分子遷移擴散系數隨壓差的變化。從圖中可以看出，隨著壓差降低，12CH4和13CH4的遷移擴散系數均減小，并且有孔喉存在時，由于孔喉的阻礙和內部的吸附作用，甲烷碳同位素分子的遷移擴散系數更小，即遷移擴散速度減慢。此外，在任一壓差下12CH4的遷移擴散系數均大于13CH4，這是由于13CH4分子的質量大，碳納米管對其吸附能力更強；此外，受驅動力影響，質量大的13CH4的擴散速度更慢，由此造成同位素分餾。

圖7 甲烷碳同位素分子遷移擴散系數隨壓差的變化

吸附和擴散是導致同位素分餾兩個主要原因，本文采用13CH4和12CH4擴散系數的比值δ表征同位素分餾程度的大小。圖8為甲烷碳同位素分子的分餾程度隨壓差的變化。由圖8可知，在無孔喉有機納米孔中，12CH4和13CH4存在分餾現象，并且分餾程度小于同位素在含孔喉有機孔中的分餾程度。本文中孔喉直徑(1 nm)與甲烷分子直徑(0.38 nm)比較接近，孔喉對甲烷分子的吸附作用明顯較大，兩種同位素分子可以自由移動的區域較小。因此當兩種粒子進入孔喉時，受范德華力的影響較大，并且由于13CH4的質量大于12CH4，所處的勢阱位置更深，更易吸附在壁面，所以與無孔喉有機孔相比，在孔喉納米孔中甲烷碳同位素分子的分餾程度更大。

圖8 甲烷碳同位素分子的分餾程度隨壓差的變化

甲烷碳同位素氣體在有無孔喉的納米孔中的分餾程度均隨壓差的降低先增大后減小，即甲烷碳同位素分子先變 “輕” 后變 “重”，模擬結果與LI等人獲得的實驗結果變化趨勢一致[20]。這是因為隨著壓差降低，兩者中吸附能力較弱的12CH4首先從壁面大量脫附，而13CH4的脫附程度較弱，由此造成兩種同位素分子分餾程度增大。繼續降低壓差，吸附能力更強的13CH4在壓差的作用下開始從壁面大量脫附，因此導致兩種分子的分餾程度逐漸減弱[21]。

4 結語

本文應用分子動力學模擬方法研究了353 K和不同壓差(10～30 MPa)下12CH4和13CH4分別在包含和不包含孔喉的頁巖有機納米孔中的徑向密度分布以及摩爾流率，分析了孔喉結構對甲烷碳同位素(12CH4和13CH4)分餾程度的影響，得出以下結論：

(1)當有機孔不含孔喉時，甲烷在孔壁面的吸附和擴散作用會引起甲烷碳同位素分餾。在相同壓差下，12CH4的質量更小，因此12CH4有更高的擴散速度，同時壁面對13CH4更強的吸附能力進一步減弱13CH4的擴散速度，導致了碳同位素氣體的分餾。

(2)當有機孔含有孔喉時，由于孔喉的阻礙作用，氣體分子會與流動方向上的壁面發生碰撞，降低通過孔喉的概率，從而導致有孔喉存在時同位素的徑向分布速度變小。

(3)由于孔喉直徑小，使得含孔喉的納米孔中絕大部分區域屬于吸附區，因此進入孔喉的粒子受壁面吸附影響，而13CH4的質量大于12CH4的質量，所處的勢阱更深，導致孔喉對13CH4的截留作用更強，進而增加了同位素分餾程度。

(4)由于同位素吸附/脫附和擴散能力的差異，隨著壓差的降低，頁巖氣甲烷碳同位素分餾程度在不同結構的納米孔內均為先增大后減小。