原油-CO2相互作用機理分子動力學模擬研究

李建山，高浩，鄢長灝，王石頭，王亮亮

（1.中國石油長慶油田分公司油氣工藝研究院，陜西西安 710018；2.中國石油長慶油田分公司低滲透油氣田勘探開發國家工程實驗室，陜西西安 710018；3.中國石油大學（華東）石油工程學院，山東青島 266580）

CO2驅提高原油采收率（CO2-EOR）是“雙碳”背景下全面推進碳捕集、利用與封存（Carbon Capture,Utilization and Storage，簡稱CCUS）過程中的重要一環[1-6]。與其他氣體介質相比，CO2在原油中的溶解度較高，在油藏條件下與原油發生劇烈的相間傳質作用，驅油效果顯著[7-8]。國內外學者針對CO2的驅油機理開展了大量的高溫高壓相態實驗、界面張力測試、細管實驗等室內物理模擬研究，認為CO2具有膨脹原油體積、降低原油黏度，CO2可以萃取輕烴、降低油氣界面張力、補充地層能量等機理[9-12]。但這些傳統的實驗研究是通過開展宏觀實驗進行數據測量或觀測去認識或反推CO2的驅油機理，無法更深層次地認識微觀分子層面上CO2與原油相互作用的內在機制。鑒于此，一些專家學者另辟蹊徑，借助計算機模擬手段研究分子層面不同組分之間的相互作用，用來彌補傳統宏觀實驗的缺點和不足。其中分子模擬技術是以實驗和基本理論為基礎建立起來的一套研究分子體系物理、化學性質模型和算法的方法[13]，該技術可以很好地模擬研究分子的可混性、內聚能、潤濕性、力學性質、擴散及阻隔、表面及孔道吸附性等分子體系的動態行為[14-16]。

因此，為了更深層次更加直觀地認識微觀分子層面上CO2的驅油機理，國內外專家利用分子動力學模擬（Molecular Dynamics，簡稱MD）方法開展了一些研究工作。FANG等[17]基于MD方法研究了不同溫度下SiO2表面CO2對原油的剝蝕作用，選取密度剖面、彎曲分布、擴散系數和徑向分布函數表征了油滴從SiO2表面微觀分離的過程，明確了低溫和高溫下原油剝離的2 種模式。LIU 等[18]研究了CO2-烷烴體系的體積膨脹的微觀機理以及溫度、壓力和烷烴結構對CO2-烷烴體系體積膨脹的影響，結果認為壓力升高、溫度降低和直鏈結構能夠提高CO2-烷烴體系的膨脹系數，環狀結構不利于CO2-烷烴體系的膨脹。張軍等[19]構建了5種不同碳數正構烷烴組成的油滴-超臨界CO2分子模型，研究了油滴在超臨界CO2中的溶解行為，表明烷烴在碳數超臨界CO2中的溶解能力隨烷烴碳數增加而逐漸降低，碳數小于C18的烷烴更易溶解在超臨界CO2中，同時明確了色散作用在溶解過程中的主導作用。

在不同油藏中，CO2在輕質原油、中質原油和重質原油中的溶解程度差異十分明顯，二者的相互作用受原油組分影響較大。此外，油藏溫度和壓力也是影響二者相互作用的重要因素。但在以往的研究中并沒有全部考慮這些因素[20]。因此，在前人的基礎上利用MD 方法構建了輕質、中輕質、中質和重質烷烴油滴-CO2分子模型，然后利用MD參數量化表征了組分、溫度、壓力對烷烴分子和CO2分子間相互作用的影響，從微觀分子層面深化認識了CO2與原油間的作用機制，研究結果為室內物理模擬研究和CO2-EOR 技術的現場應用提供良好的理論基礎和指導意義。

1 理論部分

1.1 分子動力學理論

