俞宏偉，李實，李金龍，朱韶華，孫成珍,＊

1中國石油勘探開發(fā)研究院，提高石油采收率國家重點實驗室，北京 100083

2中國石油吉林油田分公司，吉林 松原 138099

3西安交通大學，動力工程多相流國家重點實驗室，西安 710049

1 引言

提高原油采收率一直是石油工程領域的一個熱點話題。一次、二次采油僅能采集油藏中30%-40%的油，還有大量的原油難以脫附下來，推動采油技術得到了不斷發(fā)展1-5。近年來，注氣驅油(包括CO2、N2以及烴氣)6-8在二次采油以及三次采油過程中發(fā)揮了重要的作用；尤其是二氧化碳驅油技術作為將提高采收率與CO2封存同時實現(xiàn)的途徑，受到了越來越多的關注9-14。氣驅油技術提高原油采收率的主要機理是驅替相氣體可以很好地溶解于原油中14,15，原油膨脹導致界面張力、原油粘度與密度都降低，促進原油在油藏中的流動，進而提高石油的采收率。

氣驅油過程涉及到多種物理化學問題，如油氣兩相物理性質的變化、界面特性的演變、孔隙內的油氣多相流動等16-21。近年來針對于氣體在油相中的溶解行為，Zhang等人22和Yang等人23對CO2在油中的分散特性進行了研究，前者從模擬角度出發(fā)對CO2溶解度與辛烷密度降低/溶脹之間的關系提供了定量的認識，后者利用實驗手段對超臨界條件下CO2-烷烴的體積膨脹規(guī)律進行了分析。也有學者通過分子動力學(MD)模擬研究超臨界狀態(tài)下的氣體對水-油體系界面和輸運特性的影響24。同時，部分學者發(fā)現(xiàn)外部環(huán)境的改變會導致油氣兩相狀態(tài)發(fā)生變化，如Kiran等人25發(fā)現(xiàn)溫度對原油的溶解度和粘度有重要影響，Cao等人16和Zolghadr等人17研究了CO2-烷烴體系特性隨溫度的變化情況。油氣混相過程及界面特性會受到各種條件的影響，導致其特性有很大的差異性，但目前很少有學者從熱力學角度以及分子的微觀角度出發(fā)來研究不同氣體造成驅油效果不同的根本原因，進而得到統(tǒng)一的認識。在工業(yè)驅油過程中，有一項很重要的參數(shù)—最小混相壓力(MMP，Minimum Miscibility Pressure)，MMP是油氣兩相界面張力為0時所對應的最小氣相壓力26,27。當壓力達到MMP時，油氣兩相混合均勻，界面消失28。MMP是油氣達到混相的界限壓力，充分認識MMP的變化規(guī)律及其微觀機制對于氣驅油過程的控制及優(yōu)化都非常關鍵。

本文針對于吉林某油田的實際油組分，通過MD模擬研究CO2、N2及其1 : 1混合氣與油相的混相過程，獲得油氣密度分布、界面張力及其MMP，揭示了不同氣體組分對油氣混相過程的影響規(guī)律，并從油氣分子間相互作用的角度闡述了其微觀機制。研究結果對深入認識氣驅油混相過程中油氣兩相界面?zhèn)髻|過程具有重要的意義，并可指導氣驅提高原油采收率技術的優(yōu)化與設計，工程價值突出。

2 MD模型

2.1 模擬系統(tǒng)

模擬系統(tǒng)中，油的組分是根據(jù)吉林某油田的實際組分?；?，通過烷烴鏈長來劃分出短鏈、中長鏈和長鏈，并選取各個區(qū)域中的某一種烷烴來代替整個區(qū)域，如用C7來代替C7-C10。本模擬中的油組分具體見表1，系統(tǒng)溫度取為地層溫度(371.15 K)。

表1 模型油組分Table 1 The major component of oil in a simulation system.

