超临界CO2中正烷烃溶解行为的多尺度计算机模拟

杨思玉,焦贵省,吕文峰,杨永智,钱虎军,贾储源,汤 钧

(1. 中国石油勘探开发研究院,北京 100083； 2. 吉林大学 理论化学研究所,理论化学计算国家重点实验室,长春 130021； 3. 吉林大学 化学学院,长春 130012)

超临界CO2中正烷烃溶解行为的多尺度计算机模拟

杨思玉1,焦贵省2,吕文峰1,杨永智1,钱虎军2,贾储源3,汤 钧3

(1. 中国石油勘探开发研究院,北京 100083； 2. 吉林大学 理论化学研究所,理论化学计算国家重点实验室,长春 130021； 3. 吉林大学 化学学院,长春 130012)

利用分子动力学和耗散粒子动力学相结合的多尺度计算机模拟方法研究正烷烃在超临界CO2中的溶解行为. 先在微观尺度利用分子动力学模拟方法计算得到超临界CO2和正烷烃的密度及溶度等物性参数,再构造耗散粒子粗粒化模型,利用耗散粒子动力学模拟C39在超临界CO2中的溶解行为,通过直观图像及序参量对其溶解行为进行表征,并计算C39在超临界CO2中的最小混相压.

分子动力学； 耗散粒子动力学； 溶度参数； 序参量

随着绿色化学溶剂的发展,超临界流体在实验科学和工业生产等领域应用越来越广泛. 如超临界CO2[1]， 其临界温度(304.25 K)和临界压强(7.3 MPa)均较低,易于制备， 且无毒无害， 因而应用广泛.

在石油开采中,采用水驱油的方法从地下开采石油,但在石油开采后期， 水驱油的作用力不足以将油开采出来. 目前,可以采用CO2作为驱替剂代替水. CO2在油藏条件为高温高压的环境下为超流体,且超临界CO2与石油中的轻质油组分互溶,当原油组分中的轻质油组分比例较大时,CO2的驱油效果较好. 但任何成分的原油,当压强低于某个临界值(最小混相压)时,CO2不能与其混合,此时无法达到较好的驱油效果. 因此应找到混合物的最小混相压[2-4].

本文采用分子动力学模拟及耗散粒子动力学方法相结合的多尺度计算机模拟技术,以正烷烃C39与超临界CO2混合物为研究对象,通过计算混合物在不同压强下的溶度参数及模拟混合物在不同压强下的相行为,确定C39在超临界CO2中的最小混相压.

1 方法与模型

1.1经典力场分子动力学(MD)模拟

原子间的相互作用力可用Lennard-Jones作用势描述为

其中ε和σ为所需力场参数.

本文中的CO2和C39分子均采用文献[5]的TraPPE_UA力场描述,在该力场中,将H与其所在的C原子视为一个单作用点,即联合原子模型.

1.2耗散粒子动力学(DPD)模拟

耗散粒子动力学(DPD)[6-8]的本质为粗粒化分子动力学,即用粗粒化小球代表整个分子或分子片断,DPD方法在聚合物自组装行为[9-15]中应用广泛. 粗粒化粒子的运动遵循牛顿运动方程， 粒子间相互作用力[16-18]的表达式为

其中：FC为保守力;FD为耗散力;FR为随机力.

根据涨落-耗散定理[19],FD和FR应满足如下关系：

其中:wD(rij)和wR(rij)为依赖于作用粒子对i和j之间距离rij的权重函数;γ和σ为作用强度,一般σ=3.0. 在DPD模拟中,保守力强度αij是唯一需要输入的值,该参数描述了在DPD粗粒化级别作用粒子对之间的相互作用强度. 本文中CO2与C39片断之间的相互作用参数αij可从经典力场动力学模拟获得的溶度参数[20]值求得. 溶度参数δ的表达式为

其中:Evalcum和Ebulk分别为体系在理想气体状态和相应热力学状态体相的体系总能量,Evalcum可近似认为是体系中所有单个分子处在理想气态时能量的总和;V为摩尔体积.

在求得各组分的溶度参数后,通过下式计算DPD模拟作用参数α：

其中:χ为Flory-Huggins相互作用参数;Vr为参考体积,在实际计算中一般采用CO2分子的摩尔体积.

1.3模型及模拟过程

所有分子动力学模拟时间均为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 结果与讨论

2.1CO2与正烷烃C39在不同热力学状态下的密度及溶度参数

利用TraPPE_UA力场,通过经典分子动力学模拟得到的CO2密度和溶度参数及实验值[22-23]随压强的变化关系如图2所示. 由图2可见: 模拟结果和实验数据吻合较好,表明构建的模型及参数选取合理； 密度和溶度参数随压强的增加而增加， 随温度的升高而降低.

