热子气动力学模拟方法.pdf

中国工程热物理学会 传热传质学 学术会议论文 编号：113575 热子气动力学模拟方法 胡帼杰 曹炳阳 过增元（清华大学航天航空学院，热科学与动力工程教育部重点实验室，北京 100084）（Tel：[phone]，E-mail：）摘摘 要要：随着科学技术的发展，在激光、碳纳米管等前沿技术中，傅里叶导热定律不再成立，过增元等提出了热质理论，建立热质运动方程得出了普适导热定律，但要对其进行验证传统的分子动力学模拟方法有局限性。与传统分子模拟方法以静子为模拟对象不同，本文提出了一种热子气动力学模拟方法，以理想气体中的热质即热子气为对象进行直接模拟，并对热子建立了碰撞模型，以描述热子的运动规律。文中通过对热子气的平衡态系统和非平衡态系统的模拟计算，不仅验证了热子气的状态方程还正确计算了氩气的热导率。关键词关键词：热质；热子；热子气；非傅里叶导热；热子气动力学模拟 中图分类号中图分类号：TK124 文献标识码文献标识码：A 0 引言 热学是门古老的科学，它起源于人类对冷与热现象的本质的追求。关于热的本质的认识，历史上曾有过激烈的争论1-3，它们可以主要归纳为两类：一类是热质说，认为热是一种看不见、没有质量的流体，可从高温物体流向低温物体，且数量守恒；另一类是热动说，认为热是物质组成粒子无规则运动的结果，是一种特定的能量形式。17 世纪至 18 世纪，基于热质说人们取得了一系列非常重要的发现，比如布莱克发现了“比热”和“潜热”，拉普拉斯和泊松找到了声速与比热比的关系，卡诺建立了热机效率的卡诺定律等，然而由于热质说无法解释摩擦生热等现象，19 世纪 40 年代以后，随着焦耳对热功当量的精确测量和能量守恒与转化定律的建立，热质说逐渐衰落，热动说得到了普遍承认，现代的热学理论都是基于热动说的观点发展而来的4。作为传热学中核心定律，傅里叶导热定律是一个唯象定律，它避开对热的本质的假定，基于实验结果得到传热过程的规律5。在传统应用领域，傅里叶导热定律已为大量的实验和工程实际所证实并获得广泛应用。不过，随着科学技术的发展，在激光、碳纳米管等前沿技术中，面对极低温、超快速以及超高热流密度等极端条件，傅里叶导热定律不再成立，导热过程中出现非傅里叶效应，大量学者对此进行了理论分析6、数值模拟7以及实验研究8。这种傅里叶定律的缺陷吸引人们从理论上对其改进，提出了各种的修正模型9,10。然而这些模型都只是进行数学的改进而没有深入物理本质，会产生新的悖论。因此这就需要基于物理本质提出新的热学理论来研究传热过程。过增元等从探讨热的本质出发，基于爱因斯坦质能关系式，提出了热质11-17的概念，即 2hhhEMVc=，(1)其中，Mh是热质，Eh是物体的内能，V 是体积，h是物体的热质密度；具体地，对于理想气体，建立了热子和热子气11-13的概念，并采用气体分子动理论导得了热子气的压力（热质压力）和热子气的状态方程，即 222052Bhnk Tpm c=，(2)其中，n是热子气的数密度，m0是热子对应的静子的质量，T是热子气的温度，kB=1.3810-23 J/K是波尔兹曼常数，c是真空中的光速。从热质概念出发，过增元等还通过建立热质运动方程得到了固体介质中的普适导热定律13-17，即()()0vvtmhtmhhc Tc TqqTuuuKqttxxx+=，(3)其中，22tmvKc T=是特征时间，hvquc T=是热质运动速度，K是基于傅里叶定律的表现热导率，q、cv、T分别是热流密度、比热和温度。由于热量传递本质上是热质在介质中的运动，因此基于第一性原理推导而得的普适导热定律描述的就是热量传递的基本规律。基于普适导热定律，不仅对非傅里叶现象给出了本质的解释，即由热质惯性不能忽略引起，还成功预测了热拥塞18等现象。目前，关于热质和普适导热定律的研究主要是理论分析，为了进行验证，亟需相关的数值模拟和实验研究。但是，由于热质很小且在极端条件下才能有所体现，因此实验研究难度很大，人们更多的寄希望于数值模拟的方法。不过，传统的分子动力学模拟方法是以静子为对象，在运动过程中认为对象质量不变，也即没有体现热质的概念，因此无法得到关于热质的结论。本文提出了一种热子气动力学模拟方法，以理想气体中的热质即热子气为对象进行直接模拟，并通过建立热子的碰撞模型，模拟热子的运动规律，从而直接统计得到热子气的状态信息，进而验证文献11提出的热子气状态方程以及普适导热定律。1 模拟方法 图1表示了本文提出的热子气动力学模拟方法和传统分子动力学模拟方法在模拟对象选择上的不同。图1(a)表示模拟系统中的静子，也就是传统分子动力学模拟方法的模拟对象；图1(b)表示考虑热质概念时模拟系统中的全子，即热子与静子之和；图1(c)表示模拟系统中的热子，也即热子气动力学模拟方法的模拟对象。可以看到，热子气动力学模拟方法的模拟对象就是热质本身，虽然热质与静质的质量比很小，但由于在模拟对象中不包含静质，微小的热质参数就可以很好的被捕捉到，进而也就可以通过建立热子气的运动方程，直接获取热子气的状态信息和运动规律。(a)模拟系统中的静子(b)模拟系统中的全子(c)模拟系统中的热子 图 1 模拟对象示意图 根据文献11-13的结论，热子总是附着在静子上并与其一起做热运动，当热子气中全子间发生相互碰撞时，静子质量不变，热子质量会发生变化；热子与静子具有相同的体积、运动速度和数密度。对于理想气体中的热子气，热子之间与静子之间一样没有势的作用，热子之间的相互作用（热子碰撞）是随同静子之间的碰撞实现的。因此在建立热子气的运动方程时，碰撞模型的假设就是至关重要的（这里只考虑热子之间的两体碰撞）。本文对于热子的两体碰撞，如图2所示1、2两个硬球分别表示两个热子，提出了碰撞模型假设。图 2 热子两体碰撞示意图 热子碰撞满足如下的控制方程 222201 102201 102211112222m vm vm vm v+=+，(4)01 102201 1022m vm vm vm v +=+?，(5)其中，m01=m01、m02=m02分别为热子1和2对应的静子的质量，1v?、1v?和2v?、2v?分别为热子1、2碰撞前后的速度，也即静子碰撞前后的速度。这里式（4）描述碰撞过程静子动能守恒，式（5）描述碰撞前后静子动量守恒。可以发现，模型是基于热子碰撞过程应使对应的静子满足动能守恒及动量守恒的假定。对上述碰撞模型还补充假定碰撞为简单碰撞，即碰撞产生的相互作用力平行于两碰撞体的质心位矢差，有 1212r?1v?2v?2v?1v?1v?2v?xyzo1112m vm vkr =?，(6)其中，取1、01和h1分别代表全子1、静子1和热子1，12r?是1、2之间位置矢量，k是比例系数。另外，对于理想气体的热子气系统，当热子发生碰撞时，热子运动服从上述的碰撞模型；当热子不发生碰撞时，由于热子之间没有势的作用，热子保持匀速运动。因此，在热子气系统中，热子的运动方程是不连续的，那么在热子气模拟方法中，步长也不应是恒定的，应取为每两次碰撞的时间间隔，时间间隔的确定方法与硬球碰撞模型19相同。本文以理想气体氩气中的热子气为模拟对象，分别对热子气平衡态系统和非平衡态系统进行模拟，以获得系统的状态参数和导热规律对于平衡态系统，模拟在系统中的热子经过足够多的碰撞后达到的稳定状态，统计稳定系统的状态参数，分析热子气状态方程。对于状态参数的统计公式，具体地有，系统温度为 2013NiiBm vTNk=，(7)其中，N是系统中热子的个数也即静子（氩原子）的个数，m0=6.63310-26 kg是氩原子的质量，vi是静子i的速度也即热子i的速度。系统中静子（理想气体氩气）的压强为 033ijijijiij iij iArArBBBrfmrvWpnk Tnk Tnk TVVtV=+=+=+?，(8)其中，WAr是静子间碰撞产生的维里项，ijr?是热子i、j之间位置矢量，ijf?、()ijvv=?分别是i、j碰撞产生的相互作用力和速度增量，t是模拟时长，V是系统的体积。对于系统中热子气的压强，由于其是热子对边界壁面的碰撞产生的，且热质碰撞中质量会发生变化，因此应从压强的本质出发来进行统计，即 20224022126ihhhiiixiixiim vWmdIpm nv dAdtnv dAdtnvVdAdtcc=，(9)其中，mhi、ni分别是热子i的质量、数密度，Wh是热子间碰撞产生的维里项，有()()330211336hijijijhiihiiijiiij iij iij imWrfrm vm vrvvtc t=?，(10)其中，mhi、mhi和iv?、iv?分别是热子i碰撞前后的质量和速度。基于速度分布满足Maxwell分布可以进一步推导得到系统中热子气的压强为()33022222220055226ijiiij ihBBhmrvvWnk Tnk Tpm cVm cc tV=+=+?，(11)可以看到，当忽略热子碰撞产生的维里项时，式（15）即是式（2）给出的定义。对于非平衡态系统，模拟使用NEMD方法20，通过施加热浴在系统中建立温差，经过充分驰豫统计稳定系统的温度分布和热流密度，分析系统的导热规律。模拟中，系统取周期性边界条件。参数选取上，取热子的直径等于氩原子的直径，即有=3.40510-10 m。取热子的数密度等于氩气的数密度，即有 00ArArgArpNnmm R TV=，(12)其中，Ar、pAr、TAr、Rg分别为氩气的密度、压强、温度和气体常数。计算中采用m=m0=6.63310-26 kg、和=1.65310-21 J作为质量、长度和能量的基本单位进行单位约化，约化单位为：数密度1/3=2.531028 m-3；温度/kB=119.8 K；压强/3=4.187107 Pa；时间=(m2/)1/2=2.1510-12 s，速度(/m)1/2=157.9 m/s，热导率kB/()=1.8810-2 W/(mK)。文中以约化单位为单位的物理量以右上标“*”表示。2 模拟结果及分析 2.1热子气平衡态系统 模拟系统共包含1000个热子，系统取正方体结构，周期性边界条件模拟中给定热子的初始位置分布为均匀排列，如图3(a)所示，速度为高斯分布，充分碰撞之后可以得到热子气系统中热子的位置分布情况如图3(b)所示，可以看到，系统稳定之后热子分布位置随机。(a)(b)图 3 热子气平衡态系统中的热子初始分布(a)和稳定分布(b)为了验证热子气压力ph与温度T之间的关系，参照文献12本文首先给定Ar=51 kg/m3，则有n=7.68881026 m-3，n*=0.03。温度分别取T*=1.0、1.2、1.5、2.0、2.5、3.0时对稳定后的系统进行参数统计，可以分别得到对应温度下静子（即氩气）的压力pAr和热子气的压力ph，计算结果如图4(a)所示，图中还给出了理想气体的p-T线和公式（2）给定的热子气p-T线，可以看到，当前工况条件下，遵循碰撞模型一的热子气系统静子压力随温度变化的趋势与理想气体相符，热子压力随温度变化的趋势与文献给出的热子气压力公式相符，但在数值上有偏差，不过误差均在7以内。为了进一步验证压力与温度的关系，模拟中另取了一组数密度，给定Ar=1.6238 kg/m3，则有n=2.4481025 m-3，n*=0.001，温度取值与前面相同。可以得到对应的pAr和ph的分布情况如图4(b)所示，此工况下静子p-T线与理想气体以及热子p-T线与文献公式均符合很好。对比两组工况下的压力分布情况可以看出，静子压力和热子压力与温度之间的变化趋势均符合理论公式（静子对应理想气体状态方程，热子对应公式（2），即静子压力p与温度T成正比关系，热子压力ph与温度T2成正比关系，而当数密度较大时，静子和热子压力都会与理论公式在数值上有所偏离，这是由于密度较大时，静子与静子之间，热子与热子之间的间距与静子和热子的直径相比已经不能忽略，它们之间相互作用的修正会越来越大。(a)1002003004000.02.0 x10-64.0 x10-66.0 x10-68.0 x10-61.50 x1063.00 x1064.50 x106 pAr Ideal Gas ph Thermon Gasp(Pa)T(K)n*=0.03(b)1002003004000.01.0 x10-72.0 x10-73.0 x10-75.0 x1041.0 x1051.5x105n*=0.001p(Pa)T(K)pAr Ideal Gas ph Thermon Gas 图 4 n*=0.03(a)和n*=0.001(b)时静子和热子气压力随温度的变化关系 为了验证热子气压力ph与数密度n之间的关系，文中首先给定温度T=144 K，即T*=1.2，数密度分别取n*=0.001、0.003、0.008、0.01、0.03、0.04、0.08时对稳定后的系统进行参数统计，可以分别得到对应数密度下静子（即氩气）的压力pAr和热子气的压力ph，计算结果如图5(a)所示。为了确定结论的普适性，计算中另给定温度T=300 K，即T*=2.5，数密度取值与前面相同，分别得到对应的pAr和ph，结果如图5(b)所示。可以看到，静子压力和热子压力与数密度之间的变化趋势均符合理论公式，即静子压力p和热子压力ph与数密度n均成正比关系，而当数密度较大时，静子和热子压力都会与理论公式在数值上有所偏离，这同样是由于密度较大时，静子与静子之间，热子与热子之间的相互作用的无法忽略引起的。另外，从图4和图5中均可以发现，基于文中碰撞模型的热子气系统，当数密度较大时，静子压力和热子压力均略高于基于理论公式得到的压力值，即静子和热子之间的相互作用均体现为引力，这在物体上也是合理的。(a)5.0 x10261.0 x10271.5x10272.0 x10272.5x10270.01.0 x10-62.0 x10-63.0 x10-64.0 x10-61.50 x1063.00 x1064.50 x106T*=1.2 pAr Ideal Gas ph Thermon Gasp(Pa)n(m-3)(b)5.0 x10261.0 x10271.5x10272.0 x10272.5x10270.05.0 x10-61.0 x10-51.5x10-52.50 x1065.00 x1067.50 x1061.00 x107T*=2.5 pAr Ideal Gas ph Thermon Gasp(Pa)n(m-3)图 5 T*=1.2(a)和T*=2.5(b)时静子和热子气压力随数密度的变化关系 2.2 热子气非平衡态系统 文中采用非平衡方法计算系统单向热导率，模拟系统取立方体结构，周期性边界条件，如图6所示，沿x轴方向将系统分为N份，分别在第1段和第1+N/2段通过速度修正的方法施加热浴使其分别处于固定温度TL=T0-T和TH=T0+T，其中T0为系统平均温度，T为热浴温度与平均温度的差值。模拟中取T0*分别为1.0、2.5和5.0，对应T为0.08、0.16和0.32。另取数密度n*=0.008，N分别为216、512和1000，则x向的特征长度Lx*=12(N/n*)1/3对应分别为15、20和25。运动方程积分取变步长，每次的碰撞间隔时间即为步长，计算中初始1000次碰撞使系统平衡到平均温度，然后约50000次碰撞对高低温段施加热浴，使系统达到稳态导热，其后约450000次碰撞对系统进行统计平均，得到系统的温度分布和热流密度，并用傅里叶导热定律即可得到当量热导率。统计中取模拟系统截面积S=Ly*Lz*=(N/n*)2/3，其中Ly*=Lz*=(N/n*)1/3分别是x向和y向的特征长度。低温热浴高温热浴TLTHxyzo周期性边界条件周期性边界条件 图 6 模拟系统结构示意图 图7分别给出了平均温度T0*=1.0、2.5、5.0时不同长度系统的温度分布图，可以看出它们的规律一致，在热浴附近区域有温度的波动，较大的系统波动也较大，而系统中间段温度近似线性分布，模拟中取中间线性段拟合计算温度梯度。另外，热浴处的波动随温度没有规律性的变化。我们认为规模较大的系统温度分布中热浴处的波动主要是由于驰豫时间不够，系统未达到完全稳定引起的。(a)-1.0-0.50.00.51.0115120125T0=120KT(K)x/Lx Lx=5.11nm Lx=6.81nm Lx=8.51nm(b)-1.0-0.50.00.51.0290300310T(K)x/Lx Lx=5.11nm Lx=6.81nm Lx=8.51nmT0=300K(c)-1.0-0.50.00.51.0580590600610620T0=600KT(K)x/Lx Lx=5.11nm Lx=6.81nm Lx=8.51nm 图 7 平均温度T0*=1.0(a)、T0*=2.5(b)和T0*=5.0(c)时不同长度系统的温度分布 计算得到对应静子系统氩气的热导率如图8所示。可以看到，计算工况下氩气热导率值在0.014-0.055W/(mK)之间，数量级与公认的氩气热导率值相符。不过，由于文中模拟的系统规模较小，系统尺寸均在10-9m量级，得到的结果不够准确，在之后的工作中会进行进一步的研究。1002003004005006000.0150.0200.0250.0300.0350.0400.0450.0500.055(Wm-1K-1)T0 Lx=5.11nm Lx=6.81nm Lx=8.51nm 图 8 不同长度系统在不同平均温度下的的热导率值 3 结论 为了对热质理论进行数值分析方面的验证，本文提出了一种热子气动力学模拟方法，以理想气体中的热质即热子气为对象进行直接模拟，结合假定的热子碰撞模型获得热子的运动规律。基于这种新的方法，文中分别对理想气体氩中的热子气平衡态系统和非平衡态系统进行了模拟，得到了热子气压力与温度和数密度的关系，验证了文献提出的热子气状态方程，同时正确计算了氩气的热导率。参考文献 1 John T.Heat：A Mode of Motion M.NY：D.Appleton and Company,1885.2 Mendoza E.A Sketch for a History of Early Thermodynamics J.Phys Today,1961,14：32-42.3 Mason S F.自然科学史 M.周煦良,译.上海：上海译文出版社,1980.Mason S F.A history of the Sciences M.Translated by ZHOU Xuliang.Shanghai：Shanghai Translation Publishing House,1980.4 申先甲.探索热的本质 M.北京：北京出版社,1985.SHEN Xianjia.Explore the Essence of Heat M.Beijing：Beijing Publishing House Group,1985.5 Fourier J.Analytical Theory of Heat M.New York：Dover Publications,1955.6 XU Yunsheng,GUO Zengyuan.Heat Wave Phenomena in IC Chips J.Int J Heat Mass Transfer,1995,38(15):2919-2922.7 Maruyama S.A Molecular Dynamics Simulation of Heat Conduction in Finite Length SWNTs J.Physica B,2002,323(1-4):193-195.8 Brorson S D,Fujimoto J G,Ippen E P.Femtosecond Electronic Heat-transport Dynamics in Thin Gold Films J.Phys Rev Lett,1987,59(17):1962-1965.9 Vernotte P.Les Paradoxes de la Theorie Continue de lEquation de la Chaleur J.C R Acad Sci,1958,246(22):3154-3155.10 Tzou D Y.A Unified Field Approach for Heat Conduction from Macro-to Micro-Scales.J Heat Transfer,1995,117(1):8-16.11 过增元.热质的运动与传递热质与热子气J.工程热物理学报,2006,27(4):631-634.GUO Zengyuan.Motion and Transfer of Thermal MassThermal Mass and Thermon Gas J.J Eng Thermophys,2006,27(4):631-634.12 张清光，曹炳阳，过增元.热质的运动与传递热子气状态方程 J.工程热物理学报,2006,27(6):908-910.ZHANG Qingguang,CAO Bingyang,GUO Zengyang.Motion and Transfer of Thermal MassEquation of State for Thermon Gas J.J Eng Thermophys,2006,27(6):908-910.13 过增元，朱宏晔.热质的运动与传递热子气的守恒方程和傅里叶定律 J.工程热物理学报,2007,28(1):86-88.GUO Zengyang,ZHU Hongye.Motion and Transfer of Thermal MassConservation Equations of Thermon Gas and Fouriers Law J.J Eng Thermophys,2007,28(1):86-88.14 王海东，曹炳阳，过增元.金属中的热质运动电子气的热质状态方程 J.工程热物理学报,2010,31(5):817-820.WANG Haidong,CAO Bingyang,GUO Zengyuan.Motion of Thermomass in MetalsState Equation for Thermomass in Electron Gas J.J Eng Thermophys,2010,31(5):817-820.15 CAO Bingyang,GUO Zengyang.Equation of Motion of a Phonon Gas and Non-Fourier Heat Conduction.J Appl Phys,2007,102:053503.16 过增元，曹炳阳.基于热质运动的普适导热定律 J.物理学报,2008,57(07):4273-4280.GUO Zengyuan,CAO Bingyang.A General Heat Conduction Law Based on the Concept of Motion of Thermal Mass.Acta Phys Sin,2008,57(07):4273-4280.17 过增元，吴晶，曹炳阳.热质 J.机械工程学报,2009,45(3):10-13.GUO Zengyuan,WU Jing,CAO Bingyang.Thermal Mass J.J Mech Eng,2009,45(3):10-13.18 WANG Haidong,CAO Bingyang,GUO Zengyuan.Heat Flow Chocking in Carbon Nanotubes J.Int J Heat Mass Transfer,2010,53:1796-1800.19 Allen M P,Tildesley D J.Computer Simulation of Liquids M.New York：Oxford University Press,1987.20 Hoover W G.Nonequilibrium Molecular Dynamics J.Annu Rev Phys Chem,1983,34:103-127.热子气动力学模拟方法热子气动力学模拟方法作者：胡帼杰，曹炳阳，过增元作者单位：清华大学航天航空学院，热科学与动力工程教育部重点实验室，北京 100084 本文链接：http:/