超臨界CO2在巖石孔隙內流動狀態的分子模擬

高 騰，趙伶玲，李偲宇

(東南大學 能源與環境學院，南京 210096)

1 引 言

深部儲層的超臨界CO2(溫度高于304 K且壓強高于73.8 atm)泄露是地質碳封存的主要風險. 當泄漏發生時，由于巖石孔隙尺寸、溫度及壓強的影響，超臨界CO2在巖石孔隙中存在復雜的密度分布和流動特性[1-3]. 因此，研究超臨界CO2在巖石納米孔隙中的密度分布和流動特性，可以提高對巖石納米孔隙中超臨界CO2流動特性的認識.

目前，國內外學者通過分子模擬的方法對納米孔隙內不同物質的流動過程進行了廣泛研究，取得了豐碩的成果. Ghorbanian等[4]通過分子模擬的方法研究了納米流動過程中，孔隙尺寸對H2O分子密度分布、流動速度及動力粘度的影響. Markesteijn等[5]采用分子模擬研究了壓力驅動流動時，不同力場參數及溫度對H2O分子流動速度的影響，驗證了可以更加準確描述H2O分子納米流動的力場模型. Wang等[6]通過分子模擬的方法分析了超臨界辛烷在納米級巖石孔隙內的流動速度分布，發現溫度和孔隙尺寸均會對超臨界辛烷的流動速度產生影響. Jiang等[7]采用分子模擬的方法研究了孔隙大小對超臨界甲烷在硅納米通道內流動過程的影響，分析了甲烷分子在巖石孔隙內的密度分布和流動狀態. 然而，對巖石納米孔隙中超臨界CO2流動狀態的研究卻罕見報道.

本文采用分子動力學(MD)模擬的方法，選取鎂橄欖石(Mg2SiO4)晶體塊為固體相，構建了Mg2SiO4-CO2-Mg2SiO4的模擬系統，分析比較了地質封存條件下，孔隙尺寸(<37 nm)、溫度(350 K-400 K)、壓強(200 atm-300 atm)對超臨界CO2在巖石孔隙內流動特性的影響，為指導巖石孔隙內超臨界CO2泄漏提供理論指導.

2 分子模擬方法

本文采用分子模擬的方法，在計算中考慮了分子間的非鍵結作用及分子內的鍵結作用. 分子間的非鍵結作用包括范德華力和庫侖靜電力，范德華力采用L-J勢能函數描述，庫侖靜電力采用庫侖定律描述；分子內的鍵結作用包括鍵拉伸和鍵角彎曲，均采用諧波勢能函數進行模擬. 模擬中Mg2SiO4和CO2分別選擇CLAYFF[8]和EPM2力場模型[9]進行描述，示于表1.

表1 Mg2SiO4和CO2力場

注：ε代表L-J勢能曲線中的勢能阱深度，nm；σ代表勢能為零時，兩點間的距離，nm.

計算選取Mg2SiO4晶體塊為固體相，構建了如圖1所示的Mg2SiO4-CO2-Mg2SiO4的模擬系統. Mg2SiO4晶體塊沿Y方向的厚度為5 nm，系統在X和Z方向采用三維周期性邊界條件，范德華作用截距設為1.2 nm. 分子模擬所采用的軟件為LAMMPS[10]，模擬采用Leap-Frog算法，以1 fs的時間步長進行求解.

本文首先采用NPT系綜進行5 nm的平衡模擬，使系統的密度和結構保持穩定. 在平衡模擬中，溫度及壓強的耦合分別采用基于Berendsen耦合器的velocity rescaling方法[11]和semi-isotropic方法[12]. 系統平衡后，采用NVT系綜進行30 ns的非平衡模擬. 在非平衡動力學模擬中，沿+Z方向對每個原子施加一定大小的的體積力，推動流體流動[13].

圖1 超臨界CO2在巖石孔隙內的流動示意圖Fig. 1 Schematic diagram of the flow of supercritical CO2

3 模擬結果與討論

3.1 模型驗證

本文首先對鎂橄欖石的晶體結構進行研究，計算了350 K和200 atm條件下Mg-O、Si-O原子對的徑向分布函數(RDF)，示于圖2. 如圖所示，Mg-O、Si-O的RDF曲線第一峰的位置對應了Mg-O、Si-O鍵的鍵長，Mg-O、Si-O原子對的第一峰分別位于r=2.01 ?、r=1.59 ?位置，這與Kudoh等[14]人在實驗中得到的Mg-O(r=2.09 ?)、Si-O(r=1.63 ?)原子對的鍵長較為符合.