图2 用MD模拟得到的CO2密度(A)和溶度参数(B)以及实验值[22-23]随压强的变化关系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,不同温度时的密度值,并与实验数据[24]进行对比,结果如图3所示. 由图3可见,模拟结果与实验结果相符.

图3 用MD模拟得到C20正烷烃在1.38 MPa下的密度及实验值[24]随温度的变化关系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时的密度与溶度参数值Fig.4 Density and solubility parameters of n-C1—n-C40 at 23.1 MPa,370 K

C1～C40直链烷烃分子在特定驱油条件下(23.1 MPa,370 K)的密度和溶度参数值如图4所示. 由图4可见： 随着分子链长n的增长,密度不断增加,当n>10后,密度几乎不变; 溶度参数的变化趋势与密度几乎一致. CO2在该条件下的溶度参数值为8.713 MPa1/2,与C2,C3烷烃分子的溶度参数相当,而与长链(n>10)烷烃的分子溶度参数(约15 MPa1/2)相差较大,根据相似相溶原理可知,短链烷烃分子(C1～C5)易溶于CO2,而长链烷烃分子不易溶于CO2.

2.2计算CO2与正烷烃C39混合物的最小混相压

利用TraPPE_UA力场,通过经典分子动力学模拟得到的370 K时正烷烃C39在不同压强下的溶度参数如图5(A)所示. 由图5(A)可见,溶度参数随压强的增加呈线性增加. 利用式(5)和相应状态下CO2的溶度参数值及密度值可求得CO2与C39链节间的DPD相互作用参数值,结果如图5(B)所示. 由图5(B)可见,CO2与C39烷烃链中链节间的相互作用参数随压强的增加而降低.

图5 370 K时,正烷烃C39在不同压强下的溶度参数(A)和CO2与C39链节间的DPD相互作用参数(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小球的粗粒化分子链模拟,链上的每个DPD小球代表3个碳原子单元.

序参量ψ的计算公式如下：

其中δ=ρCO2-ρalkane为体系中格点位置CO2与C39正烷烃链节的数密度差. 该函数可用于描述二元共混体系中的有序度,对于出现分相行为的混合物,该函数具有较大的值,对于各处混合均匀的相容态,体系中各处的δ均接近于0,序参量[25]ψ值也接近于0.

(A)p=23.95 MPa,混相； (B) p=23.1 MPa,分相； (C) p=22.60 MPa,分相； (D) p=16.68 MPa,分相.蓝色部分为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可见: 当压强为23.95 MPa时,体系中CO2和C39正烷烃链混相均匀; 当压强低于23.1 MPa时,混合物出现分相,表明混合物互不相容,该分相行为在低压时较明显； 当压强为16.68 MPa时,CO2和C39正烷烃链分为互不相容的层状相(图6(D)). 不同压强下,V(CO2)∶V(C39)=1的混合物体系平衡后的序参量如图7所示. 由图7可见,当体系处于相溶态(高压p=23.95 MPa)时,体系的序参量值较低; 当体系完全分相(低压p=16.68 MPa)时,序参量值为0.55,表明体系中的有序度较高. 结合DPD模拟直观相图和序参量的转变可知,该体系的最小混相压为MMP≈23.1 MPa.

综上所述,本文利用分子动力学及耗散粒子动力学相结合的多尺度计算机模拟方法,通过模拟计算得到了超临界CO2及具有不同链长的正烷烃的密度和溶度参数值,所得结果与实验值相符. 以C39正烷烃链为例,利用耗散粒子动力学,研究了烷烃/CO2二元共混物在不同压强和不同温度条件下的微相行为及烷烃在超临界CO2中的溶解行为. 结合DPD模拟直观相图和序参量的转变,计算得到了正烷烃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驱油最小混相压力新方法 [J]. 重庆科技学院学报: 自然科学版,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]陈弢,刘鸿,吕中元. 三嵌段共聚物自组装结构的分子调控 [J]. 吉林大学学报: 理学版,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]唐元晖,何彦东,王晓琳. 耗散粒子动力学及其应用的新进展 [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]汪孟艳. 超临界二氧化碳体系溶解度参数的分子动力学模拟研究 [D]. 天津： 天津大学,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.

(责任编辑： 单 凝)

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—),女,汉族,博士,高级工程师， 从事油气田理论开发的研究,E-mail: yangsiy@petrochina.com.cn. 通信作者： 汤 钧(1967—),男,汉族,博士,教授,博士生导师,从事绿色聚合化学与功能材料的研究,E-mail： chemjtang@jlu.edu.cn.

国家自然科学基金(批准号： 21074042).

O641; TE319

A

1671-5489(2014)05-1049-06