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

楊思玉,焦貴省,呂文峰,楊永智,錢虎軍,賈儲源,湯 鈞

(1. 中國石油勘探開發(fā)研究院,北京 100083； 2. 吉林大學(xué) 理論化學(xué)研究所,理論化學(xué)計算國家重點實驗室,長春 130021； 3. 吉林大學(xué) 化學(xué)學(xué)院,長春 130012)

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

楊思玉1,焦貴省2,呂文峰1,楊永智1,錢虎軍2,賈儲源3,湯 鈞3

(1. 中國石油勘探開發(fā)研究院,北京 100083； 2. 吉林大學(xué) 理論化學(xué)研究所,理論化學(xué)計算國家重點實驗室,長春 130021； 3. 吉林大學(xué) 化學(xué)學(xué)院,長春 130012)

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

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

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

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

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

1 方法與模型

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

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

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

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

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

耗散粒子動力學(xué)(DPD)[6-8]的本質(zhì)為粗?；肿觿恿W(xué),即用粗粒化小球代表整個分子或分子片斷,DPD方法在聚合物自組裝行為[9-15]中應(yīng)用廣泛. 粗?；Ｗ拥倪\動遵循牛頓運動方程， 粒子間相互作用力[16-18]的表達式為

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

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

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

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

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

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

1.3模型及模擬過程

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

圖1 DPD中采用的烷烴和CO2粗粒化模型Fig.1 Coarse-grained models of n-alkane and CO2 in DPD

2 結(jié)果與討論

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

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

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

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

圖3 用MD模擬得到C20正烷烴在1.38 MPa下的密度及實驗值[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ù)值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可見： 隨著分子鏈長n的增長,密度不斷增加,當(dāng)n>10后,密度幾乎不變; 溶度參數(shù)的變化趨勢與密度幾乎一致. CO2在該條件下的溶度參數(shù)值為8.713 MPa1/2,與C2,C3烷烴分子的溶度參數(shù)相當(dāng),而與長鏈(n>10)烷烴的分子溶度參數(shù)(約15 MPa1/2)相差較大,根據(jù)相似相溶原理可知,短鏈烷烴分子(C1～C5)易溶于CO2,而長鏈烷烴分子不易溶于CO2.

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

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

圖5 370 K時,正烷烴C39在不同壓強下的溶度參數(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個CO2分子用1個DPD小球模擬,一條C39分子鏈用一條具有13個DPD小球的粗?；肿渔溎M,鏈上的每個DPD小球代表3個碳原子單元.

序參量ψ的計算公式如下：

其中δ=ρCO2-ρalkane為體系中格點位置CO2與C39正烷烴鏈節(jié)的數(shù)密度差. 該函數(shù)可用于描述二元共混體系中的有序度,對于出現(xiàn)分相行為的混合物,該函數(shù)具有較大的值,對于各處混合均勻的相容態(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時,V(CO2)∶V(C39)=1的混合物在不同壓強下的序參量 Fig.7 Order parameters of CO2 and n-C39 mixture (V(CO2)∶V(C39)=1) at 370 K and different pressures

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

綜上所述,本文利用分子動力學(xué)及耗散粒子動力學(xué)相結(jié)合的多尺度計算機模擬方法,通過模擬計算得到了超臨界CO2及具有不同鏈長的正烷烴的密度和溶度參數(shù)值,所得結(jié)果與實驗值相符. 以C39正烷烴鏈為例,利用耗散粒子動力學(xué),研究了烷烴/CO2二元共混物在不同壓強和不同溫度條件下的微相行為及烷烴在超臨界CO2中的溶解行為. 結(jié)合DPD模擬直觀相圖和序參量的轉(zhuǎn)變,計算得到了正烷烴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é)報: 自然科學(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é)報: 理學(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]唐元暉,何彥東,王曉琳. 耗散粒子動力學(xué)及其應(yīng)用的新進展 [J]. 高分子通報,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ù)的分子動力學(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—),女,漢族,博士,高級工程師， 從事油氣田理論開發(fā)的研究,E-mail: yangsiy@petrochina.com.cn. 通信作者： 湯 鈞(1967—),男,漢族,博士,教授,博士生導(dǎo)師,從事綠色聚合化學(xué)與功能材料的研究,E-mail： chemjtang@jlu.edu.cn.

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

O641; TE319

A

1671-5489(2014)05-1049-06