多环芳烃在超临界环己烷中的溶解:分子动力学模拟研究

赵向勃, 付 峰, 陈 权, 高晓明, 闫 挺

(1. 延安大学 化学与化工学院, 延安 716000; 2. 厦门海关技术中心, 厦门 361026)

1 引 言

近年来，尽管太阳能、氢能、风能等新能源发展迅速，但传统化工原料仍然在人类的社会发展中发挥着重要的作用. 原油经过长期的过度开采及使用，轻质油含量在逐渐减小. 而重油和超重油储量已占据我国现有石油储量的50%以上[1]. 重油具有粘度高、杂原子含量高的特点，传统的加氢脱碳工艺在重油加工中存在诸多局限性[2-4]. 脱碳过程会产生大量残焦，造成资源浪费[5]. 而重油中的S、N等杂原子在加氢过程中会使催化剂失活. 超临界流体(SCFs)技术是一种新型的加氢脱碳工艺，可以在不使用催化剂的情况下使得重油改质，抑制焦炭的生成[6]. 超临界流体已经广泛用于重质油轻质化、减黏等方面的研究[7-11]. 超临界水(SCW)具有低的介电常数和良好的扩散性能，使得重油在SCW中有着良好的溶解性能，再结合潜在的供氢能力，重油在超临界流体中的轻质化与减黏主要集中在SCW领域中[12-14]. 然而，近年来研究发现SCW在重油热解的过程中的供氢效果并不明显，SCW供氢能力有待深入了解[15]. 同时，Zhu等[16，17]基于量子力学的理论计算否定了SCW在涉及碳氢自由基的反应中给予氢自由基的可能性，所以越来越多的人认为SCW在重油热解过程中主要作为溶剂发挥了物理作用.

然而，水和烃类的混合体系表现出复杂的相行为. 对其复杂相行为的检测和预测一直是研究的重点[18]. 目前，水和众多低分子量碳氢化合物的相行为已经有了较大的研究[19]. 然而，重烃混合物的研究较少，因为由于混合物的不透明、乳化现象均可能干扰相边界测量. 此外，烃类高温高压下的热解过程采样困难，对这类混合物的相行为测量具有挑战性. 在之前的研究中，根据Van Konynenburg和Scott分类方案[20]，Athabasca沥青+水的伪二元混合物表现出Ⅲb型相行为[21，22]，并测量了沥青在水中的溶解度. 结果表明在含水的沥青混合物中加入溶剂时，相对于水和沥青，溶剂通常具有较低的密度和较低的临界温度，相对于沥青具有较高的水溶性. 预测富烃液相中的混合相行为和水含量会发生显著变化.

随着对重油-超临界流体体系相结构、反应机理和反应动力学的研究，研究者认识到不同的相结构及其形成的扩散环境对超临界流体中重油的非催化热解起着至关重要的作用. Morimoto等[23]基于汉森溶解度参数(HSP)提出了SCW的静态介电常数大于2.2，HSP氢键组分小于10 MPa0.5是SCW与重油良好混相的必要条件. He等[24，25]利用基团贡献法，分析了273 - 644 K、26 MPa压力下原油和水混合物的相平衡，在压力为25 MPa、温度为640 - 660 K的条件下，研究了软沥青、沥青质和水的混合物的相平衡，揭示了单相和两相区域的存在. 而Guo等[26]通过混相指数解释了SCW中重油的改质，发现混相指数的降低导致软沥青产量增加的同时沥青质和残焦的收率减少. 在恒定的混相指数下，温度增加可以导致软沥青产率的增加但是残焦的产率也会有所提升，即在不结焦的高温条件下最大限度的降低混相指数可以提高软沥青的产率. 所以重油-超临界流体体系的混相可以揭示重油-超临界流体体系中的相结构和溶解行为. 分子动力学模拟(MD)可以研究超临界状态下的复杂系统[27，28]. Jian等[29]通过MD模拟研究了SCW的结构和扩散性能. 从室温到超临界温度的升高对水中氢键的影响比沿着超临界等温线密度的变化要明显得多. 计算证实，随着温度的升高和密度的降低，每个分子的平均氢键数和团簇大小都在减少. Xin等[30]通过分子模拟研究了沥青质在超临界流体中的溶剂化反应，结果表明排斥Van der Waals相互作用对沥青质与水分子相互作用的性质和程度有主要影响.

