自1963年Gifford与Longsworth提出脉管制冷机（PTC）概念以来[1]，专家学者们对PTC的研究与探索也相继展开。PTC在液氮温区与液氦温区就有强劲的竞争力，因其冷头无运动部件，对磁场和机械的干扰少，可以实现多温区、高稳定性、低成本运行。脉管制冷机广泛应用于垂直芯片、航空航天、超导磁铁、红外线探测等领域。

对脉管制冷机的优化多集中于对几何参数[2,3]、运行工况[4,5]、多级间耦合方式[6]、声功回收及相位调节的优化[7,8]，以提高制冷性能与进一步降低极限制冷温度。其中脉管内气体交变振荡的压缩、膨胀的热力学过程依赖于脉管的几何参数、运行参数与调相方式，并伴随有二次流、涡流、射流等非线性流动状态，这些因素会导致声功传递不稳定，进而影响脉管制冷机的制冷效率。为探究脉管制冷机的制冷机理，1965年Gifford与Longsworth提出了表面泵热理论，将制冷原理归纳为气体工质与脉管壁在压缩与膨胀过程中自冷端向热端的能量泵送过程[9,10]。1984年，Mikulin等人基于已知的热力学过程建立了经典分析法，梳理了工作气体在完整循环过程中的有效制冷量的表达式。在1986年Radebaugh等人基于能量与制冷守恒定律建立的焓流调相理论为后期脉管制冷机的研究与发展提供了主要的理论支持与依据，阐明了脉管制冷机的理论制冷量即时均焓流，脉管内的制冷效应与脉管内压力波与质量流的相位差相关，为小孔、惯性管、双向进气脉管制冷机的性能优化提供了理论依据。之后，又陆续发展了向量分析法[11]、热声理论[12]、热力学非对称理论[13]用于表述脉管制冷机的制冷原理。

CFD数值模拟方法也被广泛用于研究PTC内的流动与传热过程，在脉管内沿轴向形成了温度梯度，回热器回收冷量。X.B.Zhang利用CFD研究小孔型脉管制冷机（OPTC），预测了可达到最低温度的预冷却器，证明了在脉管冷端的温度不稳定现象,并拟合出不可逆熵增表达式[14,15]。Y.P.Banjare等人运用CFD模拟GM型双向进气PTC内部氦气振荡流动，通过瞬态CFD模型求解流体动量和传热的纳维-斯托克斯方程，预测了脉管制冷机的制冷能力[16]。Z.H.Gan等人进一步通过数值模拟与实验验证的方法，将GM型单级PTC在载荷为7.5KW时的最低制冷温度降至10.6K[17]。QunteDai等人[18]利用CFD研究了三维脉管中的振荡流与热传递，发现脉管中心区域的瞬态温度类似于正弦波而冷、热两端近似于非正弦波，其中最低温度出现在距离冷端1mm的位置，最高温度出现在距离热端3mm的位置。此外，REGEN、DeltaEC等商业软件也被用于预测回热器蓄冷性能与相位调节能力[19,20,21]。通常是将脉管制冷机作为一个封闭的整体，采用节点分析法或有限元网格法，对脉管进口和出口的工况参数进行设定，基于质量守恒条件求解守恒方程，因而对脉管制冷机单侧封闭非平衡热力学变质量流动过程描述不精确，与实际流动存在差异，因此脉管单循环周期内形成轴向温度梯度的微观机理至今不明确。分子动力学仿真方法基于牛顿经典力学，对气体工质的振荡过程不做任何假设，减少CFD模拟中假设条件与严重耦合引起的误差，可以直观分析脉管内非线性交变流动的微观热力学过程[22,23,24]，MD方法现已应用于微通道内的流动传热研究[25,26,27].

1、交变振荡机理

如图1所示，脉管制冷机内的温度梯度在脉管内伴随氦气的往复振荡产生。在气体工质的非线性交变流动过程中，脉管的冷端温度可在回热器的配合下随制冷循环进一步降低。通常采用焓流理论计算系统的制冷量。

图1基本型脉管结构示意图.

1-压缩机2-主换热器3-回热器4-冷端换热器5-脉管6-副换热器

