黄桂芳，詹斯琦，黄维清

（湖南大学 物理与微电子科学学院，湖南 长沙 410082）

近年来，一维非线性晶格体系的热传导性质引起了广泛的关注．人们试图利用解析或非平衡分子动力学方法从微观动力学机理的角度去解释宏观的热传导现象［1－6］．尽管至今没能得到一个完善的关于低维材料的热传导理论，但在经典一维系统的热传导研究过程中，发现许多新奇的现象，如热整流效应，界面热阻，负微分热阻等．这些现象的存在使得热量的载体——声子能像电子那样得到控制与管理，因此开辟了热二极管［7］，热晶体管［8］，热逻辑门［9］等热学器件的研究领域．其中热二极管的研究最为广泛，人们试图建立各种模型来实现材料的热整流性质［10－12］．经过大量的理论研究后，总结出热二极管模型需要满足两个必要条件：系统的非线性（简谐）及不对称结构．此外，实验方面也取得了不错的成绩，借助纳米科技的迅速发展，Chang等人［13］首次在实验中成功制备出碳纳米管固体热二极管．随后不同形态及材料构成的热二极管也相继在实验中成功制备［14－15］，这使热二极管的实际应用成为可能．

功能梯度材料在各领域中得到广泛的应用，它是一种在结构或性能上连续／准连续变化的非均质复合材料．对于功能梯度材料的研究，目前已逐步扩大到机械、电子、声学、光学、核工程、生物等领域．相比这些性质，人们对它的热学性质的研究并不够多．文献［16］最先将功能梯度材料概念引入一维非线性热传导领域，并提出一维线性质量梯度下的Fermi－Pasta－Ulam（FPU）链，通过计算实现了热整流效应．随后线性质量梯度下的Frenkel－Kontorova（FK）链及指数型质量梯度下的FPU 链也相继被提出用以构建热二极管模型［17－18］．在实验方面，Chang等人［13］也是通过采用质量梯度沉积制备出了不对称的无定型C6H16Pt固体热二极管．由此可见，功能梯度材料在低维非线性热传导领域极具应用价值．

本文利用非平衡分子动力学模拟方法，系统地研究一维非线性功能梯度材料中不对称热传导性质．在以往的研究中，改变系统的热整流方向需要进行结构重建．然而，本文的模型能通过调节热浴温差进行热整流方向的控制，即系统出现了热整流反转现象．热整流的反转是由负微分热阻引起的，并且它还与系统的体系大小及界面耦合强度紧密相关．鉴于该模型能在原地控制热整流方向，文章进一步讨论了其实际应用的可能性．

1 模型和计算方法

图1为系统结构示意图，该系统由两段质量连续变化的FK链组成，并由一弹性系数为kint的简谐弹簧连接，我们称该弹簧为界面弹簧．系统通过引入线性递减分布的质量梯度来实现不对称性（功能梯度），为方便观察，质量大小由晶格的大小表示．最左端与最右端晶格分别与温度为TL和TR的热浴接触．

系统的哈密顿量为：

式中：N为系统总粒子数；HLFK和HRFK分别为左右两段FK 链的哈密顿量．

图1 系统结构示意图Fig．1 Schematic picture of the system

式中：pi和xi分别为第i个晶格的动量和位置．取k＝a＝1，分别表示晶格间的简谐作用强度及晶格常数．第i个粒子的质量为：

固定Mmax＝20，Mmin＝1，Mmax为最左端晶格的质量，Mmin为最右端晶格的质量．基底势强度为V＝5．

