基于深層原油物性模擬的分子力場優(yōu)選及驗證

張英男, 李汝傳, 于順昌, 郭牧之, 燕友果, 郭文躍, 張 軍

(1.中國石油大學(xué)(華東)材料科學(xué)與工程學(xué)院，山東青島 266580; 2.中國石油大學(xué)(華東)石油工程學(xué)院，山東青島 266580)

深層油氣[1-5]成為油氣勘探與增產(chǎn)增儲的重要領(lǐng)域[6-8]。盆地深層高溫高壓條件普遍，形成了由氣態(tài)烴、石油蒸汽、液態(tài)油等構(gòu)成的復(fù)雜流體系統(tǒng)，其流體物性多變、相態(tài)特征復(fù)雜[9]。分子動力學(xué)模擬方法可細致刻畫研究對象的物化性質(zhì)[10-16]，如密度、黏度和界面張力等。當(dāng)前研究中考察的溫度、壓力條件太低(185～373 K[10]; 343～373 K, 18～30 MPa[11]; 300 K, 0.1 MPa[12]; 303～323 K, 0.1 MPa[13]; 283～333 K, 10 MPa[14]; 293～373 K, 14.5～21.16 MPa[15]; 295.15 K, 0.1 MPa[16])，遠低于深層的真實溫度和壓力條件，并且研究中采用各種不同的分子力場(OPLS[11,14,15]; PCFF[12]; COMPASS[13]; CVFF[16])。筆者采用分子動力學(xué)(MD)模擬與從頭算分子動力學(xué)(AIMD)模擬相結(jié)合的計算方法，考察深層原油物性特征和分子力場成鍵性質(zhì)，優(yōu)選適用于深層高溫高壓條件下的分子力場，驗證分子力場的可靠性。

1 原油組分表征和分子力場

1.1 原油組分的實驗表征

選取塔里木盆地順北地區(qū)為研究區(qū)塊，研究區(qū)位于順托果勒低隆起北部，原油密度分布在0.792～0.831 g/cm3，平均值為0.803 g/cm3，呈現(xiàn)出低密度、低黏度的特征，屬于輕質(zhì)油。所用原油樣品取自順北1-2H井，取樣井深7 469～7 569.47 m，檢測儀器為Varian 3800氣相色譜儀。原油樣品的輕烴色譜分析(圖1)結(jié)果表明，直鏈烷烴質(zhì)量分數(shù)為35.27%，支鏈烷烴質(zhì)量分數(shù)為24.96%，環(huán)烷烴質(zhì)量分數(shù)為19.04%，芳香烴質(zhì)量分數(shù)為9.63%，原油組分的定性定量結(jié)果見表1。

在檢測所得的組分中，正庚烷、正辛烷、正壬烷和正癸烷的質(zhì)量分數(shù)均超過6.4%，支鏈烷烴中2-甲基庚烷質(zhì)量分數(shù)為2.91%，環(huán)烷烴中甲基環(huán)己烷的質(zhì)量分數(shù)為4.05%，芳香烴中間二甲苯的質(zhì)量分數(shù)為2.88%，常溫常壓下是氣體物質(zhì)的正丁烷的質(zhì)量分數(shù)為0.7%，這些物質(zhì)可以作為原油樣品中的典型組分。本文中選取正辛烷、2-甲基庚烷、甲基環(huán)己烷、間二甲苯和正丁烷5種典型組分作為研究對象，考察它們的物性特征。

圖1 順北1-2H井原油樣品輕烴色譜檢測結(jié)果Fig.1 Light hydrocarbon chromatography test results of crude oil sample from well Shunbei 1-2H

1.2 分子力場

分子動力學(xué)模擬是利用經(jīng)典牛頓第二定律求解體系內(nèi)各個粒子隨時間的運動軌跡，其中粒子所受的力由體系的勢能決定，其關(guān)鍵是得到不同體系的勢能函數(shù)，而描述體系中原子間相互作用力的勢能函數(shù)，即為分子力場[17]。分子力場可以看作各類型勢能的總和，其數(shù)學(xué)表達式為

