楊思玉,焦貴省,呂文峰,楊永智,錢(qián)虎軍,賈儲(chǔ)源,湯 鈞

(1. 中國(guó)石油勘探開(kāi)發(fā)研究院,北京 100083； 2. 吉林大學(xué) 理論化學(xué)研究所,理論化學(xué)計(jì)算國(guó)家重點(diǎn)實(shí)驗(yàn)室,長(zhǎng)春 130021； 3. 吉林大學(xué) 化學(xué)學(xué)院,長(zhǎng)春 130012)

超臨界CO2中正烷烴溶解行為的多尺度計(jì)算機(jī)模擬

楊思玉1,焦貴省2,呂文峰1,楊永智1,錢(qián)虎軍2,賈儲(chǔ)源3,湯 鈞3

(1. 中國(guó)石油勘探開(kāi)發(fā)研究院,北京 100083； 2. 吉林大學(xué) 理論化學(xué)研究所,理論化學(xué)計(jì)算國(guó)家重點(diǎn)實(shí)驗(yàn)室,長(zhǎng)春 130021； 3. 吉林大學(xué) 化學(xué)學(xué)院,長(zhǎng)春 130012)

利用分子動(dòng)力學(xué)和耗散粒子動(dòng)力學(xué)相結(jié)合的多尺度計(jì)算機(jī)模擬方法研究正烷烴在超臨界CO2中的溶解行為. 先在微觀尺度利用分子動(dòng)力學(xué)模擬方法計(jì)算得到超臨界CO2和正烷烴的密度及溶度等物性參數(shù),再構(gòu)造耗散粒子粗?；Ｐ?利用耗散粒子動(dòng)力學(xué)模擬C39在超臨界CO2中的溶解行為,通過(guò)直觀圖像及序參量對(duì)其溶解行為進(jìn)行表征,并計(jì)算C39在超臨界CO2中的最小混相壓.

分子動(dòng)力學(xué)； 耗散粒子動(dòng)力學(xué)； 溶度參數(shù)； 序參量

隨著綠色化學(xué)溶劑的發(fā)展,超臨界流體在實(shí)驗(yàn)科學(xué)和工業(yè)生產(chǎn)等領(lǐng)域應(yīng)用越來(lái)越廣泛. 如超臨界CO2[1]， 其臨界溫度(304.25 K)和臨界壓強(qiáng)(7.3 MPa)均較低,易于制備， 且無(wú)毒無(wú)害， 因而應(yīng)用廣泛.

在石油開(kāi)采中,采用水驅(qū)油的方法從地下開(kāi)采石油,但在石油開(kāi)采后期， 水驅(qū)油的作用力不足以將油開(kāi)采出來(lái). 目前,可以采用CO2作為驅(qū)替劑代替水. CO2在油藏條件為高溫高壓的環(huán)境下為超流體,且超臨界CO2與石油中的輕質(zhì)油組分互溶,當(dāng)原油組分中的輕質(zhì)油組分比例較大時(shí),CO2的驅(qū)油效果較好. 但任何成分的原油,當(dāng)壓強(qiáng)低于某個(gè)臨界值(最小混相壓)時(shí),CO2不能與其混合,此時(shí)無(wú)法達(dá)到較好的驅(qū)油效果. 因此應(yīng)找到混合物的最小混相壓[2-4].

本文采用分子動(dòng)力學(xué)模擬及耗散粒子動(dòng)力學(xué)方法相結(jié)合的多尺度計(jì)算機(jī)模擬技術(shù),以正烷烴C39與超臨界CO2混合物為研究對(duì)象,通過(guò)計(jì)算混合物在不同壓強(qiáng)下的溶度參數(shù)及模擬混合物在不同壓強(qiáng)下的相行為,確定C39在超臨界CO2中的最小混相壓.

1 方法與模型

1.1經(jīng)典力場(chǎng)分子動(dòng)力學(xué)(MD)模擬

原子間的相互作用力可用Lennard-Jones作用勢(shì)描述為

其中ε和σ為所需力場(chǎng)參數(shù).

本文中的CO2和C39分子均采用文獻(xiàn)[5]的TraPPE_UA力場(chǎng)描述,在該力場(chǎng)中,將H與其所在的C原子視為一個(gè)單作用點(diǎn),即聯(lián)合原子模型.

1.2耗散粒子動(dòng)力學(xué)(DPD)模擬