对于脉管制冷而言，一个完整的制冷循环过程由一个充气压缩过程与一个放气膨胀过程组成[28]：

其中，cp为定压比热容；τ为一个完整周期的时长；T为温度；为质量流率；为体积流率。显然，当制冷量为零时，系统可看作一个稳态循环过程。当给定脉管冷端压比时，脉管的冷端在半周期时，可获得充气压缩过程中的理论极限制冷温度，在理想无损状况下脉管内半周期的能量流与焓流值相等。

脉管内半周期时均焓流值由公式（2）（3）所示，其中为脉管制冷机的充气压力，为脉管截面积，为脉管内压力波，为脉管内压力振幅。气体工质的速度波、压力波随时间呈周期性变化，则制冷量与气体工质压力与速度波的相位差相关。基于脉管内气体的瞬时流动情况，采用分子动力学仿真模拟方法可仅设置制冷循环充气初始时刻的工况，通过迭代过程研究脉管内气体非线性流动的热力学过程，明确脉管长径比、压比、运行频率对理论极限制冷温度、焓流的影响机理。

2、分子动力学模拟方法与模型

分子动力学模拟遵循经典牛顿力学定律，被广泛应用于研究复杂、分子数量庞大的微观体系[29,30]，且开发了可适用于金属、非金属、高分子聚合物、有机生化分子等模型的力场，极大的提高了仿真结果的精确度,利用MD方法可仿真脉管制冷机中氦原子逐个时刻的运动轨迹与参数，并对模型进行模块化分析。

2.1分子动力学算法与模型建立

在研究含有大量原子的独立系统时，考虑系统总能量为系统中所有颗粒的动能（EK）与势能（U）的总和，势能分为原子间非键结范德瓦耳斯作用（VDW）与分子内势能（INT）[31]。

其中rij为原子i与j的距离。通过势能函数描述式可赋予原子初始时刻位置，计算整个系统的势能。遵循牛顿经典力学，则模型中任意一原子i受到的力,均为势能的梯度，获得原子加速度，最后对牛顿运动定律方程进行时间积分，预测i原子在经过后的原子瞬时速度与位置,为原子随时间迭代的瞬时速度。初速度满足高斯分布，N为总原子数；KB为玻尔兹曼常数；T为热力学温度。由系统平均温度，赋予下一次循环的原子初始速度，并由该时刻原子的坐标求得新的系统势能与动能，进行迭代。以下给出局部温度T、压力P、流速v、单位截面流率表达式，其中以单位时间内通过截面的原子个数表征质量流率,为t时刻第i个格子内的原子总数。

为研究He气经过层流化元件后进入脉管的交变振荡，采取将实际脉管同比例缩小的方法，对脉管段进行孤立建模（如图2所示）。脉管管壁(z维度)由两个平行铁壁面代表，单块铁壁面x、y、z方向分别由5层、2层、1000层铁晶胞组成，每块平行板包含20000个Fe原子，采用边界滑移条件，将Fe壁面设置为反弹壁面，He原子在触碰壁面时仅改变速度方向为反向，不改变数值。因此模型三维方向均采用周期性边界条件，以实现通道厚度与气库的无限大，实际脉管多选用远小于通道直径的薄金属管材，且主要制冷量来自远离边界层的中心流体振荡过程，因此分子动力学模型描述了微通道脉管中心层微元轴向的往复震荡过程。为防止高压气体在周期性边界下由高压气库由边界向脉管内倒流，在通道的热端边界放置一个与脉管宽度相同的竖直Fe壁形成脉管的封闭端。高压气库尺寸与脉管相同，脉管内气体初始压力为0.1MPa，高压气库内的压力为1.2MPa。

图2：脉管仿真模型

为研究脉管长宽比对极限制冷温度的影响，建立优化脉管长径比的对比模型，其中各模型y维度尺寸均为1.4303nm，模型名称中D表示脉管管径（z方向：100nm），L表示脉管长度（x方向：286nm），其中脉管内He原子1023个、高压气库内12420个。文中各模型尺寸与原子数列于表一中：

表1：模型尺寸（单位：nm）

2.2力场与势能函数描述式