U=Ubond+Uangle+Utorsion+Uvdw+Uel+Uχ.

(1)

式中,Ubond為鍵伸縮勢能;Uangle為鍵角彎曲勢能;Utorsion為二面角扭轉(zhuǎn)勢能;Uvdw為范德華非鍵結(jié)勢能;Uel為庫倫靜電勢能;Uχ為離平面振動勢能。

對于有機分子體系，常見的分子力場有OPLS力場[18]、CVFF力場[19]、COMPASS力場[20]、PCFF力場[21]。本文中選取上述4種典型分子力場，對其準(zhǔn)確性、可靠性進行考察，優(yōu)選出適用于深層高溫高壓模擬的分子力場。

2 分子動力學(xué)模擬

2.1 模型構(gòu)建

首先，利用Materials Studio軟件分別構(gòu)建正辛烷、2-甲基庚烷、甲基環(huán)己烷、間二甲苯和正丁烷的分子結(jié)構(gòu)。然后，利用軟件中的Amorphous Cell Construction模塊分別構(gòu)建5種典型組分的模擬盒子(圖2)，每種組分的模擬盒子尺寸均為5 nm×5 nm×5 nm。對每種組分的模型分別設(shè)置OPLS、CVFF、COMPASS和PCFF力場所對應(yīng)的力場參數(shù)及電荷參數(shù)，即完成5種組分模擬體系共計20個模型的構(gòu)建。本文的可視化結(jié)果采用VMD軟件[22]實現(xiàn)，并生成所需的界面構(gòu)型圖。

表1 順北1-2H井原油樣品輕烴組分及其質(zhì)量分數(shù)Table 1 Compound and its mass fraction of crude oil sample from well Shunbei 1-2H

圖2 樣品的模擬盒子Fig.2 Simulated boxes of samples

2.2 模擬設(shè)置

所有的分子動力學(xué)模擬均使用LAMMPS軟件[23]執(zhí)行。首先采用最速下降法優(yōu)化初始構(gòu)型，對體系能量最小化。然后設(shè)置模擬時間步長為1 fs，采用 NPT系綜進行3 ns的分子動力學(xué)模擬，模擬的前2 ns用于平衡體系，后1 ns用于統(tǒng)計分析數(shù)據(jù)。模擬時的溫度和壓力根據(jù)研究區(qū)實際條件進行設(shè)置，鉆井、測井和錄井資料表明，研究區(qū)及附近區(qū)塊的地層溫度分布在373～473 K，主體溫度在433 K，地層壓力分布在70～120 MPa，主體壓力在80 MPa(圖3)。模擬盒子采用三維周期性邊界條件，截斷半徑為1.2 nm，長程靜電力使用PPPM方法計算。

圖3 研究區(qū)溫度和壓力特征Fig.3 Characteristics of formation temperature and pressure in research area

3 分子力場優(yōu)選

3.1 密度驗證

考察5種原油組分體系的密度特征，采用AP1700物質(zhì)物性平臺[24]中提供的物性數(shù)據(jù)與模擬結(jié)果進行對比，用于驗證各個分子力場的精度。AP1700包含超過2 000種物質(zhì)的物性參數(shù)，其選取超過100萬個數(shù)據(jù)與美國國家標(biāo)準(zhǔn)與技術(shù)研究院發(fā)布的數(shù)據(jù)[25]進行對比測試，其數(shù)據(jù)幾乎完全吻合。筆者對比AP1700平臺與NIST所發(fā)布的正辛烷和正丁烷的密度及黏度數(shù)據(jù)，發(fā)現(xiàn)二者在不同溫度、壓力條件的數(shù)據(jù)完全一致，但因NIST未發(fā)布2-甲基庚烷、甲基環(huán)己烷和間二甲苯的物性數(shù)據(jù)，故本文中選擇AP1700物質(zhì)物性平臺所提供的數(shù)據(jù)作為5種原油組分物性驗證的標(biāo)準(zhǔn)。模擬中采用LAMMPS軟件的空間分層方法(ave/chunk)計算密度，模擬盒子中每間隔0.5 ?為一層，每2 ps輸出一次計算結(jié)果，統(tǒng)計模擬最后1 ns的密度平均值。基于研究區(qū)的溫度和壓力條件，對原油組分的密度特征進行兩種模式的考察，分別在不同溫度和恒定壓力(80 MPa)以及不同壓力和恒定溫度(433 K)的條件下開展模擬。