原油的組分十分復雜，要利用分子模擬技術準確研究CO2與原油間相互作用的變化相當困難，運算工作量極為龐大。因此，在模型構建中需要將原油模型簡化處理，根據前人的研究經驗，本次模擬中選用正辛烷（C8）、正十六烷（C16）、正二十四烷（C24）和正三十二烷（C32）4 種正構烷烴分別代表輕質原油、中輕質原油、中質原油和重質原油的擬組分，利用Materials Studio 2019 商業軟件中的通用力場模塊開展模擬計算，研究不同組分、溫度以及壓力條件下CO2-烷烴油滴體系間的相互作用。選用分子力場計算方法優化CO2和正構烷烴的力場參數，該力場適合于共價分子體系，能夠準確模擬和預測大多數常見的有機分子、無機共價鍵分子、聚合物、金屬、金屬氧化物和金屬鹵化物的結構和熱力學特性，其勢能函數可用式（1）表示：

式中：Et為總勢能，單位kJ/mol；Ev、Ecr兩項為鍵合能，依次代表共價鍵能和共價交叉項，單位kJ/mol；En-b為非鍵合能，單位kJ/mol。

Ev和En-b分別表示為：

式中：Eb、Ea、Eto、Eo和EUB分別為鍵伸縮能、鍵角彎曲能、二面角扭曲能、離平面相互作用項和Urey-Bradley項（圖1），單位kJ/mol；Evdw、Ec和Eh分別為短程作用的范德華作用能、長程作用的靜電（庫倫）相互作用能和氫鍵作用能，單位kJ/mol。

圖1 鍵合能和非鍵合能演示圖[21]Fig.1 Graphical demonstration of non-bonded and bonded potentials[21]

其中范德華作用能可表示為：

庫倫作用可表示為：

式中：i與j代表不同的原子；ε為相對介電常數；rij°為Lennard-Jones 半徑，單位m；rij為原子間的距離，單位m；q為離子所帶電荷，單位C；2(rij°/rij)9和3(rij°/rij)6分別為排斥作用能和色散作用能，單位kJ/mol。

1.2 模擬細節及步驟

利用綜合法優化方法對4 種不同碳數的烷烴分子和CO2分子模型進行能量優化。首先建立一個尺寸為6.152 nm×6.152 nm×6.152 nm 的CO2正方體盒子，初始狀態下CO2密度設置為0.784 g/cm3；然后在盒子中間創建一個直徑為4 nm 的球型烷烴油滴，取代中心區域的CO2分子，盒子中剩余的CO2分子數量為2 173 個（圖2）。為保持體系的性質恒定，創建的模型均采用三維周期性邊界條件。

圖2 油滴-CO2體系分子模型Fig.2 Molecular model of oil droplets-CO2 system

模擬計算過程中采用正則系統（NVT），該系統中認為模型中粒子數量和體積保持恒定不變，通過控制溫度變化使模型中分子接近實際狀態，進而研究體系中能量的變化。體系的溫度和壓力分別采用安德森控溫法和貝倫德森控壓法控制。CO2分子的初始速度根據麥克斯韋-玻爾茲曼分布函數隨機生成，然后運用耗散粒子動力學的優化修正韋爾萊速度算法求解經典牛頓力學方程，得到粒子在不同時刻的速度和位移。范德華作用和庫倫作用分別采用基于原子總量和基于電子群求和方法，截斷半徑均為1.1 nm。運行過程中首先固定體系中的CO2，對烷烴油滴進行500 ps結構優化，然后對CO2解除固定并對體系進行能量最小化。每組模擬運算的總時長為5 ns，步長為1 fs，設置每1 ps保存一次動力學軌跡文件，用于分析體系動力學和結構變化。

2 結果與分析

2.1 CO2與烷烴油滴間相互作用的微觀過程

2.1.1 碳鏈長度的影響

表1 為不同碳鏈長度的烷烴油滴與CO2相互作用的微觀過程，模擬計算在353.15 K、8 MPa 條件下運行，運行到5 000 ps 后，各體系構型和分子的運動軌跡變化不大，逐漸達到平衡狀態。模擬結果表明，隨著模擬運算時間增加，C8油滴的體積迅速膨脹，C8分子很快分散并與CO2分子混合，C16油滴外層的部分烷烴分子能夠與CO2分子充分接觸，而C24分子和C32油滴分子聚集在一起，位阻作用使得CO2分子很難進入油滴內部，這說明隨著烷烴增加，烷烴分子與CO2分子接觸程度減小，CO2進入油滴內部的能力減弱，油滴溶解CO2的能力越弱[22-24]。從整體來看，C8和C16油滴與超臨界CO2具有較好的混溶性，說明碳數小于C16的烷烴容易與超臨界CO2發生混相，且隨著烷烴碳鏈數減小，與CO2發生混相越容易。