模拟计算中，采用固定边界条件，并使两端的晶格与NOSE－HOOVER 热浴耦合，保持温度分别为TL和TR．为简单起见，定义系统的平均温度为T0＝（TL＋TR）／2，晶格链两端的温度差为ΔT＝TL－TR＝2T0Δ．当Δ＞0时，最左端晶格与温度较高的热浴耦合，最右端的晶格与温度较低的热浴耦合；当Δ＜0时，情况则相反．系统的局域化温度定义为格点i处晶格的平均动能：Ti＝m〈˙xi2〉．系统的局域热流定义为Ji＝k〈˙xi（xi－xi－1）〉，（〈〉均表示对时间的平均）．计算中，使用改进后的龙格库塔算法对晶格的运动方程进行积分，设置步长为0．1，并有足够长的计算时间（107）使系统最终达到非平衡稳态．当系统的热传导过程达到非平衡稳态以后，系统中流过不同格点的局域热流量相同且为一恒定值．我们定义从左到右流过系统的热流量的绝对值为J＋，当热浴反转，从右到左的热流量的绝对值则为J－，r＝J＋／J－ 为系统的热整流效率．

2 结果与讨论

2．1 不对称热传导

本文模型中，当界面耦合强度为kint＝1．0时，系统为一完整的FK 链．线性质量梯度下的FK 链已被证明具有热整流效应［17］．因此，我们先研究界面耦合强度较小时系统的热传导性质．图2为当N＝100，kint＝0．05时，系统在不同平均温度下的温度分布图．由图2 可见，系统的温度分布是不对称的，并且由于界面耦合弹簧的存在，所有情况中，系统都存在由界面热阻（ITR）引起的温度跃变．而且热浴反转后，界面热阻是不对称的，正是因为界面热阻的不对称性导致了系统的不对称热传导现象．此外，该现象还与系统的平均温度有关，温度过高或过低都不利于实现不对称热传导，因此我们固定系统的平均温度为T0＝0．1．

图2 不同平均温度下系统的归一化温度分布图Fig．2 Normalized temperature profile of the system for different average temperature

2．2 热整流反转

质量梯度下的FK 链，沿晶格质量增加方向的热流量要大于沿质量减小方向的热流量［17］．然而，当FK 链中引入一弹性系数较小的界面弹簧后，这种单向的不对称热传导现象将被打破．图3 为当kint＝0．05时，3种不同体系大小下热流量与热浴温差因子的关系图．由图3（a）（N＝100）可以看到，沿晶格质量减小方向的热流量J＋随｜Δ｜的增加几乎呈线性增长，而沿质量增加方向的热流在时达到最大值后随着｜Δ｜的增加而减小，导致了时J＋＜J－到J＋＞J－的热整流反转现象．因此系统仅在条件下，沿着晶格质量增加方向的热流量大于沿质量减小方向的热流量，而当｜Δ｜＞0．4时，则转为沿质量减小方向的热流量大于沿质量增加方向的热流量．可见，热整流方向的改变仅仅是通过控制热浴两端温差实现的．一般而言，热流量会随着热浴温差的增加而增大，然而在图3（a）中，当｜Δ｜＞0．5时出现了负导热现象，正因为它的出现直接导致了热整流的反转．并且负导热现象出现的区域随着体系的增大逐步移向更高的热浴温差区．因此，当体系大小为N＝1 000时，从J＋＜J－到J＋＞J－的热整流的反转现象出现在｜Δ｜≈0．7，如图3（b）所示．随着体系进一步增大，负导热现象最终消失，热整流反转也随之消失，因此，当N＝2 000时，系统中沿着晶格质量增加方向的热流量大于沿质量减小方向的热流量，即J＋＜J－ ．为了进一步研究界面耦合强度及体系大小对热传导的影响，图4给出了系统在以上3种体系下热流量J±与界面耦合强度kint的关系，固定热浴温差因子．所有的情况下，随着界面耦合强度的加强，J±迅速增加，并且在一定范围内呈现J±～k2int的关系，这种关系普遍存在于声子通过一较弱耦合界面的情况中．随着kint进一步增加，J±将逐渐达到饱和，二者的关系也趋近稳定．特别是，在图4（a）和（b）中还存在J＋＞J－到J＋＜J－的交叉现象，即当界面耦合强度较小时，系统沿晶格质量减小方向的热流量大于沿质量增加方向的热流量，随着耦合强度的增加，变为沿质量增加方向的热流量大于沿质量减小方向的热流量．而当体系较大时（图4（c）），系统的热整流方向不受界面耦合强度的影响，沿晶格质量增加方向的热流量始终大于沿质量减小方向的热流量，即J＋＜J－．可见，无论体系多大，当界面弹簧弹性系数较大时，均有J＋＜J－，这是因为当系统处于强耦合状态时，两段FK 链组成一个整体，成为一条完整的FK 链，与文献［17］中得出的热整流方向一致，且系统的热整流方向不随体系的改变而改变．而当系统处于弱耦合状态时，系统的热整流方向会随着热浴温差，体系大小的改变而改变，出现热整流反转现象．