SCW是重油改质中常用的流体，但由于多环芳烃与SCW的相互作用表现为斥力，不利于多环芳烃的溶解. 而超临界环己烷与多环芳烃之间的相互作用为引力. 因此，本文利用分子动力学方法研究了不同环数的多环芳烃以及多环芳烃混合物在超临界环己烷中的溶解过程，揭示了多环芳烃在超临界环己烷中的溶解机理. 首先采用MD模拟方法研究了纯多环芳烃或多环芳烃混合物在超临界环己烷中的溶解行为，考察了多环芳烃的尺度、多环芳烃混合物的组成对溶解过程的影响. 通过对多环芳烃在超临界环己烷中的短程溶剂结构和溶剂化自由能的研究，表征了多环芳烃与超临界环己烷环境的相互作用. 然后，通过计算多环芳烃间的内聚能密度来揭示了多环芳烃本身之间的相互作用. 本文将重油在超临界流体中改质中超临界流体的筛选，操作范围的确定提供理论指导.

2 模拟与理论方法

首先MD仿真的可靠性取决于所用的力场，COMPASS力场在烃类/SCW体系中的有效性已被证实[31，32]，因此本文所用力场采用COMPASS力场应用于多环芳烃与超临界环己烷体系，所采用的计算均基于Materials Studio 2019软件的Forcite模块. 由于重油当中所含多环芳烃的结构变化较多，因此本文只考虑没有杂环掺入和烷基取代基存在的多环芳烃. 在本次研究中，选择2环的萘(NAP)为轻多环芳烃，选择6环的苯并芘(Bghip)为重多环芳烃如图1所示. 这些多环芳烃的详细信息如表1所示.

表1 多环芳烃模型性质

图1 多环芳烃分子结构Fig. 1 Molecular structures of different PAHs.

在模拟中，首先构建充满超临界环己烷分子的立方体溶剂盒来模拟超临界环己烷环境，如图2所示，环境温度为634 - 674 K，密度为0.1 g/cm3- 0.3 g/cm3，当模拟多环芳烃的溶解时，将含有纯多环芳烃或多环芳烃混合物的油滴放在溶剂盒的中心. 在研究多环芳烃的近程溶剂结构和溶剂化自由能时，将一个多环芳烃分子放入随机溶剂盒中. 溶剂盒的剩余空间装满了特定状态下的环己烷分子.

图2 超临界环己烷溶剂盒Fig. 2 Schematic diagram ofSC-cyclohexane cube box.

多环芳烃在超临界环己烷中溶解的典型动态模拟过程如下所述：首先，构建纯多环芳烃的本体结构，对溶剂盒进行结构优化并通过NPT系综的退火操作获得了纯多环芳烃的稳定本体结构. 从多环芳烃本体结构中取出一个半径为15 Å的油滴分子，将其放置在一个立方体溶剂盒的中心，立方体溶剂盒边长为80 Å. 此外，溶剂盒中充满了环己烷分子，其密度由环己烷分子的数量来调节. 环己烷密度从0.1 g/cm3到0.3 g/cm3不等. 然后，基于NVT系综进行持续时间为500 ps的分子动力学模拟，以保证动力学平衡. 在模拟过程中，温度由NHL恒温器调节. 静电法和范德华求和法分别设置为van der Waals based和atom based.

3 结果与讨论

3.1 多环芳烃在超临界环己烷中的溶解

3.1.1溶剂密度对多环芳烃溶解度的影响

在674 K的温度下，采用MD模拟了NAP和Bghip的多环芳烃液滴在0.10 -0.30 g/cm3三种超临界环己烷密度下的溶解情况. 模拟时间在500 ps时的结果如表2所示.

表2 多环芳烃在不同密度超临界环己烷中的溶解情况