恒定壓力和變化溫度條件下，分子動力學(xué)模擬后的密度結(jié)果如圖4所示。由圖4看出，不同分子力場下的密度結(jié)果均與AP1700物性數(shù)據(jù)呈現(xiàn)相同的變化趨勢，即密度隨溫度的升高而不斷降低。對于正辛烷、2-甲基庚烷、甲基環(huán)己烷和間二甲苯4種物質(zhì)，OPLS和COMPASS力場的計算結(jié)果與AP1700數(shù)據(jù)十分接近，而CVFF和PCFF力場的計算結(jié)果存在一定的偏差。對于正丁烷，4種分子力場的計算結(jié)果與AP1700數(shù)據(jù)均存在偏差，其中CVFF和PCFF力場相對較好，OPLS 和COMPASS力場略差。恒壓變溫條件下，基于模擬結(jié)果統(tǒng)計了5種組分在4種分子力場下的平均誤差，OPLS力場的平均誤差為2.23%，COMPASS力場為2.28%，CVFF力場為2.85%，PCFF力場為3.20%。

圖4 恒定壓力和變化溫度下樣品的密度特征Fig.4 Density characteristics of samples under constant pressure and varying temperature

恒定溫度和變化壓力條件下的密度模擬結(jié)果如圖5所示。可以看出，不同分子力場下的密度結(jié)果也均與AP1700物性數(shù)據(jù)呈現(xiàn)相同的變化趨勢，即密度隨壓力的升高而不斷增大。變壓條件與變溫條件下的密度模擬結(jié)果相似，OPLS和COMPASS力場對于正辛烷、2-甲基庚烷、甲基環(huán)己烷和間二甲苯4種物質(zhì)的密度模擬更為準(zhǔn)確，而CVFF和PCFF力場對于正丁烷的密度模擬相對更好。恒溫變壓條件下，根據(jù)模擬結(jié)果統(tǒng)計了5種組分在4種分子力場下的平均誤差，OPLS力場的平均誤差為2.04%，COMPASS力場為2.08%，CVFF力場為2.71%，PCFF力場為2.98%。不同溫度和壓力條件下的密度計算結(jié)果表明，OPLS力場與AP1700物性數(shù)據(jù)間的平均誤差最小，該力場能夠較為準(zhǔn)確地反映模擬對象的密度特征。但從誤差值來看，其他3種分子力場與OPLS力場間的差異很小。

圖5 恒定溫度和變化壓力下樣品的密度特征Fig.5 Density characteristics of samples under constant temperature and varying pressure

3.2 黏度驗證

圖6 Muller-Plathe反向非平衡分子動力學(xué) 方法計算黏度Fig.6 Muller-Plathe non-equilibrium molecular dynamics method for calculating viscosity

考察不同分子力場下原油組分體系的黏度特征。利用LAMMPS軟件中的Muller-Plathe反向非平衡分子動力學(xué)[26]方法進行黏度計算，該方法的原理是動量在模擬盒子的兩個不同層中的原子之間交換，從而在體系中產(chǎn)生剪切速度，根據(jù)剪切速度曲線(圖6)獲得體系中粒子的動量，基于線性擬合求解剪切速度曲線的斜率，進而計算體系的剪切黏度。本文中沿z軸方向?qū)δM體系施加動量交換，進行6 ns的分子動力學(xué)模擬，前4 ns用于體系馳豫，后2 ns用于統(tǒng)計計算體系的剪切速度和黏度。由于非平衡分子動力學(xué)模擬基于NVT系綜開展，其無法直接控制體系的壓力，為了保持穩(wěn)定的流體狀態(tài)，在構(gòu)建模型時預(yù)先將各原油組分體系的密度設(shè)置為溫度分別為373、393、413、433、453、473 K，壓力為80 MPa條件下對應(yīng)的密度，因此黏度驗證部分只考察溫度變化對原油組分黏度特征的影響，黏度驗證標(biāo)準(zhǔn)仍然采用AP1700物質(zhì)物性平臺提供的數(shù)據(jù)。