表1 4種烷烴油滴與CO2相互作用微觀過程Table 1 Microscopic process of interaction between four kinds of alkane oil droplets and CO2

2.1.2 溫度的影響

圖3 為353.15、393.15、433.15 K 下C16烷烴油滴與CO2分子相互作用達到平衡狀態下的空間構型和主視圖，設置模擬運算過程中的壓力為8 MPa。可以直觀看出，353.15 K 時模型盒子中CO2分子最多，大量CO2分子與烷烴分子接觸并進入油滴內部。隨著溫度升高，CO2分子布朗運動增強、分子間距離變大，CO2分子逐漸從模型盒子中逸散出，與烷烴分子接觸的CO2分子數量急劇變少，僅有極少數CO2分子能夠進入油滴內部，C16分子在熱膨脹的作用下逐漸分散開，油滴內部空間變大、體積增大，說明油滴溶解CO2的能力隨溫度升高急劇下降，油滴體積增大主要由熱膨脹作用引起。

圖3 不同溫度下C16組成的烷烴油滴與CO2相互作用微觀構型Fig.3 Microscopic configuration of interaction of oil droplets composed of C16 and CO2 under different temperatures

2.1.3 壓力的影響

圖4 為4、8、12 MPa 下C16組成的烷烴油滴與CO2分子相互作用達到平衡狀態下的空間構型和主視圖，設置模擬運算過程中的溫度為353.15 K。模擬結果顯示，4 MPa時烷烴分子緊緊聚成團狀，僅有少量CO2分子進入至油滴內部，8 MPa 時CO2已經達到超臨界態，其擴散傳質能力大大增強，油滴外層烷烴分子在相互作用下向CO2中伸展，與CO2分子接觸面積增大，同時CO2分子擴散至油滴內部，油滴體積明顯增加。尤其當壓力升高至12 MPa，大量CO2分子的擴散溶解作用導致油滴體積顯著增加，CO2分子占據在烷烴分子間的空間中，二者混溶程度進一步增加。可以看出，CO2在油滴中的溶解度隨壓力升高而增加，尤其當CO2達到超臨界態后，CO2在油滴中的溶解度隨壓力升高顯著增加，促進了CO2和油滴混溶直至發生混相。

圖4 不同壓力下C16組成的烷烴油滴與CO2相互作用微觀構型Fig.4 Microscopic configuration of interaction of oil droplets composed of C16 and CO2 under different pressure

2.2 CO2與烷烴分子間的相互作用能

CO2-烷烴油滴空間構型的變化是二者相互作用的結果，通過計算二者之間相互作用能，可以量化表征組分、溫度和壓力對二者相互作用的影響。由于在作用過程中，CO2分子和烷烴分子間并無任何的斷鍵和成鍵作用，因此，只需要研究二者在相互作用過程中非鍵合能的變化。圖5給出了不同組分、溫度和壓力條件下CO2-烷烴體系間各種非鍵合相互作用能的大小和占比情況。

圖5 CO2-烷烴體系相互作用能Fig.5 Interaction energy of CO2-alkane system

模擬計算結果顯示，隨著烷烴碳數增加、體系溫度增加和壓力降低，CO2-烷烴體系的總勢能逐漸減小，說明二者的相互作用逐漸減弱。從能量占比來看，在3種模擬條件下靜電作用分別占體系總勢能的4.33%～7.54%、2.37%～5.58%、1.04%～12.17%，范德華排斥作用分別占體系總勢能的19.37%～29.35%、10.24%～23.87%、19.13%～28.49%，范德華色散作用分別占體系總勢能的63.11%～76.30%、70.55%～87.39%、59.34%～79.83%。色散作用能為CO2-烷烴相互作用貢獻大約60%以上的能量，因此，色散作用是主導CO2-烷烴體系相互作用的主要作用能。

