工程科学学报，第44卷，第X期 激发，不同的热浴产生不同的非傅立叶热输运特 使用密度来验证模型以及所选力场的合理性1,27-2 性2.由于随机性和局域性，朗之万热浴(Langevin) 我们计算了293K温度下的赤藓糖醇在驰豫过程 优于Nose-Hoover热浴Pa,因此本文采用Langevin 中的密度，结果如图3所示.在100ps的时间内， 热浴分别对热源与热沉进行控温，计算300K温度 随着驰豫时间的延长，赤藓糖醇的密度逐渐趋于 下的热导率时，热源温度设置为310K,热沉温度 稳定，因此模拟过程中设置的250ps的弛豫时间 设置为290K,整体温差为20K体系热流J=dg 可以使系统达到较为稳定的状态.在293K的温 应满足能量守恒定律，其中E为提取/注入的能量， 度下，赤藓糖醇的密度模拟值为1422.44kgm3,与 J;1为模拟时间，S 文献值1440kgm3[29的误差为1.22%，这说明模 1.3模拟过程 型合理且所选GROMOS54A7力场具有较高的计 图2为NEMD方法计算热导率示意图.在进 算精度 行模拟前首先对体系进行能量最小化处理，使体 1600 系达到能量最低状态.采用Velocity Verlet算法求 Density 1550 Average density 解运动方程，时间步长为0.5s,三个方向均采用周 期性边界条件.模型沿热流方向均分为20层，为 C1500 了防止在模拟过程中体系整体发生漂移，将z轴两 1450 293K1422.44kgm3 端的一层原子设为固定层，分别将紧邻固定层的 两层原子设置为热源与热沉.首先将体系整体置 于NPT系综（等温等压系综，Constant pressure and 1350 constant temperature)下，Nose-Hoover热浴和压浴 分别用于调节系统的温度和压力.在目标温度以 1300 0102030405060708090100 Time/ps 及1个大气压下驰豫250ps,达到稳定状态后撤掉 NPT系综.之后在NVE系综（微正则系综，Micro- 因3293K温度下赤藓糖醇密度随驰豫时间变化 Fig.3 Density of erythritol changes with a relaxation time at 293 K canonical ensemble)下分别对热源与热沉进行 Langevin控温，250ps后热源与热沉温度达到稳定 2.2碳纳米管长度对复合相变材料热导率的影响 状态后开始进行计算各层原子的温度，计算2s, 本文采用NEMD方法计算热导率，向热源/热 取最后1s的数据进行统计，获得体系的温度梯 沉中施加提取的能量如图4(a)所示，满足能量守 度和热流，进而计算得到热导率，每个算例均进行 恒定律，其斜率即为热流.体系的温度分布如图4(b) 了三次独立计算 所示，热源温度为310K,热沉温度为290K,其线 性拟合的斜率即为体系的温度梯度.为了探究碳 纳米管长度对赤藓糖醇/碳纳米管复合材料热导率 的影响，分别计算了碳纳米管长度为3.9、5.2、5.9、 7.4和8.6nm的复合结构的热导率，结果如图5所 示.随着碳纳米管轴向长度的增加，复合相变材料 热导率随之提升.纯赤藓糖醇的热导率NEMD计 算值为(0.67±0.04)WmlK-12,在本文模拟范围 内，：方向（轴向）热导率达到(4.37±0.27)Wm1K-1 (纯赤藓糖醇热导率的6.5倍)，平均热导率达到 (1.66土0.10)Wm-K-1(纯赤藓糖醇热导率的2.5 国2NEMD方法计算热导率示意图 倍)，碳纳米管的加入显著提升了热导率.碳纳米 Fig.2 Schematic model of the NEMD method for thermal conductivity 管的声子平均自由程(Phonon mean free path,MFP) calculation 与其直径成正比30，室温下约为500~700nmB1-刘， 2 结果与讨论 由于本文模拟的碳纳米管尺寸(L)显著小于其声 子平均自由程，此时声子以弹道输运为主，而尺寸 2.1模型密度验证 对这种输运模式具有显著影响)，在弹道输运模 密度是相变材料的重要热物理参数，许多研究 式下，热导率几乎随几何尺寸线性增加，当几何尺J = dE dt 激发，不同的热浴产生不同的非傅立叶热输运特 性[25] . 由于随机性和局域性，朗之万热浴（Langevin） 优于 Nosé–Hoover 热浴[26] ，因此本文采用 Langevin 热浴分别对热源与热沉进行控温，计算 300 K 温度 下的热导率时，热源温度设置为 310 K，热沉温度 设置为 290 K，整体温差为 20 K. 体系热流 ， 应满足能量守恒定律，其中 E 为提取/注入的能量， J；t 为模拟时间，s. 1.3    模拟过程 图 2 为 NEMD 方法计算热导率示意图. 在进 行模拟前首先对体系进行能量最小化处理，使体 系达到能量最低状态. 采用 Velocity Verlet 算法求 解运动方程，时间步长为 0.5 fs，三个方向均采用周 期性边界条件. 模型沿热流方向均分为 20 层，为 了防止在模拟过程中体系整体发生漂移，将 z 轴两 端的一层原子设为固定层，分别将紧邻固定层的 两层原子设置为热源与热沉. 首先将体系整体置 于 NPT 系综（等温等压系综，Constant pressure and constant temperature）下，Nosé–Hoover 热浴和压浴 分别用于调节系统的温度和压力. 在目标温度以 及 1 个大气压下驰豫 250 ps，达到稳定状态后撤掉 NPT 系综. 之后在 NVE 系综（微正则系综，Micro– canonical  ensemble） 下 分 别 对 热 源 与 热 沉 进 行 Langevin 控温，250 ps 后热源与热沉温度达到稳定 状态后开始进行计算各层原子的温度，计算 2 ns， 取最后 1 ns 的数据进行统计，获得体系的温度梯 度和热流，进而计算得到热导率，每个算例均进行 了三次独立计算. x zFixedlayer Heat source Fixedlayer Heat flux Heat sink 图 2    NEMD 方法计算热导率示意图 Fig.2    Schematic model of the NEMD method for thermal conductivity calculation 2    结果与讨论 2.1    模型密度验证 密度是相变材料的重要热物理参数，许多研究 使用密度来验证模型以及所选力场的合理性[11,27–28] . 我们计算了 293 K 温度下的赤藓糖醇在驰豫过程 中的密度，结果如图 3 所示. 在 100 ps 的时间内， 随着驰豫时间的延长，赤藓糖醇的密度逐渐趋于 稳定，因此模拟过程中设置的 250 ps 的弛豫时间 可以使系统达到较为稳定的状态. 在 293 K 的温 度下，赤藓糖醇的密度模拟值为 1422.44 kg·m–3，与 文献值 1440 kg·m–3[29] 的误差为 1.22%，这说明模 型合理且所选 GROMOS 54 A7 力场具有较高的计 算精度. 1600 Density 1550 Average density 1500 1450 1400 1350 1300 0 10 20 30 40 50 Time/ps Density/(kg·m−3 ) 60 70 80 293 K 1422.44 kg·m−3 90 100 图 3    293 K 温度下赤藓糖醇密度随驰豫时间变化 Fig.3    Density of erythritol changes with a relaxation time at 293 K 2.2    碳纳米管长度对复合相变材料热导率的影响 本文采用 NEMD 方法计算热导率，向热源/热 沉中施加/提取的能量如图 4（a）所示，满足能量守 恒定律，其斜率即为热流. 体系的温度分布如图 4（b） 所示，热源温度为 310 K，热沉温度为 290 K，其线 性拟合的斜率即为体系的温度梯度. 为了探究碳 纳米管长度对赤藓糖醇/碳纳米管复合材料热导率 的影响，分别计算了碳纳米管长度为 3.9、5.2、5.9、 7.4 和 8.6 nm 的复合结构的热导率，结果如图 5 所 示. 随着碳纳米管轴向长度的增加，复合相变材料 热导率随之提升. 纯赤藓糖醇的热导率 NEMD 计 算值为（0.67±0.04） W·m–1·K–1 [28] ，在本文模拟范围 内，z 方向（轴向）热导率达到（4.37±0.27） W·m–1·K–1 （纯赤藓糖醇热导率的 6.5 倍），平均热导率达到 （ 1.66±0.10 ）W·m–1·K–1（纯赤藓糖醇热导率的 2.5 倍），碳纳米管的加入显著提升了热导率. 碳纳米 管的声子平均自由程（Phonon mean free path, MFP） 与其直径成正比[30] ，室温下约为 500～700 nm[31–32] ， 由于本文模拟的碳纳米管尺寸（Lz）显著小于其声 子平均自由程，此时声子以弹道输运为主，而尺寸 对这种输运模式具有显著影响[33] ，在弹道输运模 式下，热导率几乎随几何尺寸线性增加，当几何尺 · 4 · 工程科学学报，第 44 卷，第 X 期