在分子动力学仿真过程中，选择合适的力场与势能函数来描述原子之间的受力情况十分关键。找到一个能同时满足不同物质间的作用方式又不牺牲精度要求的力场实属不易。基于本模型的原子类型与精度要求，采用混合力场。Fe-Fe之间采用嵌入原子方法EAM/FS描述。文中所用EAM/FS势函数参数采用2007年Mendelev与Han[32]的工作。该三体式力场的描述方法是把晶体的总势能分为原子间作用力的对势与以原子核代表多体间相互作用力镶嵌在电子云中的嵌入能[33]。而EAM/FS是Finnis和Sinclair在1984[34]年基于EAM多体式发展的将前如能函数优化为平方根形式的势函数表达式，且优化后仅为原子i、j的函数，即对于同一位置上的原子来说，不同原子对其总电子密度的贡献是不相等的。EAM/FS势函数常用于描述金属离子之间的相互作用力,如公式（11）所示：

He-He、He-Fe之间均采用Lennard-Jone势能，参数表征了L-J分子之间的相互作用强度，参数表征L-J分子直径，是长度尺寸的单位。其势能表达式以及不同原子间的计算式如式（12）所示：

2.3模拟方法

首先采用NVT系综对系统全局进行温度标定，运用Velocity-Verlet蛙跳法求解牛顿运动积分方程，控温方法选用Langevin方法，时间步长选取0.8fs，运行10万步，将全局温度设定为300K，使系统全局充分驰豫并达到温度恒定。然后采用NVE系综分别对模型内的He、Fe原子进行时间积分，时间步长选用0.4fs，采用16个CPU并行运算的方式，至内外压力平衡。

3、模拟结果及分析

模型对脉管的半周期充气压缩过程进行了仿真，脉管与高压气库的压差是气体流动的驱动力。为实现对脉管内氦原子的统计分析，将模型沿x轴一维方向每100nm划分一个微元格子，在每个统计时刻（200ps）计算格子内的局部温度、压力、速度、质量流率。设定系统的初始温度为300K，截断半径取rc=7.15Å(rc=0.5rmin)。随着NVE系综的运行，利用通道内外压差推动高压气库的He工质对管道内的气体进行压缩。随着气体工质的不断扩散，脉管截面的压力、温度、质量流量基本分布均匀，以下对轴向参数进行分析，压力波实现了轴向压力梯度、速度梯度、温度梯度的建立，在脉管的封闭端出现热端，脉管的开口端出现冷端。

3.1脉管内轴向热力学性质分布

图3、4分别为脉管内瞬时轴向速度、温度的分布曲线。伴随原子的扩散作用，脉管封闭端的温度快速升高、冷端温度迅速降低，在脉管内形成轴向温度梯度。在运行初期，脉管内的最大速度出现在脉管入口端，且速度方向由高压气库指向脉管封闭端，400ps时运动速度达到最大值随后降低，且出现最大速度的位置逐渐向热端漂移，当最大速度出现在距离脉管热端400nm处时发生反向漂移，在移动到距离脉管入口的400nm的高压气库内时再次发生漂移，这是因为脉管的气体流动是由压差驱动的缘故。脉管内的气体工质发生反流现象，这与压力变化情况相吻合，标志着充气压缩过程完全结束。因此，在脉管内的最大流速点并不位于稳定位置，而是发生持续往复移动。1200ps时获得脉管两端最大温跨（脉管冷热端温差）329.49K，但热端温度到达峰值后下降速度缓慢，至4000ps压力稳定时刻仍未降至初始温度，脉管的入口端在2400ps时获得最低值，相较封闭端极限温度的形成存在时间差。此时温跨相较1200ps时降低5.7K，且冷端在4000ps时，温度回升至初始温度。这是因为冷端受高压气库回流气体的影响温度快速回升，而热端为封闭端，仅由气体向冷端的绝热膨胀产生降温效果，因此温度下降速度较慢。因此，在压比较大、高频率工况下，微通道脉管（μm级）可在常温常压区实现制冷效应。

图3：脉管内瞬时轴向速度分布变化曲线

图4：脉管内瞬时轴向温度分布变化曲线

3.2微元内热力学性质随时间的变化