圖2 350 K、200 atm條件下Mg2SiO4中Mg-O、Si-O原子對的RDF曲線Fig. 2 RDF of Mg-O and Si-O in Mg2SiO4 under 350 K and 200 atm

本文接下來對鎂橄欖石(Mg2SiO4)晶體塊進行了拉伸模擬，晶體塊受到拉力并隨著時間的推移沿受力方向發生形變，最終使Mg2SiO4晶體塊的長度達到初始長度的1.2倍. 經數據分析，計算得到Mg2SiO4晶體塊拉伸過程的應力-應變(σ-ε)曲線，示于圖3. 從圖中可以看出，Mg2SiO4先后經歷了彈性形變(ε≤0.05)、塑性形變(0.05<ε≤0.125)和斷裂(ε>0.125)三個階段. 在彈性形變階段，Mg2SiO4所受應力(σ)與其應變(ε)基本成正比. 通過擬合計算，此階段中應力與應變的比值(即彈性模量)近似為293 GPa，較好地符合Liu等[15]人在實驗中測得的彈性模量(約為287 GPa).

此外，應用所建CO2分子的模型和力場，計算得到不同溫度條件下CO2分子的密度變化，示于圖4. 從該圖可以看出，CO2的密度隨溫度的升高而降低，模擬結果與美國標準與技術研究所(NIST)[16]提供的實驗數據[17]吻合較好，進一步驗證了本文所建立的計算模型和力場具有一定的可靠性和準確性.

圖3 Mg2SiO4的應力-應變曲線Fig. 3 Stress-strain curve of Mg2SiO4

圖4 超臨界CO2密度隨溫度的變化Fig. 4 Densities of CO2 under diffferent temperatures

3.2 孔隙尺寸的影響

3.2.1孔隙尺寸對速度的影響

本文首先比較了超臨界CO2在不同大小巖石孔隙內的速度分布，示于圖5. 該圖以孔隙中間位置為0點.

從圖中可以看出，當孔隙尺寸H大于5 nm時，超臨界CO2的速度曲線呈“拋物線”狀，符合Poiseuille流動的運動規律，即由巖石表面至孔隙中心，CO2分子受到巖石表面分子的吸引力逐漸減小，阻礙CO2運動的作用力由壁面到中心逐漸減小，呈現出壁面速度低而中心速度高的分布特征. 當孔隙尺寸H小于5 nm時，超臨界CO2分子的速度在巖石孔隙內無明顯的“拋物線”狀，類似于平推流. 這與Botan等[18]研究的水在不同孔隙尺寸的粘土納米通道內的流動規律相符. 當孔隙尺寸H小于5 nm時，巖石壁面與CO2分子間的相互作用會強化CO2分子間的相互作用力，從而導致CO2分子以類似平推流的流動速度運動，而當孔隙尺寸H大于5 nm時，在遠離壁面的CO2分子間的相互作用較弱，使得孔隙中間區域的超臨界CO2分子可以在巖石孔隙內呈Poiseuille流動.

圖5 不同孔隙尺寸下超臨界CO2的速度分布Fig.5 Velocities of supercritical CO2 under different pore sizes

3.2.2孔隙尺寸對密度的影響

超臨界CO2分子的平均密度隨孔隙尺寸的變化，示于圖6. 從圖中可以看出，當孔隙尺寸H小于15 nm時，超臨界CO2的平均密度隨孔隙的減小而增大；而當孔隙尺寸H大于15 nm時，超臨界CO2的平均密度不再受孔隙尺寸的影響，趨于620 kg/m3，這也與美國標準與技術研究所(NIST)提供的實驗數據吻合較好.

圖6 不同孔隙尺寸下超臨界CO2的平均密度Fig.6 Average density of supercritical CO2 under different pore sizes

為分析這種現象產生的原因，本文比較了孔隙尺寸H為7 nm和15 nm時，超臨界CO2分子在巖石孔隙內的密度分布，示于圖7. 該圖以孔隙中間位置為0點. 從圖中可以看出，在鎂橄欖石的近壁面處，超臨界CO2分子的密度分布呈震蕩狀態，有兩個明顯的峰，這主要是因為鎂橄欖石表面與超臨界CO2分子間具有強烈的相互作用，使CO2分子在鎂橄欖石近壁面形成一種“類固體”[19]，具有較高的密度，這與Ghorbanian等[20]人發現的單相H2O分子在碳納米管中的密度分布規律相似.

