正十四烷烃相变过程的分子动力学模拟

章学来，王绪哲，王继芬，徐笑锋，华维三，房满庭

（1上海海事大学蓄冷技术研究所，上海 201306；2上海第二工业大学应用化学系，上海 201209）

近年来，相变储能技术已经在我们的日常生活中应用的十分广泛，如太阳能的热利用、建筑节能[1]、冷链运输[2]及电子器件与电池的热管理中[3]等。相变材料（PCM，phase change material）的研究又是相变储能技术中的核心部分。常用的相变储能材料一般由有机物与无机物两大部分构成，无机物通常在相变的时候出现较大的过冷度与相分离现象，十分影响储能效果，而有机物在相变时表现出一定的优势，如过冷度小、热稳定性能较好等。

JEREMY等[4]提出了一个求解正十四烷相变储能时局部熔化速率与凝固速率的理论模型。MAHDAOUI等[5]通过实验和数值研究了不同浓度的十四烷和十六烷二元混合体系的结晶过程，研究发现在不同的冷却速率下，能够从DSC曲线预测不同初始浓度十四烷和十六烷二元混合体系的固相线和液相线的温度。KARIZNOVI等[6]采用压力-体积-温度（PVT）仪，测量了甲烷+正癸烷+正十四烷三元体系在1～8 MPa的环境温度下和各种压力下的汽液平衡（VLE）数据。JIANG等[7]将Cu纳米粒子分散在十四烷中制备纳米流体，研究了纳米颗粒的体积分数与温度对导热系数的影响。FU等[8]制备了正十四烷@聚苯乙烯-二氧化硅（Tet@PS-SiO2）纳米包覆相变材料（NEPCM），实验测量该纳米包覆相变材料具有较高的比热容与导热性。SCHALBART等[9]制备了十四烷混合于水中的纳米乳液，实验验证了该纳米乳液具有较好的热稳定性，能够满足热量的储存与运输过程的需求。DORAZIO等[10]通过计算动力学研究，分析了自热重整过程中十四烷转化为氢的机理。周孙希等[11]采用低共熔法将十四烷与正辛酸按照一定比例混合研制了一种二元有机储能材料，对材料进行100次循环后发现具有较好的热物性与热稳定性能，具有很好的应用价值。段远源等[12]研究了影响十四烷纳米流体导热系数的因素，发现碳纳米管能够显著的增强导热系数且稳定性较好。王瑞杰等[13]研制了基于十四烷为代表的相变材料的微胶囊的包覆技术及温控装置，大大的提高了相变材料应用过程中的稳定性。

上述研究表明，大部分学者多集中在如何提高正十四烷烃的导热性能及热稳定性上，而对于十四烷烃发生相变过程的微观结构的变化很少有研究，因此本文采用分子动力学研究方法，构建正十四烷烃体系的无定形结构对正十四烷烃的扩散性、末端距与径向分布函数进行了研究，从微观的角度出发，研究所得到的成果将为正十四烷烃在储能中的应用提供一定的参考。

1 模拟方法

本文采用符合牛顿运动方程的经典分子动力学模拟。用于原子模拟研究的COMPASS力场是一个从头算起的力场，COMPASS力场已成功地应用于有机材料中研究正十六烷体系[14]的凝固过程。COMPASS力场势的函数形式可以表示为

式 中，Ebond、Eangle、Etorsion、Eoop、Ecross、Eelec和ELJ分别代表着键的拉伸能、弯曲能、二面角扭曲能、离平面振动能、交叉能量项、库仑静电力与范德瓦尔相互作用，我们一般用Lennard-Jones势来描述范德瓦尔相互作用，其中

式中，b为键长；q为原子电荷；rij为价角；θ为扭转角；Ø为面外角；x和ε为介电常数；H、k、V和F为系综的相关参数。

基于上述模型建立了由正十四烷所构成的纯PCM体系，由50条正十四烷烃分子链构成的无定形结构，正十四烷烃分子链初始长度为16.364 Å，采用25.4 Å×25.4 Å×25.4 Å立方体盒子，初始密度1 g/cm3，如图1所示。

图1 本文所建立的模型Fig.1 Established system model

首先对体系进行200 ps能量最小化的结构优化。接着进行退火处理，在NPT系综下，体系从240 K被加热到440 K，然后又冷却至240 K，温度间隔为20 K，在一个标准大气压1.01×10-4GPa下循环5次，进行总时间为100 ps的退火处理，退火结果如图2所示。紧接着在NVE系综下进行100 ps的动力学模拟。

经过上述不断的松弛，可以认为体系得到充分稳定，可以进行分子动力学模拟，在NPT系综下，温度由270 K升高至290 K，中间温隔2 K。每个温度下进行时间步长0.2 fs，进行500000步，100 ps的模拟。其中采用Nosé-Hoover法控温、Berendsen法控压对于相变材料的模拟具有较好的效果[15]。

图2 正十四烷退火处理Fig.2 N-tetradecane annealing

2 结果与讨论

扩散系数D是表征物质传递与输送过程的重要参数，对于一般的实验方法很难准确的测量。根据爱因斯坦方程[16]，我们可以计算自扩散系数

式中，t表示所计算的时间；N代表原子的个数；ri代表第i个原子的位置；＜＞代表系综的平均。