由图5、6所示，3D×3L模型内质量流率一个完整周期为3000ps，微元内压力随时间的变化曲线较平缓亦呈类似正弦分布，且振荡周期与质量流率相同。2D*3L模型的周期为3500ps，1D*3L模型的周期为3750ps，由此可见小管径的脉管压力波振荡频率更低。冷端、热端的温度曲线近似于正弦变化，存在温度扰动。

图5：脉管内各截面质量流率随时间的变化曲线

图6：微元内平均温度随时间的变化曲线

如图7、8所示，脉管冷端的压力波与速度波均与时间呈正弦线性变化，这与CDF模拟结果吻合良好[35,36]。脉管内温度波随时间的变化呈半正弦线性变化，这与压力波、速度波随时间的变化曲线不同，脉管热端温度波的振幅较冷端更大，造成了制冷机效率低下，理论比卡诺效率仅为24%。

图7：脉管冷端压力波与速度波随时间变化曲线

图8：脉管冷端温度波与质量流随时间变化曲线

3.3脉管管径的影响

如图9所示，随脉管直径的增大，热端100nm处温升相差小于10K；冷端700nm处温降相差小于12K，温跨基本维持不变，脉管两端的温跨不受脉管直径的影响。热端在600～800ps时出现最高温度；冷端在1200时出现无负荷最低制冷温度。热端出现最高温度的时间比冷端出现最低温度的时间要早。这是由于压力波在充气压缩过程的半周期时刻首先到达脉管的封闭端，形成最高温度；而最低温度则是在压力波返回时，即放气膨胀的终了在开口端形成。当脉管长度为2L固定值，管径从1D至4D变化过程中时，高低温形成的时间差均为600ps，由此可见，脉管内温度波周期与脉管的直径无关。

图9：产生极限温度时刻的脉管轴向温度分布

图10：脉管热端温度随时间的变化曲线

图10所示为脉管热端100nm处温度随运行时间的变化。脉管热端的温度在一个压缩周期内随时间的推移先增大，后减小。当脉管长度不变时，脉管直径越大，且温度波出现第一个波峰的时间也越早，且当脉管管径较细时（100～200nm），时间差越明显；而小直径的脉管的温度波曲线较为平滑，而随着脉管直径的增大，热端的温度达到最高值后下降速率越快。脉管热端温度的变化速率比冷端快，导致冷端的低温气团在未完全释放冷量时受到热端反流的高温气团的扰动，导致冷量损耗。遵循焓流调相理论的表述，当脉管冷端质量流与压力波的相位差为30°时，制冷机的效率较高[37]。

3.4脉管长度的影响

由图11可知，随着脉管直径的增加，质量流率的振幅随之增大，但相位不变。这是由于截面积增大，制冷工质量增多引起的；压力波的振幅随直径的增大而降低，导致脉管入口端压比减小，不利于冷端的制冷。且伴随脉管直径的增大，压力波的第一个波峰相位出现了微弱的滞后性，但在第二个波峰处自动修正为同相，说明脉管管径对压力波与质量流的相位差没有影响，即在2D×2L的脉管冷端质量流领先压力波54°，这与现有的实验数据吻合良好[38]。

图11：脉管冷端压力波与质量流随时间的变化曲线

图12：脉管内轴向温度、压力分布曲线

图13：脉管冷端质量流率、压力随时间变化曲线

图12所示为在脉管热端出现最高温度时轴向的温度、压力分布、质量流率。此时，随脉管长度的增加，热端的压力值降低，冷端的压力值升高，轴向压力梯度增大，导致压力波的振幅变大，这表明更长的脉管将会在热端对更多的能量，根据能量守恒原理，则冷端的温度会进一步降低。伴随压力分布的变化，在管径为200nm时，2L、3L、4L的温跨分别为252K、257K与259K，因此，在前半周期内，更长的脉管长度显示出了更大的温差，但在计算模型范围之内其增长率降低，说明存在一个最优的长径比可使脉管两端的温跨达到最大值。