温度674 K下不同密度的压力在0.005 GPa - 0.008 GPa之间，根据观察，在不同密度的环己烷中，NAP油滴分子能够迅速地均匀分散在溶剂盒中. 随着环己烷密度的增大，多环芳烃与环己烷的体系的混相指数逐渐减小. 同时，环己烷密度的变化对NAP和Bghip的溶解有重要影响. 当环己烷密度为0.10 g/cm3时，NAP分子快速溶解于环己烷相中，球形油滴形态严重偏离最初的形态. 增加环己烷密度至0.30 g/cm3时，油滴扩散的效果较好，但是相比于0.10 g/cm3的超临界环己烷体系仍然存在NAP分子的聚集. 与此同时，对于Bghip，在0.10 g/cm3时，Bghip分子迅速溶解于超临界环己烷中. 但是相比于NAP分子溶解性稍差，在0.30 g/cm3模拟结束时，溶剂盒中只发现少量解离的Bghip分子.

从图中可以看出NAP油滴可在超临界环己烷中完全分离，但仔细观察表2的图像可以发现，在溶剂箱中仍存在NAP分子的区域富集，表明NAP发生了团聚. 由于富集区域的位置在整个模拟过程中随机变化，因此在超临界环己烷中形成的NAP团聚体是不稳定的. 而当重多环芳烃溶解到环己烷相时，Bghip在溶剂盒中出现稳定的团聚，其主要以二聚体、三聚体等类似于沥青质π-π堆积的形式存在.

3.1.2温度对多环芳烃溶解度的影响

在温度为634 - 674 K、超临界环己烷密度为0.10-0.30 g/cm3时，采用MD模拟方法研究了NAP和Bghip液滴在超临界环己烷中的溶解情况，模拟时间500 ps情况下的结果如表3和表4所示.

表4 苯并芘在不同温度下的溶解行为

当环己烷密度为0.30 g/cm3时，温度升高20/40 K对NAP液滴溶解的影响有轻微的提高. NAP分子从油滴中溶解在环己烷相中的数量基本相同. 在环己烷密度为0.10 g/cm3时，在634和654 K温度下也能观察到类似的现象，但随着温度的进一步升高，情况发生了较大的变化. 在674 K温度下，大量的NAP分子会以单体的形式溶解到环己烷相中.

当环己烷密度为0.30 g/cm3时，温度升高20/40 K对Bghip液滴的溶解与NAP分子同样具有明显的影响. 大量的Bghip分子从油滴中解离并溶解在环己烷相中. 在环己烷密度为0.10 g/cm3时，在634 K和654 K温度下也能观察到类似的现象. 但随着温度的进一步升高，发生了的变化. 在674 K温度下，大量的Bghip分子以单体、二聚体或更大规模的聚集体的形式溶解到环己烷相中，此时Bghip液滴的剩余部分可以看作是一个大尺度的团簇，类似于重油热解过程中的结焦前驱体.

从表3和表4的图像可以看出，较低环己烷密度与适宜的温度有利于NAP、Bghip在超临界环己烷中的溶解.

3.1.3多环芳烃混合物在超临界环己烷中的溶解

上文对纯NAP分子和Bghip分子的溶解行为进行了研究，而多环芳烃具有不同尺度，组成差异较大的特点，采用MD模拟方法研究了多环芳烃混合物即Bghip80NAP20(下标代表油滴中苯并芘与萘的比例)混合油滴在超临界环己烷中的溶解行为，如表5所示.

表5 Bghip80NAP20油滴在不同密度超临界环己烷中的溶解行为

相比较于纯Bghip与NAP体系，NAP与Bghip混合物体系在超临界环己烷环境下的溶解有重要影响. 模拟500 ps后，Bghip80NAP20的油滴在溶剂箱中有明显的聚集成团簇状. 在不同温度下，温度越高Bghip80NAP20的含量相较于0 ps时的含量有着明显的降低，说明温度升高会促进多环芳烃混合物在环己烷中的溶解，Bghip80NAP20油滴大部分完全分离，只有小部分的聚集. Bghip初始含量有了明显的降低，促进了多环芳烃混合物在环己烷相中的溶解，油滴发生完全解离. 环己烷相中以NAP为主要成分的多环芳烃分子数量增加，从而减小了油滴的尺寸. 导致NAP和Bghip在溶剂箱中均以单体或聚集体的形式扩散. 通过对比表456的图像，轻PAHs比重PAHs在超临界环己烷中的溶解性更好.