2．3 热整流机理

针对以上不对称热传导现象，我们根据有效声子理论研究了其热整流机理：当FK 链中的晶格处于高温极限时，其动能足以克服基底势能的束缚，因此基底势可忽略不计并且系统可近似地看做为简谐晶格链，晶格中声子的振动频率范围为0＜ω／2π＜；当晶格处于低温极限时则被困在基底势垒中，通过线性化运动方程可得到晶格中声子的振动频率范围为

因此，当系统处于强耦合状态时，系统回归为一条完整的FK 链，不存在界面热阻，我们主要研究热浴两端晶格的声子频率：当系统左端质量较大的晶格与高温热浴接触时，晶格处于高温极限中，声子的振动频率为0＜ω／2π＜0．071．而与低温热浴接触的右端质量较小的晶格中，声子的频率为0．35＜ω／2π＜0．47．可见左右两端晶格中声子的频率不重合，不能匹配，这就意味着声子不能较顺利地从左往右输运热量，因此系统中的热流量较小．

图3 不同体系下热流量J 与｜Δ｜的关系图Fig．3 Dependence of the heat flux Jon｜Δ｜for different system size

图4 不同体系下热流量J 与kint的关系图Fig．4 Dependence of the heat flux J on kintfor different system size

当系统左端晶格与低温热浴接触时，其声子频率为0．079＜ω／2π＜0．106．而右端晶格处于高温极限，其声子频率为0＜ω／2π＜0．308．可见左右两端晶格中声子的频率能较好地匹配，更有利于声子的热输运，因此，J＋＜J－，这与我们的结果吻合．然而，当系统处于弱耦合状态时，以上结论无法解释热整流的反转．但可以很直观地看出，负导热现象直接导致了热整流反转，而负导热现象被认为是由负微分热阻所引起，因此需要考虑界面效应．正是由于弹性系数较小的界面弹簧存在，当系统右端质量较小粒子与低温热浴接触时，导致系统出现了负微分热阻，即随着系统右端的温度进一步减小，界面两端晶格中声子频率呈现由匹配到分离的现象．于是J－随着热浴温差的增加而减小，热整流方向由最初的J＋＜J－变为J＋＞J－，因此热整流反转了．由此可见，只有界面耦合强度较小情况下，才能通过调节热浴温差控制热整流的方向．

此外，热整流反转除了与耦合强度有关外，还与体系的大小密切相关．由图3和图4可知，只有在体系较小的条件下才能出现热整流反转，当体系大于一定的范围后，无论界面耦合强度如何，系统都呈现出和强耦合状态下一样的热整流性质，没有热整流反转现象．因此，图3（c）中热整流反转的消失可理解如下：界面弹簧对系统的影响随着体系的增大而减小，在较大的体系中，界面弹簧对系统的作用可以忽略不计，系统的两部分仍可被看做强耦合在一起的整体．因此，可以总结得出，热整流反转的根本原因是界面弹簧的引入，并且只有在特定的体系下，当系统为弱耦合状态时，才能通过调节热浴温差控制热整流方向．

2．4 热整流效率