分子動力學(xué)模擬的黏度結(jié)果見圖7。相比于密度模擬，黏度模擬在非平衡態(tài)下進行，因此模擬結(jié)果存在一定的波動，但黏度總體上仍有隨溫度升高而下降的趨勢。通過對比各個分子力場下的模擬結(jié)果可以發(fā)現(xiàn)，OPLS力場的計算結(jié)果與AP1700物性數(shù)據(jù)最為接近，其他3種分子力場則存在較大的偏差。根據(jù)模擬結(jié)果統(tǒng)計了5種原油組分在4種分子力場下的平均誤差，OPLS力場為15.96%，COMPASS力場為48.53%，CVFF力場為48.61%，PCFF力場為50.15%。黏度計算結(jié)果表明，OPLS力場可以較為準(zhǔn)確地反映模擬對象的黏度特征。

通過密度和黏度兩個典型物性參數(shù)的對比可以看出，OPLS力場的模擬結(jié)果與參考數(shù)據(jù)間的平均誤差最小、吻合度最高，其不僅適用于平衡態(tài)的密度模擬，也較為適用于非平衡態(tài)的黏度模擬，并且對不同類型的原油分子均表現(xiàn)出良好的適用性。因此OPLS力場可以被確認為是最適用于深層高溫高壓條件下進行原油物性模擬的分子力場。

圖7 樣品的黏度特征Fig.7 Viscosity characteristics of samples

4 基于從頭算分子動力學(xué)的OPLS力場驗證

根據(jù)分子力場的勢能表達式，分子力場由分子間相互作用和分子內(nèi)相互作用兩大部分構(gòu)成。其中分子間相互作用取決于范德華相互作用(原子間距)和靜電相互作用(原子電荷)，這一項可以較為準(zhǔn)確地被勢函數(shù)表達；分子內(nèi)相互作用的精確描述是決定力場準(zhǔn)確性的關(guān)鍵，其包括原子間的鍵伸縮勢能(鍵長)、鍵彎曲勢能(鍵角)和二面角扭轉(zhuǎn)勢能(二面角)。雖然OPLS力場在計算典型原油組分的密度和黏度中表現(xiàn)理想，可能表明OPLS力場在模擬中所表達的分子內(nèi)作用對應(yīng)最優(yōu)的鍵長、鍵角及二面角特征，但是該力場對分子內(nèi)相互作用表達的準(zhǔn)確性仍需從更深層次進行驗證，即從更準(zhǔn)確的電子層次進行驗證。本文中采用從頭算分子動力學(xué)(AIMD)模擬對OPLS力場成鍵特征的準(zhǔn)確性進行驗證。

圖8 正辛烷分子內(nèi)碳原子間的成鍵及其特征Fig.8 Bonding and its characteristics between carbon atoms in n-octane molecule