構建初始模型時，在特定大小盒子里隨機生成氣體分子和烷烴分子，模型總尺寸為8 nm × 8 nm ×19.5 nm，采用三明治模型以保證體系的對稱性12,24，沿著z方向依次為氣相、油相、氣相，如圖1所示。本模擬研究了不同氣體壓力下的油氣混相過程，先利用Materials Studio軟件的不定型建模，在定溫度371.15 K、定體積確定密度和壓力的基礎上確定分子的數(shù)目，然后再將盒子按照“氣-油-氣”疊加起來。表2展示了不同壓力下的氣體分子數(shù)目(僅針對一個氣相盒子，整個體系的氣體分子數(shù)是該分子數(shù)目的2倍)。MD模擬計算采用LAMMPS (大型原子/分子大規(guī)模并行模擬器)軟件來開展。x-y方向采用周期性邊界條件，z方向為反射邊界條件。氣體分子的數(shù)目通過371.15 K下不同壓力來確定，以此模擬采油過程中向油藏中注入不同壓力的氣體。二氧化碳分子采用全原子EPM2模型，C－O鍵鍵長為0.116 nm，鍵角為180°12；N2分子采用全原子模型，N－N鍵鍵長為0.1112 nm，鍵角為180°28；烷烴分子采用聯(lián)合原子模型(United-atom model)，把CH3和CH2等效成一個原子，利用這種模型在不影響模擬結果精度的基礎上可以大大降低計算量29。鍵角采用Harmonic來描述其相互作用，二面角則選取opls模型來表示。采用PPPM方法來考慮長程靜電力，通過帶有極性項的12-6 Lennard-Jones勢能模型充分考慮了范德華力和庫侖力，如式(1)所示。作用力的截斷半徑取為10 ? (1 ? = 0.1 nm)。用Lorentz-Berthelot混合法則30來計算體系中交叉原子間的Lennard-Jones勢能參數(shù)，如式(2)所示。各個原子的勢能參數(shù)及所帶電荷值見表3。

表2 不同壓力下的氣體分子數(shù)目Table 2 The number of gas molecules with different pressure.

表3 原子的Lennard-Jones勢能參數(shù)以及電荷值Table 1 Lennard-Jones potential parameters of molecules and atomic charges.

圖1 模擬初始構型圖Fig. 1 Initial configuration of simulation system.

式(1)中，ε、σ、qi、qj、C、χ、rcut分別代表了能量參數(shù)、長度參數(shù)、原子i和原子j的電荷、靜電常數(shù)、介電常數(shù)、范德華力及短程靜電力的截斷半徑；式(2)中，ε和σ分別代表了Lennard-Jones模型中的能量參數(shù)和長度參數(shù)。

在整個模擬過程中，為了保持體系的溫度體積保持不變，采用NVT系綜，利用Nose-Hoover熱浴維持溫度，防止不同工況下油相的密度發(fā)生改變。時間步長選用1 fs，每隔10000步輸出一次原子坐標信息，總模擬時長為20 ns。模擬系統(tǒng)在10 ns左右可達到穩(wěn)定狀態(tài)，形成穩(wěn)定的界面。本文對最后5 ns的數(shù)據(jù)進行統(tǒng)計計算與分析，得到油相與氣相的密度分布、界面位置以及界面厚度，然后計算不同工況下的界面張力，進而確定不同氣體對應的MMP。

3 結果與討論

3.1 密度分布

模擬每隔5 ns進行密度分布統(tǒng)計，發(fā)現(xiàn)體系在10 ns時開始達到穩(wěn)定，15-20 ns密度分布曲線波動不大，并與10-15 ns統(tǒng)計出來的密度分布曲線對比發(fā)現(xiàn)幾乎沒有差異。因此取15-20 ns數(shù)據(jù)，計算分析出了氣體分子和油分子沿垂直于界面方向(z方向)的密度分布。圖2a，b分析了不同壓力的純CO2驅油下油氣兩相密度(ρ)在z方向(Z)上的變化情況。不難看出對于CO2氣體驅替的情況，有一部分CO2分子會在油氣界面處形成一個吸附層，吸附層處CO2密度達到最高。隨著氣體壓力的升高，CO2密度增大，而油的主相密度隨之減小，說明在驅油過程中油分子由于驅替氣體壓力的升高而膨脹，油氣兩相混合更均勻，油分子被更好的驅替出來。這一變化與Zhang等人22得到的CO2氣體與辛烷相互擴散規(guī)律是吻合的。同時，在較低壓力情況下系統(tǒng)達到穩(wěn)定時，僅有一小部分油分子與氣體混合，大部分油與CO2仍有非常明顯的界面存在；當壓力較高，油氣兩相混合程度已經(jīng)很高，界面相對比較模糊(圖2c-f)，說明此時的氣驅油效果非常好。

圖2 純CO2氣驅油下油與CO2的密度分布曲線與云圖Fig. 2 The density distribution of oil and gas phase in CO2-oil system.

圖3和圖4分別給出了CO2與N2以1 : 1比例混合和純N2驅油情況下油氣的密度分布曲線。可以看出，隨著壓力的升高油氣兩相的密度變化和純CO2驅油是一致的，氣體壓力升高導致氣相密度增大而油相密度減小，說明在驅油過程中油分子由于驅替氣體壓力的升高而膨脹，氣液兩相混合更均勻，油分子被更好的驅替出來。對于1 : 1混合氣，CO2的密度曲線在油氣界面處密度出現(xiàn)峰值，而N2沒有明顯的峰值出現(xiàn)，這是CO2在界面處受到強的吸引作用力造成的，吸引力的強弱將在下文的分子平均作用力勢(PMF，Potential of mean force)分析中得到具體解釋。