耗散粒子動(dòng)力學(xué)(DPD)[6-8]的本質(zhì)為粗?；肿觿?dòng)力學(xué),即用粗?；∏虼碚麄€(gè)分子或分子片斷,DPD方法在聚合物自組裝行為[9-15]中應(yīng)用廣泛. 粗粒化粒子的運(yùn)動(dòng)遵循牛頓運(yùn)動(dòng)方程， 粒子間相互作用力[16-18]的表達(dá)式為

其中：FC為保守力;FD為耗散力;FR為隨機(jī)力.

根據(jù)漲落-耗散定理[19],FD和FR應(yīng)滿(mǎn)足如下關(guān)系：

其中:wD(rij)和wR(rij)為依賴(lài)于作用粒子對(duì)i和j之間距離rij的權(quán)重函數(shù);γ和σ為作用強(qiáng)度,一般σ=3.0. 在DPD模擬中,保守力強(qiáng)度αij是唯一需要輸入的值,該參數(shù)描述了在DPD粗粒化級(jí)別作用粒子對(duì)之間的相互作用強(qiáng)度. 本文中CO2與C39片斷之間的相互作用參數(shù)αij可從經(jīng)典力場(chǎng)動(dòng)力學(xué)模擬獲得的溶度參數(shù)[20]值求得. 溶度參數(shù)δ的表達(dá)式為

其中:Evalcum和Ebulk分別為體系在理想氣體狀態(tài)和相應(yīng)熱力學(xué)狀態(tài)體相的體系總能量,Evalcum可近似認(rèn)為是體系中所有單個(gè)分子處在理想氣態(tài)時(shí)能量的總和;V為摩爾體積.

在求得各組分的溶度參數(shù)后,通過(guò)下式計(jì)算DPD模擬作用參數(shù)α：

其中:χ為Flory-Huggins相互作用參數(shù);Vr為參考體積,在實(shí)際計(jì)算中一般采用CO2分子的摩爾體積.

1.3模型及模擬過(guò)程

所有分子動(dòng)力學(xué)模擬時(shí)間均為5 ns,時(shí)間步長(zhǎng)為1 fs,模擬采用DL_POLY2.20程序包進(jìn)行. DPD模擬采用如圖1所示的模型,即用1個(gè)DPD粒子代表1個(gè)CO2分子或烷烴鏈上的3個(gè)甲基或亞甲基單元,模擬時(shí)間為300萬(wàn)步,積分步長(zhǎng)為0.02,采用GPU加速HOOMD程序包[21]中的DPD模塊進(jìn)行模擬.

圖1 DPD中采用的烷烴和CO2粗?；Ｐ虵ig.1 Coarse-grained models of n-alkane and CO2 in DPD

2 結(jié)果與討論

2.1CO2與正烷烴C39在不同熱力學(xué)狀態(tài)下的密度及溶度參數(shù)

利用TraPPE_UA力場(chǎng),通過(guò)經(jīng)典分子動(dòng)力學(xué)模擬得到的CO2密度和溶度參數(shù)及實(shí)驗(yàn)值[22-23]隨壓強(qiáng)的變化關(guān)系如圖2所示. 由圖2可見(jiàn): 模擬結(jié)果和實(shí)驗(yàn)數(shù)據(jù)吻合較好,表明構(gòu)建的模型及參數(shù)選取合理； 密度和溶度參數(shù)隨壓強(qiáng)的增加而增加， 隨溫度的升高而降低.

圖2 用MD模擬得到的CO2密度(A)和溶度參數(shù)(B)以及實(shí)驗(yàn)值[22-23]隨壓強(qiáng)的變化關(guān)系Fig.2 CO2 density (A) and solubility parameter (B) in different pressures calculated by MD with experimental values[22-23]

利用TraPPE_UA力場(chǎng)描述不同烷烴鏈分子,計(jì)算C20正烷烴分子在壓強(qiáng)為1.38 MPa,不同溫度時(shí)的密度值,并與實(shí)驗(yàn)數(shù)據(jù)[24]進(jìn)行對(duì)比,結(jié)果如圖3所示. 由圖3可見(jiàn),模擬結(jié)果與實(shí)驗(yàn)結(jié)果相符.

圖3 用MD模擬得到C20正烷烴在1.38 MPa下的密度及實(shí)驗(yàn)值[24]隨溫度的變化關(guān)系Fig.3 n-C20 density calculated by MD at 1.38 MPa and different temperatures with experimental values[24]

圖4 C1～C40直鏈烷烴鏈分子在23.1 MPa,370 K時(shí)的密度與溶度參數(shù)值Fig.4 Density and solubility parameters of n-C1—n-C40 at 23.1 MPa,370 K