3.1.4多环芳烃在超临界环己烷中的溶解过程

从表6中可以看出，环己烷密度的增加对多环芳烃混合物在超临界环己烷中的溶解有不利的影响. 当环己烷密度为0.10 g/cm3时，溶剂箱中仅存多环芳烃单体和团聚体，而当环己烷密度为0.20 g/cm3或0.30 g/cm3时，溶剂箱中虽然也存在多环芳烃单体和团聚体，但是相比于0.10 g/cm3来说，溶解性似乎更差. 并且在低环己烷密度下还是在高环己烷密度下进行模拟，都能观察到类似的现象，即轻多环芳烃优先溶解到环己烷相中. 以环己烷密度为0.10 g/cm3时的模拟为例. 仅经过125 ps的模拟，就有相当数量的Bghip80NAP20分子优先溶解到环己烷相中. 而0.20 g/cm3和0.30 g/cm3在125 ps明显可以看出Bghip80NAP20分子存在明显的聚集现象.

表6 674 K下Bghip80NAP20在超临界环己烷中的溶解过程

3.2 多环芳烃与超临界环己烷相互作用

3.2.1径向分配函数

多环芳烃与超临界环己烷的混溶是由溶质-溶剂与溶质-溶质相互作用决定的. 超临界环己烷中多环芳烃的近程溶剂结构和溶剂化自由能反映了体系中溶质-溶剂间的相互作用. 为了解释多环芳烃的近程溶剂结构，计算了674 K温度下、0.10 - 0.30 g/cm3环己烷密度中gBghip- H(r)和gNAP- H(r)的径向分布函数(RDF)如图3和图4所示.

图3 674 K下0.10 - 0.30 g/cm3密度范围中超临界环己烷中Bghip径向分配函数Fig. 3 gBghip-H(r)in SC-cyclohexane at temperature of 674 K and cyclohexane densities from 0.10 to 0.30 g/cm3

图4 674 K下0.10 - 0.30 g/cm3密度范围中超临界环己烷中NAP径向分配函数Fig. 4 gNAP-H(r)in SC-cyclohexane at temperature of 674 K and cyclohexane densities from 0.10 to 0.30 g/cm3

在三种环己烷密度下，0.10 g/cm3的gBghip-H(r)与其他两者分布不同. 在Bghip分子的油滴半径内，gBghip-H(r)的值保持为0. 在2.1 Å以上半径处，gBghip-H(r)的值逐渐增大，而在0.20 g/cm3和0.30 g/cm3时在10.0 Å附近开始接近1，而后持续保持. 而Fujisawa和Kajimoto等[33，34]根据光谱表征，得出如果溶质和溶剂之间存在强吸引的非键相互作用，则容易在溶质周围形成溶剂化壳. 在0.10 g/cm3时候gBghip-H(r)发现了对应溶剂壳存在的峰，则可以认为多环芳烃与超临界环己烷间具有较强的相互作用，0.10 g/cm3的环己烷可作为催化剂促进Bghip的溶解，其详细的机理值得进一步研究. 而0.20 g/cm3和0.30 g/cm3则没有产生溶剂峰，说明高密度下重多环芳烃与超临界环己烷相互作用较弱. 而在三种密度下的gNAP- H(r)的RDF分布一致，在2.3Å处开始增加，而后在10.0 Å附近开始趋向于1，继而保持. 说明轻多环芳烃与超临界环己烷的相互作用较弱，其结果与SCW有着相似的结果.

3.2.2溶剂化自由能