圖3 CO2與N2 1 : 1混合氣驅油下油與氣相的密度分布曲線Fig. 3 The density distribution of oil and gas phase in CO2/N2-oil system.

圖4 純N2氣驅油下油與氣相的密度分布曲線Fig. 4 The density distribution of oil and gas phase in N2-oil system.

3.2 界面厚度

根據(jù)上文分析得到的油氣兩相的密度分布曲線，以及吉布斯切面的定義(油的密度為其主相密度的10%-90%)可確定出油氣界面位置及界面厚度(δ)，如圖5a所示；對于不同的驅替氣體，可統(tǒng)計出界面厚度與壓力的關系，見圖5b。對于過高的氣體壓力未給出界面厚度，因為此時氣相與油相幾乎混合均勻，界面已經(jīng)變得很模糊，無法進行界定。從圖5b可以發(fā)現(xiàn)隨著壓力的升高，界面厚度隨之增大，這是因為更多的氣體分子從氣體主相擴散到油相中，與此同時油相中向氣相移動的分子也增多，油氣兩相混合程度增強，界面變得更寬。極端情況便是超過MMP時，油氣兩相完全混合均勻，界面不存在。另外可以看出，當驅替氣體處于同一壓力狀態(tài)時，隨著CO2組分的增加，界面變得更厚。這是因為同等氣體壓力下CO2的驅替效果好于N2，油氣的混合程度更明顯，界面更寬。但是由于不同氣體的MMP不同，他們達到接近混相狀態(tài)壓力范圍也不一樣；例如CO2在較小的壓力下就能混相，所以其界面厚度在較小壓力下就會變得很大。在接近混相的狀態(tài)下，不同氣體的界面厚度都差不多。

圖5 界面示意圖以及界面厚度與氣體壓力的關系Fig.5 Schematic diagram of interface and relationship between interface thickness and gas pressure.

3.3 界面張力與MMP

油氣界面的界面張力使用吉布斯公式可以計算得出，基于壓力張量的界面張力(δ)計算P公式如下24,31：

式(3)中，PN和PT分別代表法向與切向的壓力，Pii(i=x,y,z)是不同方向上的壓力張量，Lz為模擬系統(tǒng)中盒子沿著z方向上的尺寸；式(4)中，k為原子編號，ν為該原子的速度，N為體系中所有原子的數(shù)量，m為原子的質量，r為原子不同時刻的位移，f為原子受到的力，V為整個體系的體積。

通過對最后5 ns得到的Pii取值并進行平均運算，結合計算γ的公式，得到了不同壓力狀態(tài)下平衡后5 ns的平均界面張力，如圖6所示。在計算過程中沒有在界面張力接近于零的壓力區(qū)間進行取值計算，是因為這個時候油氣基本達到混相狀態(tài)，統(tǒng)計出來的參數(shù)精度不高。從圖6可以看出隨著壓力的升高，界面張力隨之降低，這是因為驅替氣體壓力升高，油相膨脹密度降低，氣相與油相的混合程度更強，故而界面張力隨之減小，這與Zolghadr等人17利用實驗測得的CO2與烷烴混相的規(guī)律是一致的。除此之外，Lara等人8與Ayirala等人26的結論也符合該規(guī)律。理論上，達到MMP時油氣界面的張力為零；在實際計算過程中，無法統(tǒng)計出完全混相時的界面張力。Orr等人32明確指出利用油氣界面消失可以精確地計算MMP。因此，本文通過已經(jīng)確定出的不同壓力時的界面張力，進行線性擬合并延伸17,29,33，得到零張力對應的壓力，即為MMP33。圖中有些點相對發(fā)散，這是由MD模擬中初始時刻隨機速度引起的誤差造成的，但這對MMP的計算結果影響不大。通過上述方法得到不同氣體對應的MMP，如圖7所示?？梢钥闯?，純CO2驅替時MMP為23.0 MPa，CO2與N2摩爾比為1 : 1時得到的MMP為50.7 MPa，純N2驅替油時獲得的MMP為119.0 MPa。吉林某油田實際測得的CO2MMP為22.3 MPa，模擬值與實際值吻合得很好，說明通過MD模擬完全可以精確地描述油氣混相過程及預測MMP。此外，從圖7中可以看出純CO2下的擬合線斜率最大，這說明利用純CO2驅油并逐漸升高氣相壓力時體系能夠很快達到油氣混相狀態(tài)，而另外兩條線明顯比純CO2驅油下的斜率小，意味著當體系中N2比例升高后氣相需要增加更大的壓力才能使油氣兩相混合均勻。該差異體現(xiàn)了不同氣體種類對油氣混合程度的影響，為了解釋這一現(xiàn)象，下文將對不同分子與油相之間的作用力進行詳細分析與解釋。