一般来说，自扩散系数是温度的单值函数，随着温度的增加自扩散系数也随之增大。根据RAO等[17]研究所得，当扩散性随着温度发生较大变化时可以认为是体系发生了相变行为。对于模拟相对独立原子或者粒子而言只具有唯一的扩散系数，而本文所研究的正十四烷烃体系由C与H两种原子成键所构成的，应选择正十四烷烃分子链的扩散系数的变化来判断体系是否发生相变。如图3所示为正十四烷烃体系在温度由270 K逐步升温至290 K过程中自扩散系数随着温度的变化情况，根据数据显示当温度在278 K及以下时，正十四烷烃分子链的自扩散系数D皆小于3×1011m2/s，而当温度加热到280 K时自扩散系数D达到4.32×1011m2/s，温度高于280 K时自扩散系数皆大于3×1011m2/s，并且随着温度的逐步上升自扩散系数整体上也随之增加。自扩散系数在278 K加热到280 K时发生了较大变化，说明此时体系发生了相变行为，体系由固相转变为液相，此时液态的扩散性强于固态。由于本文所设置升温过程中的温差较小，自扩散系数的大小存在一定的波动，但整体趋势为自扩散系数随着温度的升高而增大。图4是体系加热过程中正十四烷烃的密度值随着温度升高变化情况，通过实验测量正十四烷在环境温度为284 K的常压下的密度与熔点情况，结果与模拟所得数据比较见表1，模拟值与实验值进行的比较，发现误差较小，说明模拟较为准确。

体系发生相变的前后，具体到分子结构上主要是分子链的拉伸与扭转、键的拉伸与扭转等。末端距指的是分子两末端之间的直线距离，是表征分子链形态和尺寸的一个重要参数[18]。径向分布函数为其他粒子相距参考粒子的密度，可以解释为系统的区域密度与平均密度的比，其表达式为[19]

式中，N为体系原子的数目；T为模拟的总步数；δr为设定的距离差；ΔN是介于r→r+δr间原子数目。

图3 正十四烷自扩散系数与温度的关系Fig.3 Relationship between self-diffusion coefficient and temperature of n-tetradecane

图4 正十四烷密度与温度的关系Fig.4 Relationship between density and temperature of n-tetradecane

表1 单质十四烷烃熔点、密度模拟值与实验值比较Table 1 Comparison of melting point,density simulation values of elemental tetradecane with experimental values

图5 体系中末端距分布与温度的关系Fig.5 Relationship between end-to-end distance and temperature in the system

图5、图6分别是体系末端距分布与温度的关系以及径向分布函数图。图5中的数据对比可以看出，随着温度的升高，正十四烷烃分子链的末端距由初始状态下的16.364 Å逐渐呈现很大的变化，分子链的末端距分布近似于高斯分布。图5(a)、5(b)和5(c)分别是体系在270 K、274 K和278 K情况下的末端距分布情况，此时体系为固态未发生相变，可以看出正十四烷烃分子链的末端距在14～15 Å之间分布较多，并且随着温度的升高峰值变得更小，而末端距在8～10 Å之间的分布相对有所增多，此部分随着温度的升高变得更加密集，相反对于正十四烷分子链末端距在16 Å以上的分布随着温度的升高变得更加稀疏，说明在相变发生之前正十四烷烃分子末端距离整体上呈现弯曲状态，随着温度的增加分子链的构象由最初的直链逐步转变成弯曲状态。图5(d)、5(e)和5(f)分别是体系在282 K、286 K和290 K下的末端距分布情况，此时体系已经发生相变呈现液态，变化最明显的是正十四烷烃分子链的末端距在16 Å以上的分布，随着温度的升高，末端距在16 Å以上的分布逐渐密集，可见分子链存在着拉伸的趋势，参照图4密度与温度的关系，可以看出在液态状况下正十四烷烃的密度变小，即比体积增大，而比体积增大的原因便为此时正十四烷烃分子链随着温度的升高由相变前的弯曲状态逐渐被拉伸。可得正十四烷烃分子链在固态时随着温度的升高逐渐由直链状态趋向于弯曲状态，而在相变过后分子链随着温度的升高又重新被拉伸趋向于直链状态。

从图6的径向分布函数图中可以发现，径向分布函数的最大峰出现在1.11 Å处，这表明原子1.11 Å处呈现其他分子或者原子的可能性是最大的，当r＜0.99 Å时的径向分布函数为0，此时表明原子与原子最近的距离最小的距离为0.99 Å。

图6 径向分布函数Fig.6 Radial distribution function

3 结 论

本文利用分子动力学模拟方法，采用COMPASS力场对正十四烷烃无定形结构进行模拟得出如下结论。

（1）采用NPT系综，模拟了十四烷体系从270 K升温到290 K时自扩散系数的变化情况，通过正十四烷烃分子的扩散系数随着温度出现转折可以判断出体系的相变温度在278.5 K。通过实验测量，相变温度的模拟值和密度的模拟值均与实验所得值误差较小，可以认为模拟相对准确。研究发现，随着正十四烷烃体系的自扩散系数整体上随着温度的升高呈现上升趋势，符合自扩散系数是温度的单值函数理论，并且发现在固态时正十四烷烃分子链的扩散系数均小于3×1011m2/s，而体系发生相变后，分子链的扩散系数大于3×1011m2/s。

（2）模拟分析了相变前后的十四烷烃分子链末端距分布情况与径向分布函数，研究发现正十四烷烃分子链在发生相变之前随着温度的升高，正十四烷烃的分子链构象逐渐由直链趋向于弯曲，而在相变发生后，液态十四烷烃随着温度的升高分子链的构象又逐渐由弯曲状态伸展开，趋向于直链，此时正十四烷烃体系的密度减小，即比体积增大的原因便是分子链随着温度的升高逐渐伸展。并且由径向分布函数图分析了原子周围其他分子或者原子的可能性。

（3）通过模拟分子的扩散性、末端距分布与径向分布函数等微观性质，可以为后续研究其他添加剂对于该分子体系微观结构的改变进行对比，为寻找更高效实用的相变储能材料提供一定的参照。