2.3 CO2與烷烴油滴間的內聚能密度和溶解度參數

為了進一步厘清不同組分、溫度和壓力條件下引起CO2與烷烴之間相互作用差異的原因，引入內聚能密度和溶解度參數對二者的相互作用進行表征。內聚能密度是評價分子間作用力大小的物理量，主要表征物質在聚集時分子間相互作用的強度；溶解度參數則可以反映不同物質間的溶解特性。內聚能密度和溶解度參數之間的關系可以由式（6）計算[19]：

式中：δ為溶解度參數，單位MPa1/2；δd、δv、δh分別為溶解度參數中的色散力、偶極矩（極性力）及氫鍵作用力參數，單位MPa1/2；E為內聚能，單位kJ/mol；V為物質的摩爾體積，單位mL/mol；E/V為物質的內聚能密度（Cohesive Energy Density，簡稱CED），單位104kPa。

由于研究中的4 種烷烴分子和CO2分子皆為非極性分子，分子間不存在極性力和氫鍵作用力，因此，內聚能密度和溶解度參數主要受色散力的影響。

表2 給出了平衡狀態下烷烴油滴內聚能密度和CO2-烷烴油滴間的溶解度參數值。模擬結果顯示，烷烴碳數由C8增加至C32，油滴分子內聚能密度由1.80×104kPa增加至2.24×104kPa，溶解度參數由6.24 MPa1/2減小至2.28 MPa1/2，這主要是烷烴分子鏈長度增加使得分子間位阻作用增強，導致分子間內聚能密度增大。在相互作用過程中，CO2分子進入烷烴分子內部空間要受到烷烴分子間內聚能的約束，烷烴分子間的內聚能密度越小，CO2分子越容易進入烷烴分子所構成的油滴中，CO2分子與烷烴分子間的混溶程度越高，溶解度參數越大。對于同一種組分C16來說，當溫度由353.15 K升高至433.15K，油滴的內聚能密度由2.02×104kPa減小至0.65×104kPa，溶解度系數由4.49 MPa1/2減小至2.56 MPa1/2，這是因為溫度升高引起烷烴分子布朗運動增強，導致烷烴分子空間變大、內部結構膨脹松散，烷烴分子間色散作用弱，內聚能密度減小，同時溫度升高CO2氣體分子的布朗運動更加劇烈，氣體分子空間急劇膨脹，CO2分子傾向于從烷烴油滴中向外擴散逃逸，CO2分子與烷烴分子間的混溶程度變差，因而溶解度系數隨溫度升高而顯著減小。體系壓力由4 MPa增加至12 MPa，烷烴分子空間受到擠壓，分子間距小、色散作用強，內聚能密度由0.442×104kPa 增加至9.372×104kPa。在低壓（4 MPa）下，在油滴分子間內聚能的約束下，僅有極少量的CO2能夠溶解進油滴，溶解度參數僅為2.09 MPa1/2。當壓力超過8 MPa，此時CO2已經達到超臨界態，其擴散能力顯著增強，超臨界CO2分子能夠克服油滴內聚能增加產生的約束作用，大量CO2擴散進入至烷烴分子間，12 MPa 下溶解度參數提高至19.68 MPa1/2，溶解能力顯著增強。

表2 烷烴油滴內聚能密度和CO2-烷烴體系溶解度參數Table 2 Cohesive energy density of alkane oil droplets and solubility parameters of CO2-alkane system

2.4 CO2在烷烴分子周圍的徑向分布函數和配位數