從頭算分子動力學(xué)模擬[27]方法將密度泛函理論和分子動力學(xué)相結(jié)合，從電子系統(tǒng)動力學(xué)出發(fā)計算原子間相互作用勢和電子系統(tǒng)總能量，基于電子波函數(shù)正交化產(chǎn)生的“虛擬力”求解原子和電子的運動方程。AIMD方法無需預(yù)先附加力場參數(shù)，可以方便地進行分子動力學(xué)模擬，并獲得準(zhǔn)確的分子內(nèi)成鍵信息。利用Materials Studio軟件中的Dmol3模塊進行從頭算分子動力學(xué)模擬。由于直鏈烷烴、支鏈烷烴和環(huán)烷烴均屬于飽和烴，其OPLS力場的參數(shù)一致，因此只選擇正辛烷作為三者的代表進行研究，同時考察間二甲苯。模擬時的研究對象為20個分子，設(shè)置模擬時間步長為1 fs，為了更加全面地考察分子內(nèi)成鍵特征對溫度的敏感性，將模擬溫度設(shè)置在273～773 K。在NVT系綜下進行100 ps的AIMD模擬，統(tǒng)計平均最后10 ps內(nèi)正辛烷和間二甲苯分子內(nèi)碳原子間的成鍵參數(shù)，并與選用OPLS力場進行分子動力學(xué)模擬的統(tǒng)計結(jié)果進行對比，兩種分子的成鍵信息如圖8(a)和圖9(a)所示。正辛烷分子內(nèi)碳原子間的成鍵特征見圖8(b)～(d)，可以發(fā)現(xiàn)AIMD和MD模擬所得的鍵長、鍵角和二面角均沒有隨溫度的改變而產(chǎn)生明顯的變化，表明正辛烷分子的成鍵特征與溫度條件并不敏感，其成鍵特征在深層環(huán)境條件下較為穩(wěn)定。對比MD和AIMD的模擬結(jié)果，MD模擬下的鍵長和鍵角特征與AIMD模擬所得的結(jié)果很接近，吻合度較高，但二者的二面角特征存在些許偏差，AIMD模擬的平均二面角角度為156.15°，MD模擬的平均二面角角度為152.07°。主要原因在于AIMD模擬時的對象個數(shù)較MD模擬更少，正辛烷分子的直鏈結(jié)構(gòu)相對保持較好，而在MD模擬體系中包含了數(shù)百個正辛烷分子，模擬過程中分子間存在更多的碰撞和纏繞，導(dǎo)致二面角的平均值有所降低。

間二甲苯分子內(nèi)碳原子間的成鍵特征見圖9(b)～(d)。與正辛烷分子的成鍵特征類似，間二甲苯在AIMD和MD模擬后所得的鍵長、鍵角和二面角特征也沒有隨溫度的改變產(chǎn)生明顯的變化，并且兩種模擬所得的成鍵特征極為接近。由圖9(d)看出，間二甲苯分子的二面角特征并未像正辛烷分子那樣在兩種模擬下產(chǎn)生一定的偏差，間二甲苯的碳原子間二面角均值在173°，這是由于間二甲苯分子中的苯環(huán)結(jié)構(gòu)剛性較強，相鄰碳原子間只存在小幅振動，不會產(chǎn)生較大的扭轉(zhuǎn)。

圖9 間二甲苯分子內(nèi)碳原子間的成鍵及其特征Fig.9 Bonding and its characteristics between carbon atoms in m-xylene molecule

綜上所述，在當(dāng)前的深層油藏條件下，AIMD與MD模擬所得的成鍵特征沒有隨溫度改變而產(chǎn)生顯著變化，說明分子的成鍵性質(zhì)對溫度條件不敏感。OPLS力場設(shè)置下，兩種模擬所得的成鍵特征吻合良好，證實OPLS力場在分子動力學(xué)模擬中具有較高的可靠性。

5 結(jié) 論

(1)采用OPLS力場模擬得到的密度和黏度結(jié)果與AP1700參考數(shù)據(jù)最為接近，且該力場對不同類型的烷烴分子均表現(xiàn)出較好的適用性，OPLS力場可以被優(yōu)選為最適宜在深層高溫高壓條件下進行模擬的分子力場。

(2)原油分子的成鍵特征對溫度不敏感，并且OPLS力場下MD模擬的成鍵特征能夠較好地吻合AIMD的模擬結(jié)果，證實OPLS力場的可靠性，是最適用于深層高溫高壓條件下進行原油模擬的分子力場。