张智奇 钱胜 王瑞金 朱泽飞
(杭州电子科技大学机械工程学院, 杭州 310018)
纳米流体中悬浮的纳米颗粒可以增强其导热性能已经得到广泛认可, 然而纳米流体颗粒增强传热的机理目前尚不清楚. 研究表明, 纳米颗粒的聚集是纳米流体导热系数增大的重要机制, 而且纳米颗粒聚集的形态对纳米流体的导热系数有重要影响, 但是目前的导热系数模型大多是建立在Maxwell有效介质理论的“静态”和“均匀分散”假设基础上. 本文用平衡分子动力学模拟Cu-Ar纳米流体, 采用Green-Kubo公式计算导热系数, 采用Schmidt-Ott关系式计算不同聚集形态下的分形维数. 对比导热系数与分形维数可以发现: 在相同体积分数下, 较低的分形维数会有更高的导热系数, 分析了分形维数与导热系数的定量关系. 此外, 通过径向分布函数可以看出纳米颗粒紧密聚集与松散聚集的差异, 基液分子在纳米颗粒附近的纳米薄层中处于动态平衡状态. 研究结果有助于理解纳米颗粒聚集形态对导热系数的影响机理.
近20年来, 纳米流体作为传热介质能极大地提高传热性能, 引起了众多学者的兴趣. 1995年,Choi[1,2]发现在液体中加入纳米颗粒大大提高了液体的导热性能, 尤其是在低纳米颗粒体积分数的情况下(一般小于0.01%)尤为明显. 所谓纳米流体就是纳米颗粒稳定悬浮在液体中形成的胶体. 通常使用有效导热系数来评估纳米流体的导热性能, 有效导热系数模型最初是根据Maxwell有效介质理论建立的[3], 模型假设纳米颗粒在液体中是静态和均匀分散的. 此后, Hamilton 和 Crosser[4], Xuan等[5], 王补宣等[6]以及Prasher等[7]分别考虑了非球形纳米颗粒、纳米颗粒的布朗运动、纳米颗粒外的纳米薄层以及纳米颗粒的聚集, 建立了新的有效导热系数模型. 但是, 以上这些模型都是在Maxwell模型的框架基础上得到的, 并不能摆脱纳米颗粒在液体中是静态和均匀分散的假设. 目前的实验结果、理论结果和数值模拟结果之间尙有不少矛盾之处, 原因是迄今为止纳米流体强化传热的微观机理尚不清楚. Lin等[8]研究了球形颗粒和柱状颗粒纳米流体的导热, 发现柱状颗粒的纳米流体导热强于球形颗粒的纳米流体. 李屹同等[9]针对ZnO-水纳米流体, 研究了纳米颗粒周围液体分子的密度分布, 解释了纳米层提高纳米流体导热系数的机理;但是Xue等[10]的结论是纳米颗粒周围的纳米层对热传递没有影响, 如果考虑界面热阻, 纳米层甚至会降低热导率. Keblinsk等[11]对一系列实验结果的综合分析表明, 布朗运动、纳米层和近场辐射都不是提高纳米流体热导率的关键原因, 最根本的是纳米颗粒的聚集.
在实验中很容易观察到纳米颗粒的聚集[12,13],颗粒聚集对于导热系数的影响方面已经有不少成果[14-25]. Sedighi和 Mohebbi[14]通过分子动力学研究得出, 在一定的纳米颗粒浓度下, 纳米颗粒的聚集可以提高纳米流体的导热系数. Hong和Kim[15]的研究也表明纳米颗粒的凝胶状聚合可以大大增加导热系数. 最令人惊讶的是含有Fe3O4纳米颗粒的磁流体在外磁场的作用下形成直链状聚集, 该方向上的导热系数是基液的三倍[16], 比其他聚集状态下的导热系数大得多. 尽管如此, 一些研究工作表明, 由于聚集形态各异, 纳米粒子的聚集并不一定能提高纳米流体的导热系数[17]. Thaseem和Christopher[18]也表示, 在不同聚集形态下, 5% 体积分数下的纳米流体的导热系数可以增加51%或减少32%. 为了探索颗粒聚集对有效导热系数的影响规律, 研究者已经建立了一些导热系数模型, 如Xuan等[19]较早在Maxwell有效介质理论[3]基础上建立了考虑颗粒聚集的有效导热系数模型;Feng等[20]建立了一个考虑纳米颗粒聚集的新模型, 研究表明纳米粒子越小, 聚集越多, 有效热导系数越大; Wang等[21]在 Hamilton[4]模型的基础上, 用聚集体等效半径和分形维数建立了有效导热系 数 模 型; Evans等[22]在 Bruggeman模 型 和Maxwell-Garnett理论的基础上, 建立了一个考虑纳米颗粒聚集的新模型, 并用蒙特卡罗直接模拟得到了验证; Prasher等[23]将聚集的颗粒分成在骨架和在死端两类, 分别考虑对导热的影响后进行叠加, 得到了与聚集体惯性半径、分形维数和化学维数相关的导热系数模型, 最后与模拟结果比较时也只考虑了分形维数 df= 1.8 一种情况; Xiao 研究组[24,25]以分形理论为基础建立了有效导热系数的分形模型, 模型包含静态部分和动态部分, 静态部分根据Maxwell-Garnett模型得到, 动态部分考虑了旋转椭球的Brownian运动引起的对流传热, 其中旋转引起的对流传热和尺寸分布都与分形维数相关. 总之, 目前已经建立的考虑颗粒聚集的导热系数模型, 依然以Maxwell有效介质理论为基础,再考虑聚集对尺寸分布、纳米层和布朗运动等的影响, 没有涉及纳米流体中聚集体的分形维数对导热系数的影响的定量分析, 也没有验证静态部分和动态部分简单叠加的正确性.
聚集类型有两种[26], 分别为克服第一势能壁垒的紧密聚集 (reaction limited aggregation, RLA)和克服第二势能壁垒的松散聚集(diffusion limited aggregation, DLA). Gaganpreet[27]从介观角度分析了纳米流体的传热, 结果表明: 松散聚集情况下纳米颗粒聚集体分形维数较低、规模较大、导热系数较高, 且松散聚集体生长迅速, 但容易被剪切稀化. Kang 等[28]和 Duan 等[29]的研究也表明, 聚集形态对导热系数和黏度的影响很大, 更大规模的聚集会导致更大的导热系数和更大的黏度, 产生非牛流特征, 而超声处理后, 聚集结构被打破, 导热系数减小, 黏度也会减小, 非牛顿流现象就会消失. Lin等[30]研究了Al2O3纳米颗粒在湍流作用下聚集和破碎对导热的影响, 得到了努赛尔数与聚集体平均尺寸的关系表达式.
分子动力学模拟是研究纳米流体传热强化的有效方法, 但由于计算量庞大, 只适合于模拟小系统. 最近几年有使用分子动力学模拟来探究纳米粒子聚集行为对纳米流体导热系数影响的报道[14], 在分子动力学模拟中, 计算纳米流体导热系数的方法有: 平衡态分子动力学模拟和非平衡态分子动力学模拟. 平衡态分子动力学模拟基于Green-Kubo公式, 非平衡态分子动力学模拟基于导热系数的定义. 此外, 不同分子之间的势函数难以确定和评估,Cu-Ar纳米流体的势函数有经验值, Kang等[31]利用Lennard-Jones (L-J)势描述了氩原子之间的相互作用, 并采用嵌入原子势 (embedded atom method,EAM)描述了铜原子之间的相互作用. Lee等[32]对纳米流体中颗粒聚集和非聚集状态下的热导系数进行了分子动力学模拟, 结果表明纳米颗粒聚集时的热导系数比非聚集的情况下提高了35%.李凌等[33]研究了Cu-Ar纳米流体时也得到了类似结论.
综上可知, 直链状聚集时导热系数比其他聚集状态大得多, 说明聚集形态会对导热系数产生重要影响. 然而, 迄今为止还没有纳米颗粒聚集形态对导热系数的影响的定量研究, 也未见导热系数模型采用静态部分和动态部分简单叠加的正确性验证.本文针对Cu-Ar纳米流体, 用分子动力学模拟的方法, 定量分析分形维数对纳米流体导热系数的影响, 揭示添加纳米颗粒增强纳米流体传热系数的内在机制.
与文献 [14, 34]一样, 本文选用 Cu-Ar纳米流体来研究导热系数, 这是采用分子动力学手段研究不同聚集形态下纳米流体导热系数的理想选择, 因为在Cu-Ar纳米流体中仅需考虑两体势, 而不用像以水为基液时还需要考虑点电荷等因素, 大大减小了计算量. 本文模拟中采用L-J势函数, 其与液氩的实验数据非常吻合. 分子动力学是基于牛顿第二定律来计算模型中每个原子的运动状态的一种数值模拟方法, 即
在模拟中选取Ar作为基液, Cu纳米颗粒为面心立方结构, 其晶格常数为 3.165 Å[34]. 为了简便起见, 假设Cu纳米颗粒中的所有Cu原子具有相同的速度和加速度[32]. Ar-Ar之间采用L-J势函数:
对于Cu原子之间的相互作用, 可以采用更加精准的 EAM势函数[35]. 经过试算比较, 发现Cu原子之间用EAM势和用L-J势时结果相差无几, 但用L-J势时计算量明显减小, 故采用用LJ势. Ar原子与Cu原子之间的作用可以用Berthlot混合法则求得[36]:
此外, 三种原子之间的作用关系如图1所示,所有颗粒的截断半径统一设置为10 Å, 具体结果列于表1.
图1 不同原子间的相互作用势Fig. 1. Potential functions between various particles interaction.
表1 势函数参数Table 1. Parameters of potential functions.
平衡分子动力学与非平衡分子动力学都可以计算纳米流体的导热系数, 本文选择平衡分子动力学, 原因有: 1)在平衡分子动力学中, 可以用单次模拟来计算纳米流体的导热系数, 而在非平衡分子动力学中, 需要模拟多个体系求得多个温度梯度;2) 当系统处于非平衡状态时, 非平衡分子动力学的计算稳定性不如平衡分子动力学, 而且当热通量过大时可能会引起模拟体系的崩溃; 3)与非平衡分子动力学方法相比, 平衡分子动力学方法计算纳米流体所需要的模拟体系的体积更小, 可大大减少计算量[37].
颗粒聚集体通常是由一定数量的基础颗粒聚集成随机的微结构. 由于聚集过程复杂, 纳米流体中的大多数聚集体具有较广的尺寸分布[38]. 尽管它们的大小、形状、回转半径和颗粒密度各不相同,但这些具有复杂几何形状的聚集体可以用分形维数去衡量[39].颗直径为的颗粒组成的聚集体,其回转半径为, 则Schmidt-Ott公式可表示为
在模拟中先使用NVT进行弛豫, 再在微正则系综(NVE)下收集数据, 并在三个方向上均使用周期边界条件[41]. 每个时间步各原子的位置与速度采用Verlet算法, 这一算法在Ar和Cu体系中具有良好的稳定性[42]. 在分子动力学模拟中, 时间步长是影响计算量、准确性和稳定性的关键, 太小的时间步长会导致计算量巨大, 而太大的时间步长会导致系统崩溃, 对Cu-Ar纳米流体在85 K的恒温条件下进行经过一系列试验计算后, 确定时间步长设置为 10 fs.
在计算Cu-Ar纳米流体的导热系数之前, 先计算液氩的导热系数, 以验证计算方法和程序.Sarkar和 Selvam[43]研究表明, 当颗粒数目大于500时, 系统的导热系数与系统中的颗粒数无关.本文中, 4000颗Ar原子以面心立方结构(FCC)分布在 5.72 nm × 5.72 nm × 5.72 nm 的模拟盒子中. 在弛豫阶段, 采用Nose-Hoover热浴使得系统的温度保持在85 K. 在收集数据计算导热系数阶段, 撤去热浴, 以免影响热流的计算[39]. 热流计算公式如下:
由于模拟是针对的离散时间的, 故(3)式可以离散为
如图2所示, 归一化的热流自关联函数很快收敛到0, 这样就可以进行数据采样了. 图3所示为每次计算得出的导热系数, 图中每根细线代表一次计算, 粗线代表最终的平均结果. 由于氩的熔点和沸点非常接近, 分别为84和87.5 K, 为此本文模拟85 K 下的液氩的导热系数, 数值为 0.131 W/(m·K),与文献[30]中液氩在85 K下实验得到的0.129 W/(m·K)相比, 误差仅为 1.55%.
图2 归一化的热流自关联函数Fig. 2. Normalized heat current autocorrelation function.
由于Cu-Ar纳米流体没有实验得到的导热系数, 故通过与已发表的数值模拟结果对比来判断模型的有效性. 如图4所示, 在模拟盒子的中心位置,用一颗由327颗铜原子组成的粒径为1.926 nm的铜纳米颗粒代替原本位置上共计79颗的Ar原子,则纳米颗粒的体积分数为2%. 铜原子堆积方式为面心立方堆积. 图5为Cu-Ar纳米流体的模型. 采用与计算液氩导热系数相同的方式计算前文所描述的2%体积分数的Cu-Ar纳米流体, 其结果为0.165 W/(m·K), 与文献 [30]中得到的结果相吻合, 验证了计算模型和程序的有效性.
图3 液氩在 85 K 时的导热系数Fig. 3. Thermal conductivity of Ar at 85 K.
图4 铜纳米颗粒模型Fig. 4. Model of Cu nanoparticle.
图5 Cu-Ar纳米流体模型Fig. 5. Cu-Ar nanofluid model.
纳米颗粒周围的纳米薄层被认为是纳米流体导热系数的重要影响因素之一, 但仍有不少研究人员提出相反的观点[9,10,14,34]. 为了揭示纳米颗粒周围纳米薄层增强导热的内在机理, 采用分子动力学模拟来获取在能量最小化后运行10个时间步的Cu纳米颗粒周围Ar原子的分布, 如图6所示. 由于Cu原子与Ar原子之间更强的相互作用, 部分Ar原子被吸引聚集在铜纳米颗粒周围. 纳米颗粒周围0.5 nm的Ar原子分布如图7(a)所示, 在该时间点2000时间步后, 这些原子的位置如图7(b)所示. 由图可见部分Ar原子会离开原来的纳米薄层, 但是会有大约相同数量的Ar原子被吸附, 纳米薄层达到一种动态平衡. 这可能是纳米薄层增强传热的关键机理, 因为这种动态平衡下纳米颗粒与Ar原子间可以不断产生能量交换, 也因此用界面热阻的观点解释显得不合理[10].
图6 能量最小化后的纳米薄层结构Fig. 6. Nanolayer at minimizational system energy.
图7 铜纳米颗粒周围的 Ar原子初始分布 (a)和 2000 步时的分布(b)Fig. 7. Initial distribution (a) and at 2000th time-step (b) of Ar-atoms around Cu-nanoparticle.
当纳米颗粒间的斥力较小而引力较大时, 颗粒会产生聚集. 如引言中所述, 聚集方式有紧密聚集RLA和松散聚集DLA两种[26]. 当电子浓度足够高时, 静电引力克服颗粒之间的斥力时会产生紧密聚集, 而当电子浓度中等情况下, 不能完全克服颗粒之间的斥力则会产生松散聚集. 在本文的模拟体系中, 采用两颗相同大小的纳米铜颗粒进行模拟.图8(a)和图8(b)分别显示了紧密聚集和松散聚集两种聚集, 从图8可以清楚地看出: 松散结构中的两个铜纳米颗粒之间存在一层Ar原子, 而在紧密聚集中的两个纳米颗粒之间有直接接触.
径向分布函数可以解释为局部密度(local density)与总体密度 (bulk density)之比, 写作:
图8 纳米颗粒聚集 (a)紧密聚集; (b)松散聚集Fig. 8. Aggregation of nanoparticles: (a) Compact aggregation; (b) loose aggregation.
当模拟体系中存在两颗纳米颗粒时, 情况发生了变化. 图10(a)所示为两颗纳米颗粒未发生聚集时的对分布函数, 图10(b)所示为两颗纳米颗粒发生松散聚集时的对分布函数, 观察可见两者几乎相同, 而且与图10(b)中存在单颗纳米颗粒的情形十分类似. 在图10(b)中, 两颗纳米颗粒中心的距离经计算为22.358 Å, 而图中第一个峰出现的位置距纳米颗粒的中心为11.142 Å. 前者略大于后者的2倍. 这一现象表明, 松散聚集可以理解为纳米颗粒及其周围的纳米薄层作为一个整体产生聚集,其中纳米薄层占6.358 Å, 纳米颗粒的粒径为16 Å.而在如图10(c)所示的紧密聚集中, 两颗纳米颗粒质心的距离为15.237 Å, 略小于纳米颗粒的粒径,这与图8(a)和图8(b)所示的情形一致.
一些研究表明, 在分散状态下纳米颗粒的形状会对于强化传热产生影响[7,44-46]. 然而, 在聚集状态下纳米颗粒形态对传热的影响尚无定量的研究.通过分子动力学模拟计算了Cu-Ar纳米流体的导热系数, 并根据模拟得到的纳米颗粒的质心位置计算分形维数. 为了减少计算量, 本文选用包含3—6颗纳米颗粒的小系统. 在模拟中, 纳米颗粒的大小根据其体积分数计算确定, 以便在只改变体积分数的情况下, 模拟盒中存在相同数量的纳米粒子. 计算结果列于表2和表3中.
图9 对分布函数 (a)液氩体系; (b) Cu-Ar 纳米流体Fig. 9. Pair correlation function of fluid Ar (a) and nanofluid Cu/Ar (b).
图10 两颗纳米颗粒下纳米流体的对分布函数 (a)分散; (b)松聚集; (c)紧聚集Fig. 10. Pair correlation function of nanofluid with two nanoparticles: (a) Dispersed nanoparticle; (b) loose aggregation; (c) compact aggregation.
从表2中可以看出, 分形维数在1—2之间,这与大多数实验结果不符(大部分为1.4—2.7范围), 其主要原因是本文模拟体系中的纳米颗粒数少(最多6颗), 而实验测量中是有成百上千的纳米颗粒聚集在一起. 另外从表2还可以看出, 体积分数越大, 导热系数越大, 但其增量随体积分数的增大而减小, 这与已发表的研究成果相一致[47]. 另外,分形维数越小, 对应的导热系数越大, 但由于纳米颗粒大小的不同, 会存在一些例外, 如对比表2与表3中4%体积分数下的结果, 可以发现分形维数为1.60的6个颗粒的导热系数为0.278 W/(m·K),而分形维数为1.59的4个颗粒的导热系数为0.246 W/(m·K), 其原因是同体积分数下, 6 个纳米颗粒的粒径更小. 从具有相同体积分数的3个颗粒和4个颗粒的比较中也可以得出同样的结论. 从图11中可更清楚地看出, 在体积分数和纳米颗粒尺寸一定的情况下, 随着分形维数的增加, 热导系数几乎线性减小, 而且这种关系随着体系中颗粒数的增大更加明显.
表2 3—4 颗纳米颗粒下的导热系数与分形维数Table 2. Thermal conductivity for various aggregation morphologies with 3−4 particles.
表3 6颗粒4%体积分数时导热系数与分形维数Table 3. Thermal conductivity for various aggregation morphologies with 6 particles.
图11 分形维数与导热系数对比 (3 颗为实线, 4 颗为虚线)Fig. 11. Thermal conductivity at various fractal dimensions(3-p solid line, 4-p dash line).
现有的研究成果表明, 纳米颗粒的聚集对纳米流体的导热性有着很大的影响. 采用平衡态分子动力学模拟并通过Green-Kubo公式计算热导系数;随后, 从模拟中获取纳米颗粒的质心位置并通过Schmidt-Ott公式计算出聚集体的分形维数; 此外,还分析了单颗纳米颗粒和两颗纳米颗粒情况下的对分布函数. 本文得出以下结论:
1)纳米颗粒的聚集方式可分为紧凑聚集(RLA)和松散聚集(DLA)两种, 松散聚集情况下, 两颗纳米颗粒之间存在着纳米薄层, 而紧密聚集情况下两颗纳米颗粒直接接触;
2)纳米颗粒周围纳米薄层的厚度可以从径向分布函数中读出, 此外, 纳米薄层中的基液分子达到动态平衡, 这应该是纳米层增强传热的内在机理之一, 也说明界面热阻来解释边界层传热是不合理的;
3)纳米颗粒体积分数越大, 纳米流体的导热系数就越大, 但其增量会随体积分数的增大而减小,如果纳米颗粒的体积分数和粒径均固定, 导热系数随聚集体分形维数的减小而几乎线性增大, 且较小尺寸的纳米颗粒聚集更有利于提高纳米流体的导热系数.