C1～C40直鏈烷烴分子在特定驅(qū)油條件下(23.1 MPa,370 K)的密度和溶度參數(shù)值如圖4所示. 由圖4可見(jiàn)： 隨著分子鏈長(zhǎng)n的增長(zhǎng),密度不斷增加,當(dāng)n>10后,密度幾乎不變; 溶度參數(shù)的變化趨勢(shì)與密度幾乎一致. CO2在該條件下的溶度參數(shù)值為8.713 MPa1/2,與C2,C3烷烴分子的溶度參數(shù)相當(dāng),而與長(zhǎng)鏈(n>10)烷烴的分子溶度參數(shù)(約15 MPa1/2)相差較大,根據(jù)相似相溶原理可知,短鏈烷烴分子(C1～C5)易溶于CO2,而長(zhǎng)鏈烷烴分子不易溶于CO2.

2.2計(jì)算CO2與正烷烴C39混合物的最小混相壓

利用TraPPE_UA力場(chǎng),通過(guò)經(jīng)典分子動(dòng)力學(xué)模擬得到的370 K時(shí)正烷烴C39在不同壓強(qiáng)下的溶度參數(shù)如圖5(A)所示. 由圖5(A)可見(jiàn),溶度參數(shù)隨壓強(qiáng)的增加呈線性增加. 利用式(5)和相應(yīng)狀態(tài)下CO2的溶度參數(shù)值及密度值可求得CO2與C39鏈節(jié)間的DPD相互作用參數(shù)值,結(jié)果如圖5(B)所示. 由圖5(B)可見(jiàn),CO2與C39烷烴鏈中鏈節(jié)間的相互作用參數(shù)隨壓強(qiáng)的增加而降低.

圖5 370 K時(shí),正烷烴C39在不同壓強(qiáng)下的溶度參數(shù)(A)和CO2與C39鏈節(jié)間的DPD相互作用參數(shù)(B)Fig.5 Solubility parameters of n-C39 at 370 K,different pressures (A) and DPD interaction parameters of CO2 and n-C39 at 370 K (B)

在DPD模擬中,1個(gè)CO2分子用1個(gè)DPD小球模擬,一條C39分子鏈用一條具有13個(gè)DPD小球的粗?；肿渔溎M,鏈上的每個(gè)DPD小球代表3個(gè)碳原子單元.

序參量ψ的計(jì)算公式如下：

其中δ=ρCO2-ρa(bǔ)lkane為體系中格點(diǎn)位置CO2與C39正烷烴鏈節(jié)的數(shù)密度差. 該函數(shù)可用于描述二元共混體系中的有序度,對(duì)于出現(xiàn)分相行為的混合物,該函數(shù)具有較大的值,對(duì)于各處混合均勻的相容態(tài),體系中各處的δ均接近于0,序參量[25]ψ值也接近于0.

(A)p=23.95 MPa,混相； (B) p=23.1 MPa,分相； (C) p=22.60 MPa,分相； (D) p=16.68 MPa,分相.藍(lán)色部分為CO2,紅色部分為烷烴.

圖7 370 K時(shí),V(CO2)∶V(C39)=1的混合物在不同壓強(qiáng)下的序參量 Fig.7 Order parameters of CO2 and n-C39 mixture (V(CO2)∶V(C39)=1) at 370 K and different pressures

不同壓強(qiáng)下,V(CO2)∶V(C39)=1混合物體系的DPD模擬直觀相圖如圖6所示. 由圖6可見(jiàn): 當(dāng)壓強(qiáng)為23.95 MPa時(shí),體系中CO2和C39正烷烴鏈混相均勻; 當(dāng)壓強(qiáng)低于23.1 MPa時(shí),混合物出現(xiàn)分相,表明混合物互不相容,該分相行為在低壓時(shí)較明顯； 當(dāng)壓強(qiáng)為16.68 MPa時(shí),CO2和C39正烷烴鏈分為互不相容的層狀相(圖6(D)). 不同壓強(qiáng)下,V(CO2)∶V(C39)=1的混合物體系平衡后的序參量如圖7所示. 由圖7可見(jiàn),當(dāng)體系處于相溶態(tài)(高壓p=23.95 MPa)時(shí),體系的序參量值較低; 當(dāng)體系完全分相(低壓p=16.68 MPa)時(shí),序參量值為0.55,表明體系中的有序度較高. 結(jié)合DPD模擬直觀相圖和序參量的轉(zhuǎn)變可知,該體系的最小混相壓為MMP≈23.1 MPa.