進一步的分析發現：在不同大小的巖石孔隙內均發現了密度震蕩，近壁面處的密度震蕩現象并不會隨巖石孔隙尺寸的增大或減小而消失；同時，密度震蕩的距離L并不受孔隙尺寸大小的影響，在不同大小的巖石孔隙內，密度震蕩的距離均持續到距離巖石表面L=1.5 nm處.

因此，在孔隙尺寸H小于15 nm時，孔隙尺寸越小，近壁面處密度震蕩的影響作用越明顯，孔隙內超臨界CO2分子的平均密度越大. 雖然隨著孔隙尺寸的不斷增大，近壁面處的密度震蕩現象依然存在，但是由于其震蕩距離L保持不變，其對超臨界CO2平均密度的影響不斷減小，最終使得孔隙尺寸H大于15 nm時，密度震蕩對超臨界CO2平均密度的影響不再重要.

3.3 溫度和壓強的影響

3.3.1溫度和壓強對速度的影響

考慮到孔隙尺寸H小于15 nm時，密度震蕩現象的影響，本文選取孔隙尺寸H為15nm，研究不同溫度和壓強下超臨界CO2在鎂橄欖石表面的速度分布，示于圖8.

圖7 不同孔隙尺寸下超臨界CO2的密度分布Fig.7 Density profile of supercritical CO2 under different pore sizes

圖8 不同溫度和壓強下超臨界CO2的速度變化Fig.8 Velocity profile of supercritical CO2 under different temperatures and pressures

從圖中可以看出，隨著溫度的升高和壓強的降低，巖石孔隙內超臨界CO2分子的速度逐漸增加，其速度曲線仍然呈“拋物線”狀，但是速度分布的“拋物線”曲率卻呈現出不同. 這種行為在宏觀上反映為流體的動力粘度變化，即隨著溫度的升高，CO2分子的流動速度不斷增加，動力粘度不斷減小；隨著壓強升高，CO2的流動速度不斷減小，動力粘度卻不斷增加. 而Bao等[21]人發現：H2O在納米孔隙中的動力粘度隨溫度的升高而增加，卻不受壓強直接影響.

因此，超臨界CO2在鎂橄欖石孔隙內流動時，采用降低溫度、增加壓強的方法可以減弱鎂橄欖石表面對CO2分子的吸引力，使得CO2分子的流動速度降低，防止超臨界CO2的泄漏.

3.3.2溫度和壓強對密度的影響

孔隙尺寸H為15 nm時，不同溫度和壓強下超臨界CO2在鎂橄欖石表面的密度分布，示于圖9. 從圖中可以看出，不同溫度和壓強條件下，超臨界CO2分子在鎂橄欖石表面的密度震蕩依然存在，而且兩個峰的峰值對應的位置不受溫度和壓強的影響. 其中，第一個峰的峰值相對較大，決定于巖石表面與超臨界CO2分子間的相互作用力，受溫度和壓強影響影響較大. 隨溫度的升高而降低，隨壓強的增大而增加，即增大壓強可導致鎂橄欖石表面對超臨界CO2分子的吸引力增強，而升高溫度反而會減弱這種吸引力，增加CO2分子的流動速度.

圖9 不同溫度和壓強下近壁面處超臨界CO2的密度分布Fig.9 Density profile of supercritical CO2 near the surface under different temperatures and pressures

4 結 論

本文采用非平衡動力學模擬的方法研究了超臨界CO2在巖石孔隙內的流動特性，分析了孔隙大小、溫度及壓強對CO2流動特性的影響，得到以下結論.

(1) 超臨界CO2在巖石孔隙內流動時，當巖石孔隙小于5.0 nm時，速度分布呈平推流的運動規律，只有當巖石孔隙大于5.0 nm時，超臨界CO2分子在巖石孔隙內的流動才符合Poiseuille流動的運動規律.

(2) 超臨界CO2在巖石近壁面處具有較大的分子密度，呈現一種密度震蕩現象，單側壁面的震蕩距離在1.5 nm左右，不受孔隙尺寸大小的影響，當巖石孔隙小于15 nm時，這種密度震蕩就會顯著影響超臨界CO2在巖石孔隙內的平均密度.

(3)巖石壁面處CO2分子的密度震蕩受溫度和壓強的影響，采用降低溫度和增大降低壓強的方法，可以使得密度震蕩的第一峰值增大，造成巖石表面對CO2分子的吸引力增強，減小CO2分子的流動速度.