徑向分布函數（Radial Distribution Function，簡稱RDF）可以用來描述某一基團周圍其他基團在殼層局部分布的情況，是指其他粒子出現在距離質心粒子r處空間上的概率密度同該粒子在整個空間上平均分布概率密度的比值。圖6 為平衡狀態下油滴質心處烷烴分子與周圍CO2分子間C-C 原子的徑向分布函數。隨著碳鏈長度、體系溫度和壓力的變化，徑向分布函數第一波峰（r=0.212 nm）和第一波谷（r=0.406 nm）的位置并未發生變化，但是第一波峰和第一波谷的振幅發生了明顯變化，即隨著碳鏈長度增加、溫度升高和壓力降低，第一波峰和第一波谷的振幅減小。對RDF 曲線由r=0.001 nm 至r=0.406 nm與x軸圍成的面積進行積分，可得到烷烴分子周圍CO2分子的配位數，結果見表3。顯然，CO2配位數隨碳鏈長度和溫度的增加而減小，隨壓力的升高而增大。體系中烷烴分子周圍CO2分子數量增加，表明二者之間的色散作用增強、接觸面積增加，一方面油滴溶解CO2的能力增強，另一方面大量CO2分子占據油滴內部空間，分散油滴內部烷烴分子的能力顯著增強，這有利于增加油滴體積，提高油滴彈性能量，減小油滴分子間的摩擦力，降低油滴黏度。

表3 烷烴分子周圍CO2分子的配位數Table 3 Coordination number of CO2 molecule around alkane molecule

圖6 C-C徑向分布函數Fig.6 Radial distribution function（RDF）of C-C

2.5 烷烴分子在CO2中的彎曲度

彎曲度可以用來表征和量化分子鏈伸展和彎曲的程度。從模型中提取油滴表層烷烴分子自由伸展時和彎曲時的長度，彎曲度可由式7計算：

式中：χ為烷烴的彎曲度；l為烷烴在CO2相中的分子鏈末端距，單位nm；l0為烷烴自由伸展時分子鏈末端距，單位nm。

彎曲度越接近0，這說明烷烴分子在CO2中越舒展。

圖7 為平衡狀態下油滴表層烷烴分子在CO2相中彎曲度的概率分布。結果表明：①烷烴分子彎曲度的概率峰隨鏈長增加而增加，說明短烷烴分子間的纏結程度較低，容易在CO2相中伸展和溶解，油氣界面張力較低，容易發生混相，而長鏈烷烴由于空間位阻作用使得分子鏈與CO2的接觸面積變小，烷烴分子與CO2分子間的色散作用較弱，烷烴不易向CO2中伸展，油氣界面張力較高；②隨著系統溫度增加，烷烴分子彎曲度的概率峰向右移，說明在低溫下烷烴分子更容易向CO2中伸展，393.15 K 和433.15 K 下的概率峰幾乎重疊，主要是溫度增加使得烷烴分子因熱運動而自發地產生伸展，抵消了因溫度升高導致溶解作用減小產生的烷烴分子彎曲度減小現象；③隨著系統壓力增加，烷烴分子彎曲度的概率峰向左移，說明高壓下烷烴分子鏈彎曲度降低，烷烴分子在CO2中越容易伸展，CO2與烷烴分子鏈充分接觸，有利于提高CO2對油滴外層分子的萃取能力，降低油滴和CO2的油氣界面張力，促進油氣發生混相。

圖7 平衡構型中烷烴分子彎曲度的概率分布Fig.7 Probability distribution of tortuosity for alkane molecules in equilibrium structure configuration

3 結論

1）隨著烷烴分子鏈長減小、體系溫度降低和體系壓力增加，烷烴分子與CO2分子間的相互作用增強，油滴溶解度參數和CO2配位數增加，油滴外層分子鏈的彎曲度減小，CO2向油滴內部擴散以及對油滴表層分子的吸引萃取能力提高。色散力為CO2-烷烴分子相互作用貢獻60%以上的能量，是主導二者相互作用的主要能量。

2） 應用CO2驅油時，在溫度低于353.15 K（80 ℃）、壓力高于12 MPa 的輕質油藏和中輕質油藏（原油的等效鏈長小于C16）中要充分發揮CO2所有的驅油優勢，盡可能地實現混相驅或近混相驅。在溫度較高、壓力較低的中質和重質油藏中，CO2主要以非混相驅的形式驅油，應當充分發揮CO2溶解降黏、膨脹原油體積和補充能量的優勢。