綜上所述,本文利用分子動(dòng)力學(xué)及耗散粒子動(dòng)力學(xué)相結(jié)合的多尺度計(jì)算機(jī)模擬方法,通過(guò)模擬計(jì)算得到了超臨界CO2及具有不同鏈長(zhǎng)的正烷烴的密度和溶度參數(shù)值,所得結(jié)果與實(shí)驗(yàn)值相符. 以C39正烷烴鏈為例,利用耗散粒子動(dòng)力學(xué),研究了烷烴/CO2二元共混物在不同壓強(qiáng)和不同溫度條件下的微相行為及烷烴在超臨界CO2中的溶解行為. 結(jié)合DPD模擬直觀相圖和序參量的轉(zhuǎn)變,計(jì)算得到了正烷烴C39在超臨界CO2中的最小混相壓.

[1]Girard E,Tassaing T,Camy S,et al. Enhancement of Poly(vinyl ester) Solubility in Supercritical CO2by Partial Fluorination: The Key Role of Polymer-Polymer Interactions [J]. J Am Chem Soc,2012,134(29): 11920-11923.

[2]Sarbu T,Styranec T,Beckman E J. Non-fluorous Polymers with Very High Solubility in Supercritical CO2down to Low Pressures [J]. Nature,2012,405: 165-168.

[3]Sarbu T,Styranec T,Beckman E J. Design and Synthesis of Low Cost,Sustainable CO2-Philes [J]. Ind Eng Chem Res,2000,39(12): 4678-4683.

[4]彭超,劉建儀,張廣東,等. 降低CO2驅(qū)油最小混相壓力新方法 [J]. 重慶科技學(xué)院學(xué)報(bào): 自然科學(xué)版,2012,14(1): 48-51. (PENG Chao,LIU Jianyi,ZHANG Guangdong,et al. A New Method of Reducing the Miscible-Phase Pressure of CO2Flooding [J]. Journal of Chongqing University of Science and Technology: Natural Science Edition,2012,14(1): 48-51.)

[5]ZHANG Ling,Siepmann J I. Pressure Dependence of the Vapor-Liquid-Liquid Phase Behavior in Ternary Mixtures Consisting ofn-Alkanes,n-Perfluoroalkanes,and Carbon Dioxide [J]. J Phys Chem B,2005,109(7): 2911-2919.

[6]Hoogerbrugge P J,Koelman J M V A. Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics [J]. Europhys Lett,1992,19(3): 155-160.

[7]Koelman J M V A,Hoogerbrugge P J. Dynamic Simulations of Hard-Sphere Suspensions under Steady Shear [J]. Europhys Lett,1993,21(3): 363-368.

[8]Groot R D,Warren P B. Dissipative Particle Dynamics: Bridging the Gap between Atomistic and Mesoscopic Simulation [J]. J Chem Phys,1997,107(11): 4423-4435.

[9]LI Xuejin,Pivkin I V,LIANG Haojun,et al. Shape Transformations of Membrane Vesicles from Amphiphilic Triblock Copolymers: A Dissipative Particle Dynamics Simulation Study [J]. Macromolecules,2009,42(8): 3195-3200.

[10]HE Pengtao,LI Xuejin,DENG Mingge,et al. Complex Micelles from the Self-assembly of Coil-Rod-Coil Amphiphilic Triblock Copolymers in Selective Solvents [J]. Soft Matter,2010,6: 1539-1546.

[11]LI Xuejin,LIU Yuan,WANG Lei,et al. Fusion and Fission Pathways of Vesicles from Amphiphilic Triblock Copolymers: A Dissipative Particle Dynamics Simulation Study [J]. Phys Chem Chem Phys,2009,11(20): 4051-4059.

[12]ZHANG Zunmin,GUO Hongxia. The Phase Behavior,Structure,and Dynamics of Rodlike Mesogens with Various Flexibility Using Dissipative Particle Dynamics Simulation [J]. J Chem Phys,133: 144911.

[13]HUANG Manxia,LI Ziqi,GUO Hongxia. The Effect of Janus Nanospheres on the Phase Separation of Immiscible Polymer Blends via Dissipative Particle Dynamics Simulations [J]. Soft Matter,2012,8: 6834-6845.

[14]陳弢,劉鴻,呂中元. 三嵌段共聚物自組裝結(jié)構(gòu)的分子調(diào)控 [J]. 吉林大學(xué)學(xué)報(bào): 理學(xué)版,2012,50(4): 798-804. (CHEN Tao,LIU Hong,Lü Zhongyuan. Molecular Control of Self-assembly Structure of Triblock Copolymers [J]. Journal of Jilin University: Science Edition,2012,50(4): 798-804.)