计算了NAP和Bghip的溶剂化自由能，如图5所示. 634 K时，NAP溶解于超临界环己烷相中，其溶剂化自由能绝对值为4.7989 kJ/mol. Bghip分子的溶剂化自由能绝对值为21.0463 kJ/mol，其绝对值明显大于NAP分子. 当环己烷密度为0.30 g/cm3时，NAP的溶剂化自由能绝对值为7.1729 kJ/mol，高浓度的环己烷密度对多环芳烃分子在超临界环己烷环境中的溶解度有消极影响. 将Bghip和NAP在三种不同环己烷密度的超临界环己烷中溶解，结果表明，环己烷密度的变化影响多环芳烃的溶解度. 随着环己烷密度的增加，溶剂化自由能的绝对值显著增加，在0.30 g/cm3的环己烷密度下，溶剂化自由能绝对值约为7.729 kJ/mol和16.8161 kJ/mol. 另外，在相同温度不同环己烷密度下，Bghip中的溶剂化自由能绝对值相对于超临界环己烷环境中的轻多环芳烃的溶剂化自由能绝对值较大. 而溶剂化自由能绝对值越小说明溶解性越好，所以在低密度下和高温下超临界环己烷对NAP和Bghip有着良好的溶解性.

(a)不同温度下的溶剂化自由能(a)Solvation free energy at different temperatures

尽管多环芳烃的尺度和环己烷的热力学状态发生了变化，但溶剂化自由能的计算值始终为负值. Bghip分子的溶剂化自由能的绝对值总是比NAP分子更大，证实了在超临界状态下重质多环芳烃比轻多环芳烃更加稳定，同样也意味着NAP在超临界环己烷中的溶解度应该比Bghip效果好.

3.2.3内聚能密度

由前文的MD模拟结果可知，多环芳烃与超临界环己烷环境下的混相指数同时由另一个因素决定，即多环芳烃分子间的相互作用. 为了了解多环芳烃与超临界环己烷环境之间的相互作用，本文计算了NAP和Bghip在超临界环己烷环境下对应温度和压力下的内聚能密度(Cohesive energy density，CED)，计算得到的CED如图6所示.

(a)不同密度下多环芳烃内聚能密度(a)CED of PAHs at different densities

超临界环己烷环境下中多环芳烃的CED越高，代表着多环芳烃的溶解越困难. 在相同环己烷密度下6环的Bghip的CED和2环NAP的CED相差不大，并且伴随温度的升高两者的CED在逐渐减小，说明在超临界环己烷环境下对NAP和Bghip都有着不错的溶解能力，并且温度的升高对于NAP和Bghip的溶解性有所提高，这与其前文的MD结果相对应. 而对于不同密度来说，0.30 g/cm3的Bghip分子和NAP分子是0.10 g/cm3的Bghip分子和NAP分子的8倍左右. 这也证实Bghip分子比NAP分子更加稳定且不容易被溶解.

随着环己烷温度的增加，Bghip在超临界环己烷环境之间的吸引作用增强. 相比之下，Bghip的内聚能密度对环己烷密度变化更为敏感. 在674 K温度下，环己烷密度从0.10 g/cm3增加到0.30 g/cm3，导致Bghip的CED大大增加了2.64E7 j/cm3. 当环己烷密度为0.30 g/cm3时，随着温度从634 K升高到674 K，Bghip的CED只降低了2.2E5 j/cm3，表明多环芳烃之间的相互作用增加. 并且随着环己烷密度的增大超临界环己烷对NAP和Bghip的溶解能力减弱，而温度的增加，会对超临界环己烷环境下NAP和Bghip的溶解有着积极的影响.

4 结 论

本文采用分子动力学模拟方法研究了多环芳烃在超临界环己烷下的溶解度行为，有助于了解重油在超临界环己烷中的溶解行为. 结果表明在不同温度与密度下超临界环己烷对于具有2个芳香环的NAP(轻多环芳烃)均具有良好的溶解性. 而温度对于含有6个芳香环的BghiP的油滴(重多环芳烃)在超临界环己烷中的溶解影响较小. 在由多环芳烃混合物组成的油滴溶解过程中，轻环芳烃优先溶解到环己烷相，重环芳烃则集中在粒径减小的油滴中. 同时，温度升高以及密度的将降低会导致多环芳烃的内聚能密度降低，增加了多环芳烃与超临界环己烷之间的相互作用. 因此，环己烷密度越小以及温度的升高可以促进多环芳烃或多环芳烃混合物在超临界环己烷中的溶解.