高 騰，趙伶玲，李偲宇

(東南大學(xué) 能源與環(huán)境學(xué)院，南京 210096)

1 引 言

深部?jī)?chǔ)層的超臨界CO2(溫度高于304 K且壓強(qiáng)高于73.8 atm)泄露是地質(zhì)碳封存的主要風(fēng)險(xiǎn). 當(dāng)泄漏發(fā)生時(shí)，由于巖石孔隙尺寸、溫度及壓強(qiáng)的影響，超臨界CO2在巖石孔隙中存在復(fù)雜的密度分布和流動(dòng)特性[1-3]. 因此，研究超臨界CO2在巖石納米孔隙中的密度分布和流動(dòng)特性，可以提高對(duì)巖石納米孔隙中超臨界CO2流動(dòng)特性的認(rèn)識(shí).

目前，國(guó)內(nèi)外學(xué)者通過分子模擬的方法對(duì)納米孔隙內(nèi)不同物質(zhì)的流動(dòng)過程進(jìn)行了廣泛研究，取得了豐碩的成果. Ghorbanian等[4]通過分子模擬的方法研究了納米流動(dòng)過程中，孔隙尺寸對(duì)H2O分子密度分布、流動(dòng)速度及動(dòng)力粘度的影響. Markesteijn等[5]采用分子模擬研究了壓力驅(qū)動(dòng)流動(dòng)時(shí)，不同力場(chǎng)參數(shù)及溫度對(duì)H2O分子流動(dòng)速度的影響，驗(yàn)證了可以更加準(zhǔn)確描述H2O分子納米流動(dòng)的力場(chǎng)模型. Wang等[6]通過分子模擬的方法分析了超臨界辛烷在納米級(jí)巖石孔隙內(nèi)的流動(dòng)速度分布，發(fā)現(xiàn)溫度和孔隙尺寸均會(huì)對(duì)超臨界辛烷的流動(dòng)速度產(chǎn)生影響. Jiang等[7]采用分子模擬的方法研究了孔隙大小對(duì)超臨界甲烷在硅納米通道內(nèi)流動(dòng)過程的影響，分析了甲烷分子在巖石孔隙內(nèi)的密度分布和流動(dòng)狀態(tài). 然而，對(duì)巖石納米孔隙中超臨界CO2流動(dòng)狀態(tài)的研究卻罕見報(bào)道.

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

2 分子模擬方法

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

表1 Mg2SiO4和CO2力場(chǎng)

注：ε代表L-J勢(shì)能曲線中的勢(shì)能阱深度，nm；σ代表勢(shì)能為零時(shí)，兩點(diǎn)間的距離，nm.

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

本文首先采用NPT系綜進(jìn)行5 nm的平衡模擬，使系統(tǒng)的密度和結(jié)構(gòu)保持穩(wěn)定. 在平衡模擬中，溫度及壓強(qiáng)的耦合分別采用基于Berendsen耦合器的velocity rescaling方法[11]和semi-isotropic方法[12]. 系統(tǒng)平衡后，采用NVT系綜進(jìn)行30 ns的非平衡模擬. 在非平衡動(dòng)力學(xué)模擬中，沿+Z方向?qū)γ總€(gè)原子施加一定大小的的體積力，推動(dòng)流體流動(dòng)[13].

圖1 超臨界CO2在巖石孔隙內(nèi)的流動(dòng)示意圖Fig. 1 Schematic diagram of the flow of supercritical CO2

3 模擬結(jié)果與討論

3.1 模型驗(yàn)證

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

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

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

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

圖3 Mg2SiO4的應(yīng)力-應(yīng)變曲線Fig. 3 Stress-strain curve of Mg2SiO4

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

3.2 孔隙尺寸的影響

3.2.1孔隙尺寸對(duì)速度的影響

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

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

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

3.2.2孔隙尺寸對(duì)密度的影響

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

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

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

進(jìn)一步的分析發(fā)現(xiàn)：在不同大小的巖石孔隙內(nèi)均發(fā)現(xiàn)了密度震蕩，近壁面處的密度震蕩現(xiàn)象并不會(huì)隨巖石孔隙尺寸的增大或減小而消失；同時(shí)，密度震蕩的距離L并不受孔隙尺寸大小的影響，在不同大小的巖石孔隙內(nèi)，密度震蕩的距離均持續(xù)到距離巖石表面L=1.5 nm處.

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

3.3 溫度和壓強(qiáng)的影響

3.3.1溫度和壓強(qiáng)對(duì)速度的影響

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

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

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

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

因此，超臨界CO2在鎂橄欖石孔隙內(nèi)流動(dòng)時(shí)，采用降低溫度、增加壓強(qiáng)的方法可以減弱鎂橄欖石表面對(duì)CO2分子的吸引力，使得CO2分子的流動(dòng)速度降低，防止超臨界CO2的泄漏.

3.3.2溫度和壓強(qiáng)對(duì)密度的影響

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

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

4 結(jié) 論

本文采用非平衡動(dòng)力學(xué)模擬的方法研究了超臨界CO2在巖石孔隙內(nèi)的流動(dòng)特性，分析了孔隙大小、溫度及壓強(qiáng)對(duì)CO2流動(dòng)特性的影響，得到以下結(jié)論.

(1) 超臨界CO2在巖石孔隙內(nèi)流動(dòng)時(shí)，當(dāng)巖石孔隙小于5.0 nm時(shí)，速度分布呈平推流的運(yùn)動(dòng)規(guī)律，只有當(dāng)巖石孔隙大于5.0 nm時(shí)，超臨界CO2分子在巖石孔隙內(nèi)的流動(dòng)才符合Poiseuille流動(dòng)的運(yùn)動(dòng)規(guī)律.

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

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