值得一提的是，对于特定的体系（如N＝100），系统的热整流效率r＝J＋／J－ 可从r→0．94变为r→2．1，与目前实验中制备的固体热二极管的热整流效率较为一致．这意味着，该模型在热整流方向控制领域具有潜在的应用价值．此外，由于负微分热阻的存在以及功能梯度材料的广泛应用，该模型也是构建热晶体管，热逻辑门等热学器件不错的候选者．

3 结 论

本文研究了质量梯度下一维非线性功能梯度材料的热整流性质．该系统在特定的界面耦合强度及特定的体系大小范围内，可以仅通过控制热浴温差实现由J＋＜J－到J＋＞J－的热整流反转．引起热整流反转的根本原因是系统中存在的界面弹簧．鉴于该系统无需进行结构重建便能在原地进行热整流方向的控制，讨论了它在实验制备及应用方面的可能性．

［1］LEPRI S，LIVI R，POLITI A．Heat conduction in chains of nonlinear oscillators［J］．Physical Review Letters，1997，78（10）：1896－1899．

［2］HU Bam－bi，LI Bao－wen，ZHAO Hong．Heat conduction in one－dimensional chains［J］．Physical Review E，1998，57（3）：2992－2995．

［3］LEPRI S，LIVI R，POLITI A．Thermal conduction in classical low dimensional lattices［J］．Physics Reports－Review Section of Physics Letters，2003，377（1）：1－80．

［4］PROSEN T，CAMPBELL D K．Normal and anomalous heat transport in one－dimensional classical lattices ［J］．Chaos，2005，15（1）：015117．

［5］GAUL C，BUETTNER H．Quantum mechanical heat transport in disordered harmonic chains［J］．Physical Review E，2007，76（1）：011111．

［6］LI Nian－bei，REN Jie，WANG Lei，etal．Colloquium：Phononics：Manipulating heat flow with electronic analogs and beyond［J］．Reviews of Modern Physics，2012，84（3）：1045－1066．

［7］TERRANEO M，PEYRARD M，CASATI G．Controlling the energy flow in nonlinear lattices：A model for a thermal rectifier［J］．Physical Review Letters，2002，88（9）：094302．

［8］LI Bao－wen，WANG Lei，CASATI G．Negative differential thermal resistance and thermal transistor［J］．Applied Physics Letters，2006，88（14）：143501．

［9］WANG Lei，LI Bao－wen．Thermal logic gates：Computation with phonons［J］．Physical Review Letters，2007，99（17）：177208．

［10］LI Bao－wen，WANG Lei，CASATI G．Thermal diode：Rectification of heat flux［J］．Physical Review Letters，2004，93（18）：184301．

［11］YANG Nuo，ZHANG Gang，LI Bao－wen．Thermal rectification in asymmetric graphene ribbons［J］．Applied Physics Letters，2009，95（3）：033107．

［12］YANG Nuo，ZHANG Gang，LI Bao－wen．Carbon nanocone：A promising thermal rectifier［J］．Applied Physics Letters，2008，93（24）：243111．

［13］CHANG C W，OKAWA D，MAJUMDAR A，etal．Solidstate thermal rectifier［J］．Science，2006，314（5802）：1121－1124．

［14］KOBAYASHI W，TERAOKA Y，TERASAKI I，An oxide thermal rectifier［J］．Applied Physics Letters，2009，95（17）：171905．

［15］TIAN He，XIE Dan，YANG Yi，etal．A novel solid－state thermal rectifier based on reduced graphene oxide［J］．Scientific Reports，2012，2：523．

［16］YANG Nuo，LI Nian－bei，WANG Lei，etal．Thermal rectification and negative differential thermal resistance in lattices with mass gradient［J］．Physical Review B，2007，76（2）：020301．

［17］HU Tao，HU Ke，TANG Yi．Anharmonic on－site potential implies asymmetry thermal transport in 1D mass－graded harmonic lattice［J］．Physica B－Condensed Matter，2010，405（21）：4407－4412．

［18］SHAH TEJAL N，GAJJAR P N．Study of thermal conductivity and thermal rectification in exponential mass graded lattices［J］．Physics Letters A，2012，376（8）：438－441．