[15]唐元暉,何彥東,王曉琳. 耗散粒子動(dòng)力學(xué)及其應(yīng)用的新進(jìn)展 [J]. 高分子通報(bào),2012(1)： 8-14. (TANG Yuanhui,HE Yandong,WANG Xiaolin. The New Developments of Dissipative Particle Dynamics and Its Applications [J]. Polymer Bulletin,2012(1): 8-14.)

[16]QIAN Hujun,CHEN Lijun,Lü Zhongyuan,et al. Surface Diffusion Dynamics of a Single Polymer Chain in Dilute Solution [J]. Phys Rev Lett,2007,99(6): 068301.

[17]YOU Liyan,CHEN Lijun,QIAN Hujun,et al. Microphase Transitions of Perforated Lamellae of Cyclic Diblock Copolymers under Steady Shear [J]. Macromolecules,2007,40(14): 5222-5227.

[18]ZHU Youliang,Lü Zhongyuan. Phase Diagram of Spherical Particles Interacted with Harmonic Repulsions [J]. J Chem Phys,2011,134(4): 044903.

[20]Patel S,Lavasanifar A,Choi P. Application of Molecular Dynamics Simulation to Predict the Compatability between Water-Insoluble Drugs and Self-associating Poly(ethylene oxide)-b-Poly(ε-caprolactone) Block Copolymers [J]. Biomacromolecules,2008,9(11): 3014-3023.

[21]Anderson J A,Lorenz C D,Travesset A. General Purpose Molecular Dynamics Simulations Fully Implemented on Graphics Processing Units [J]. J Comput Phys,2008,227(10): 5342-5359.

[22]Span R,Wagner W. A New Equation of State for Carbon Dioxide Covering the Fluid Region from the Triple-Point Temperature to 1 100 K at Pressures up to 800 MPa [J]. J Phys Chem Ref Data,1996,25(6): 1509-1595.

[23]汪孟艷. 超臨界二氧化碳體系溶解度參數(shù)的分子動(dòng)力學(xué)模擬研究 [D]. 天津： 天津大學(xué),2007. (WANG Mengyan. Supercritical Carbon Dioxide by Molecular Dynamics Simulation [D]. Tianjin: Tianjin University,2007.)

[24]Rodden J B,Erkey C,Akgerman A. High-Temperature Diffusion,Viscosity,and Density Measurements inn-Eicosane [J]. J Chem Eng Data,1988,33(3): 344-347.

[25]Beardsley T M,Matsen W M. Effects of Polydispersity on the Order-Disorder Transition of Diblock Copolymer Melts [J]. Eur Phys J E,2008,27(3): 323-333.

(責(zé)任編輯： 單 凝)

AMulti-scaleComputerSimulationoftheSolvationBehaviorofn-AlkaneinSupercriticalCarbonDioxide

YANG Siyu1,JIAO Guisheng2,Lü Wenfeng1,YANG Yongzhi1,
QIAN Hujun2,JIA Chuyuan3,TANG Jun3
(1.ResearchInstituteofPetroleumExploration&Development,Beijing100083,China;
2.InstituteofTheoreticalChemistry,StateKeyLaboratoryofTheoreticalandComputationalChemistry,
JilinUniversity,Changchun130021,China; 3.CollegeofChemistry,JilinUniversity,Changchun130012,China)

The solvation behavior ofn-alkane in supercritical CO2was studied via the multi-scale computer simulations combining atomistic molecular dynamics (MD) and dissipative particle dynamics (DPD) techniques. Firstly,the physical properties such as density,solubility parameter were calculated from atomistic molecular dynamics simulations. Next,a dissipative particle dynamics model was constructed with the solubility parameters to simulate the solvation behavior ofn-C39 in supercritical CO2. The solubility ofn-C39 in supercritical CO2at different pressures was measured by visual images and order parameters. The minimum miscible pressure ofn-C39 in supercritical CO2was also estimated based on the order parameter values at different pressures.

molecular dynamics; dissipative particle dynamics; solubility parameter; order parameter

2013-11-06.

楊思玉(1972—),女,漢族,博士,高級(jí)工程師， 從事油氣田理論開(kāi)發(fā)的研究,E-mail: yangsiy@petrochina.com.cn. 通信作者： 湯 鈞(1967—),男,漢族,博士,教授,博士生導(dǎo)師,從事綠色聚合化學(xué)與功能材料的研究,E-mail： chemjtang@jlu.edu.cn.

國(guó)家自然科學(xué)基金(批準(zhǔn)號(hào)： 21074042).

O641; TE319

A

1671-5489(2014)05-1049-06