图13所示脉管内气体振荡频率随长度的增加而降低，而压力波的相位随长度的增加显示出了滞后性，因此合理设计脉管的长度可以在保证单位时间质量流量的情况下，调节脉管内质量流与压力波的相位差以提高制冷效率。以上讨论均针对脉管制冷循环的压缩充气过程，根据能量守恒定理可知，放气膨胀过程可在冷端获得相较压缩过程更低的温度，温跨进一步增大，证明脉管冷端的极限制冷温度受脉管长度的影响，温度波的周期频率随脉管长径比的增大而降低。因此，对一定长度的脉管制冷机存在压力波的共振频率使制冷机可获得最优制冷效率。脉管内压力波的振荡频率与冷头长度存在线性正比变化趋势，与压比呈线性反比变化趋势。分子动力学的计算结果基于绝热、滑移边界层模型，忽略了壁面对王府振荡气体摩擦与温度的影响，可通过计算微通道内气体自振荡的周期与频率预测指导微通道脉管制冷机的最优运行频率，例如：根据3D*3L模型的计算结果（图10），压力波周期为4600ps时,对于冷头长度为220mm，压比1.2的斯特林脉管制冷机冷头最优频率为69.92Hz，与实验测得冷头最优频率误差约2.9%，与实验结果吻合良好[39,40]，为脉管制冷机预测冷头共振频率提供了数值模型。

4、结论

本文基于分子动力学模拟计算方法，对孤立脉管内气体工质的压缩、膨胀过程进行了仿真，并对脉管长径变化对轴向温度分布、气体振荡频率的影响进行了分析。脉管内的气体工质伴随充气压缩过程出现了轴向温度梯度、速度梯度。当脉管长度一定时，热端达到最高温度的时刻与冷端获得最低温度的时刻存在固定时间差，且与脉管直径无关，但随脉管长度的增加而延长。冷热两端的最大温度跨度亦不受脉管直径影响，但对应不同直径的脉管，存在一个最佳长径比可获得最大温度跨度，且温度梯度的变化趋势与压差驱动力相同但存在相位差。

在一个完整的振荡过程中，脉管的热端、冷端分别出现了随时间呈正弦线性变化的压力波、速度波与质量流率曲线与呈半正弦变化的温度波。随脉管直径的增加，热端在达到最高温度后降温速度更快，即热端的气体微团运动到冷端的时间更短，对冷端的扰动更大，且气体在脉管内的振荡周期更短。随脉管长度的增加，脉管的温度跨度增大，交变振荡周期延长并呈线性变化，最终对冷头的最优共振频率给出了数值模型。

参考文献：

[2]陈爽,杨森,黄琦,蒋珍华,吴亦农.脉管制冷机脉管热端层流化丝网的实验研究[J].低温与超导,2019,47(5):12-16..

[3]吴明,何雅玲,陶文铨,陈钟颀,潘雁频,王小军.最佳脉管长径比优化设计的理论分析及数值模拟[J].西安交通大学学报,2002(3):230-232.

[7]杨森,黄琦,蒋珍华,吴亦农.脉管内射流对脉管制冷机性能影响的仿真和实验研究[J].低温工程,2018(2):30-34+68.

[19]刘少帅,蒋珍华,张安阔,唐振刚,丁磊,吴亦农.空间红外探测用30K单级脉管制冷机高能效研究[J].红外与毫米波学报,2018,37(4):403-410.

[21]刘少帅,张华,张安阔,陈曦,吴亦农.80K脉管制冷机惯性管调相机理及优化研究[J].制冷学报,2016,37(5):100-105.

[22]王禹贺,祁影霞,车闫瑾,陈曦,张华.基本型脉管内气体振荡制冷机理的分子动力学模拟[J].制冷学报,2019,40(1):71-78.

[33]韩悌信,曾祥国,黄了,肖雨果.管道钢微观变形失效机理的分子动力学研究[J].四川大学学报.2016,53(4):865-874.

[39]邓伟峰,蒋珍华,刘少帅,张安阔,吴亦农.高温区大冷量脉管制冷机优化设计与实验特性[J].化工学报.2019,70(1):107-115.

[40]王乃亮,赵密广,陈厚磊,梁惊涛.72HZ大冷量空间脉冲管制冷机[J].科学通报.2018,63(36):3968-3972.

车闫瑾,祁影霞,潘帅,王禹贺,张华.微通道脉管非线性交变振荡的分子动力学微观机理研究[J/OL].计算物理:1-11[2020-06-28].

基金:上海市动力工程多相流动与传热重点实验室项目(13DZ2260900).