圖6 不同氣體壓力對界面張力以及MMP的影響Fig. 6 Effects of gas pressure on interfacial tension and MMPs.

圖7 不同氣體驅油下的MMPFig. 7 Effects of gas components on MMP.

3.4 油氣混相分子機制

從上文中可以得知，同一氣體在不同壓力下驅油效果不同，不同氣體驅油得到的MMP也存在一定的差異。為了探究影響氣驅油過程中油氣混相的分子機制，首先以N2驅油為例，計算分析了不同壓力下整個系統(tǒng)的總能量(E)隨著時間(T)變化情況，如圖8a所示?？梢钥闯觯S著氣體壓力的升高，系統(tǒng)總能量逐漸升高，油氣混合程度高。這是因為壓力越高，更多的氣體分子與油分子相互作用碰撞，系統(tǒng)的勢能與動能升高，說明油氣混合程度高。這是因為壓力越高，更多的氣體分子與油分子相互作用碰撞，系統(tǒng)的勢能與動能升高，油氣混合程度加劇。這是針對同一氣體在不同壓力下，從系統(tǒng)能量角度對油氣混合機制進行了分析。

圖8 N2驅油的能量變化曲線與不同氣體的PMF分布曲線Fig. 8 The time-varying energy profiles in N2 flooding and PMF profiles of different gas molecules.

為了探究CO2和N2兩種不同氣體驅油效果顯著差異的原因，本文從微觀作用力角度對這兩種氣體分子與油相之間的相互作用能進行了分析。利用PMF分布曲線來研究相互作用能的差異。PMF是利用LAMMPS中通過柔性諧波抑制偏置的Adaptive Biasing Force方法來進行采樣計算的。計算過程中固定油相，氣體分子依然按照NVT系綜進行運動；沿著模擬系統(tǒng)的z方向劃分網(wǎng)格，網(wǎng)格間距為0.02 nm，每個網(wǎng)格空間中的樣本數(shù)設置為40。圖8b所示為系統(tǒng)歸一化高度z/L在0.1-0.9范圍內的PMF分布曲線。從圖中可以清楚看出，在油相所在的位置處，CO2和N2的PMF分布曲線都表現(xiàn)出了吸引力，但CO2的凹陷程度要高于N2，表明了油相對CO2的吸引力明顯大于N2。所以在驅油過程中，同一壓力下CO2與油相的混合程度更高。要獲得同等的驅油效果，N2分子需要的壓力要遠高于CO2，N2驅油的MMP要遠大于CO2。這就從微觀作用力角度更深入地解釋了油氣的混合機制。不難看出，油氣混合程度的強弱受到了體系中總能量以及氣體分子與油相作用強度的影響，如圖9所示。當整個體系處于較高的能量并且驅替相與油相的吸引作用較強時，油氣分子向另一相擴散得更多，此時油相膨脹程度增加，整個系統(tǒng)更容易混合。

圖9 氣驅油過程中油氣混相分子機制示意圖Fig. 9 Schematic diagram of molecular mechanism of miscibility process in gas flooding(green parts represent miscible areas).

4 結論

本文采用MD模擬手段對吉林某油田氣驅油過程進行了模擬計算，研究了不同驅替氣體種類以及壓力工況下，油氣兩相的密度分布、界面特性以及MMP，并對氣驅油過程中油氣混相的分子機制進行了分析。結果表明，隨著驅替氣體壓力的升高，油氣兩相的界面厚度增加，油氣兩相的混合程度增加。值得注意的是，當氣體壓力足夠高時，兩相的混合程度變得很高，這時沒有清晰的界面存在。純CO2驅油的MMP遠遠小于純N2驅油，當這兩種氣體以摩爾比為1 : 1混合時，得到的MMP介于兩種純氣體驅油之間，說明在相同壓力下CO2組分越多油氣混合越容易。最后，本文通過對系統(tǒng)的總能量和不同氣體分子PMF的分析解釋了油氣混相的分子機制。同一氣體在壓力升高時，系統(tǒng)的總能量隨之升高，油氣分子間相互碰撞增強，油氣兩相混合程度增加；對于不同氣體分子，油相對CO2分子的吸引強而對N2的吸引弱，導致CO2對應的MMP較小?？傊?，本文從分子層面上獲得的氣驅油油氣混相過程的界面?zhèn)髻|特性及分子機制對氣驅油技術具有重要理論指導。