几何因素对微通道脉管制冷性能影响的分子动力学分析

2019-12-20祁影霞车闫瑾淡一埔刘雅丽

真空与低温 2019年6期

潘帅，祁影霞，车闫瑾，淡一埔，刘雅丽，张华

（上海理工大学能源与动力工程学院，上海 200093）

0 引言

自1964年Gifford等[1]发明脉管制冷机（PTR）以来，脉管制冷机以体积小、冷端无运动部件、运行寿命长、可靠性高、振动小等优异特性，被广泛用于低温电子器件、航空航天技术、红外超导技术、低温医学等领域。脉管制冷机的发展也从最开始的基本型逐渐发展到小孔型脉管制冷机（OPTR）[2]，双向进气型脉管制冷机（DPDR）[3]以及惯性管加气库调相形式的产生[4]。同时对于脉管制冷机内部机理的研究也在不断深入。

吴明等[5]采用数值模拟的方法对脉管的长径比进行了研究，发现脉管存在最佳长径比，在此长径比下制冷量最大。同时最佳长径比随着脉管的结构及运行参数的不同而变化。Zhi等[6]模拟了脉冲管低温冷却器在不同冷端相角下的制冷机理，研究了其关键部件气团的热力学循环。通过比较冷端换热器两侧气团的传热特性，揭示了冷端产生冷却功率的机理。结果表明，冷端相角决定了脉冲管内气体工质循环过程中的传热方向和大小。刘少帅等[7]建立阻抗匹配模型，从冷指、压缩机方面优化脉管制冷机的性能。同时，建立了30 K温区脉管冷指的一维DeltaEC理论模型，优化惯性管脉管制冷机冷指与压缩机耦合运行参数，分析回热器、双段惯性管长度和气库体积对冷指声学阻抗及压缩机声功转化效率的耦合关系。迄今关于脉管内部的仿真研究主要以计算流体力学（CFD）和数值模拟为主，均采用节点分析技术对PTR特定稳态进行研究。将PTR作为封闭体进行单元划分；每个控制单元采用能量守恒，质量守恒和动量守恒方程，用于模拟某个PTR系统的稳态。因此，不可能模拟脉管中气体的热力学不对称交替流动。

对于非平衡过程，分子动力学（MD）模拟是一种行之有效的方法。随着计算机技术的快速发展，MD方法被用于计算复杂的系统。Cao等[8]运用非平衡分子动力学模拟研究表面粗糙度对亚微米铂通道中气态氩气滑流的影响，得出速度滑移的边界条件不仅取决于克努森数，还取决于表面粗糙度。摩擦系数不仅随着克努森数的减少而增加，而且随着表面粗糙度的增加而增加。Kazemi等[9]对有机纳米通道中气体流动进行了非平衡分子动力学模拟研究。研究表明，气体分子通过有机纳米级管道的传输方式是扩散。因此建立一个可以预测回收非常规流体的实际模型，在有机微纳米孔中的传输方程应该用扩散传输方程代替。本文采用分子动力学模拟，建立了不同几何条件下的微通道脉管模型，获得脉管内分子运动的微观细节，研究几何参数变化对脉管制冷性能的影响。

1 分子动力学模拟方法

分子动力学模拟遵循经典牛顿力学定律。在研究含有大量原子的独立系统时，考虑系统总能量为系统中所有颗粒的动能（EK）与势能（U）的总和，势能分为原子间非键合范德瓦尔斯作用（VDW）与分子内势能（INT）如式（1）所示。

模型中所有粒子均为独立原子，所以忽略分子内势能。VDW近似为原子对之间VDW的加成。

式中：rij为i、j两原子之间的距离。

基于模型的精度要求，选择采用混合势力场。Fe-Fe之间采用嵌入原子方法EAM/FS（Embedded Atom Method）描述。该多体式力场的描述方法是把晶体的总势能分为原子间作用力的短程对势和原子镶嵌在电子云中的嵌入能[10]。EAM/FS势函数常用于描述金属离子之间的相互作用力，如式（3）所示：

式中：Ei为原子间的总势能；Fα为原子嵌入电子云的嵌入能；ραβ为电子密度函数；φαβ为原子核间的短程对势；rij为i、j两原子之间的距离。

文中采用的流体原子为氦原子，壁面原子为铁原子。固体、气体分子之间存在相互作用，当两个粒子间距非常小时产生排斥，当粒子相距一定距离时又会产生吸引力。He-He、He-Fe之间均采用Lennard-Jone势能，其势能表达式如式（4）所示。

式中：ε为势能参数，其值为0.000 607 098；σ为势能参数，其值为2.103；rij为i、j两原子之间的距离。

运用牛顿经典力学获得各原子的受力及加速度：

模拟运用Leap-frog方法进行积分预测获得不同时间i原子的运动位移、速度等参数：

2 模型建立及模拟方法

建立矩形截面直通道脉冲管模型，通道外为高压He、通道内为低压He，模拟脉管的充气过程；相反，将通道外设置为低压He，通道内设置为高压气体，模拟脉管的放气过程。模型建立过程：首先建立高低压He原子为流体介质的稳态模型，同时引入周期性边界条件，之后分别让高低压稳态模型在NTV（粒子数N、体积V、温度T）正则系综下运行，分别使两个稳态模型内部气体分子均匀混合；其次分别建立轴向Fe壁面和纵向Fe壁面。基于上述方法，建立了基本型脉管的充放气模型。对于充气模型，初始通道外高压1 200 kPa，管内低压100 kPa。基本型脉管模型从左向右分为两部分，分别是微通道脉管及内部低压He和通道外高压He。图1为基本脉管充气模型是（1B-2L）初始时刻的原子分布，其中1B代表模型宽573.26 nm，2L代表模型长1 145.41 nm。具体参数如表1所列。

表1 各脉管充气模型尺寸信息Tab.1

图1 基本型脉管充气模型及初始原子分布Fig.1 The basic inflation model and the initial atom distribution

放气模型是将高低压He空间对调位置，与充气模型类似。对系统在NTV正则系综下运行温度的标定，将初始温度设定为300 K，仅对通道内的He气体粒子进行NVE（原子数N、气体体积V、气体分子总能量E）计算，时间步长选用0.4 fs，每50 000步输出一次结果，运行至通道内外压力平衡。

3 模拟结果与分析

3.1 模型轴向温度、压力速度分布

基于建立的基本型脉管1B-2L充气模型，在轴向方向，从左到右进行网格划分，然后在不同时段内对每一个网格内的He原子压力、速度和个数进行统计平均得出模型的轴向速度、压力以及温度随时间的变化分布，如图2～4所示。将坐标原点定义到模型最左端，从左向右为正方向。

图2 基本型脉管充气过程轴向压力分布曲线Fig.2 Axial pressure distribution during inflation of basic type pulse tube

图3 基本型脉管充气过程轴向平均速度分布曲线Fig.3 Axial average velocity distribution of basic vessel during inflation

图4 基本型脉管充气过程轴向温度分布曲线Fig.4 Axial temperature distribution of basic pulse tube during inflation

图2为基本型脉管的充气过程轴向压力分布，初始管外压力为1 200 kPa，管内压力为100 kPa，随时间推移管外压力逐渐减小，脉管内压力逐渐增大，在1 200 ps时脉管内压力高于脉管外压力。这是由于内外的压力差，使得通道外He原子流向内部，通道内He原子增多，压力逐渐增大。通道外情况则相反。当内外压差为0时，由于惯性作用，He原子继续向通道内流动，出现通道内压力高于通道外的情况。1 200 ps之后，管内压力下降，管外压力上升，直至内外压力再平衡（平衡压力650 kPa）。图3为基本型脉管充气过程轴向平均速度分布，初始状态（管外压力1 200 kPa，管内压力100 kPa）时He原子在压差的作用下速度迅速提升，在200 ps出现最大值（434 m/s），之后峰值逐渐减小同时向管内移动。这是因为初始状态通道进口截面处（600 nm）压差最大速度提升最为迅速，同时速度的变化较压力变化会有一定的延迟，所以速度最大值是在脉管进口靠左（500 nm）处而不是在进口截面（600 nm）处，随着压力向管内传递，内外压差逐渐减小，速度在向管内传递的同时也在减小。1 400 ps出现速度负值，是因为此时通道内部压力梯度反向，使得一部分He原子反向流动到通道外直至内外压力平衡。图4为基本型脉管充气过程轴向温度分布，初始时刻设定系统温度300 K，从初始时刻到200 ps，在脉管进口截面附近，通道外温度下降，通道内温度上升，而模型两端温度还基本维持在300 K。之后系统最高温度逐渐左移，最低温逐渐右移。从原子做功角度分析，运行初期通道外He原子在压力差作用下向通道内运动，通道外原子对临近的通道内原子做功，所以通道内温度上升，通道外温度下降，而两端由于功的作用还没有传递过去，所以还维持在初始温度300 K。随着过程进行，通道外原子持续对通道内做功，功的作用持续向两端传递，使两端温差逐渐扩大，在1 200 ps时热端温度达到最高498 K，在远离封闭端处温度降低至217 K。之后热端温度下降冷端温度上升，这是由于He原子流动转向，致使轴向气体混流加剧，导致通道内外的温度差减小。微通道脉管的放气过程与充气过程类似，只是参数的变化方向相反。

3.2 微通道长度变化影响

图5和图6表示模型中脉管宽度保持不变，长度改变后，轴向温度分布的变化。可以看出，对于微通道脉管模型，随着长度增大，冷端温度逐渐降低而热端温度逐步上升，同时发现热端最高温度随着长度增大而增大的幅度减小。冷端最低温度从289.5 K（1B-1L）逐渐降低到265.2 K（1B-3L），热端的最高温度从461.8 K（1B-1L）上升到505.8 K（1B-3L），需要说明的是冷端为模型气体进口，热端为模型脉管末端。对比冷热端温度，取逐时温差，获得冷热端温度差的逐时变化如图7所示。

图5 模型中脉管长度改变后冷端温度分布曲线Fig.5 Temperature distribution of cold end in the model with the change of length of pulse tube

图6 模型中脉管长度改变后热端温度分布曲线Fig.6 Temperature distribution at the hot end in the model with the change of length of the vessel

图7 脉管长度改变后冷热端温度差的逐时变化曲线Fig.7 Time by time change of temperature difference between hot and cold ends after pulse length change

由图7可以发现，随着长度增加，最大温差增大，但增大的幅度在减小，表明当压比一定时（12∶1），存在一最佳长径比。由于脉管长度的增长，微通道脉管内部气体分子在压差作用下，运行时间增长，当脉管内外压差变为零后，由惯性力主导的脉管管外分子对脉管封闭端做功的持续位移也越长，做功增多，所以脉管冷端温度降低，热端温度增高。然而这个过程不是无限的，不能持续地增加脉管长度，脉管内气体分子也无法在惯性作用下持续对脉管封闭端做功。所以适当增加微通道脉管长度可以提高制冷性能，但也存在一个极值。

3.3 微通道宽度变化的影响

模型中脉管长度保持不变，宽度改变后，轴向温度分布变化如图8～10所示。可以看出，随着模型脉管宽度的增大，冷端温度逐渐降低而热端温度逐步上升，同时可以发现，热端最高温度随着模型宽度增大而增大的幅度在减小，冷端的最低温度从265.3 K（1B-3L）逐渐降低到244.2 K（3B-3L），热端的最高温度从505.9 K（1B-1L）上升到550.6 K（3B-3L）。对比冷热端温度，取逐时温差，获得温差逐时变化如图10所示。可知随着宽度增大，冷热端最大温差增大但增大的幅度在减小，表明对于当压力比一定时（12∶1），存在一最佳长径比。由于保持长度不变增加微通道脉管宽度，同时也增加了脉管中分子数量，相应的分子对微通道脉管热端做功越多，在脉管冷端产生的冷量也越多，脉管制冷机性能提高。然而这个促进作用也有限度，由于脉管宽度的增加，也增大了脉管内部分子运动的混乱程度，增加了微通道脉管冷热端的热传导，抑制了微通道脉管制冷机性能的提高。所以增加微通道宽度对脉管制冷机的制冷性能有促进作用，但也存在一个最佳值。

图8 模型中脉管宽度改变后冷端温度分布曲线Fig.8 Temperature distribution of cold end in the model with the change of pulse width

图9 模型中脉管宽度改变后热端温度分布Fig.9 Temperature distribution at the hot end in the model with the change of the width of the vessel

图10 脉管宽度改变后冷热端温度差的逐时变化曲线Fig.10 Time by time change of temperature difference between hot and cold ends after pulse width change

3.4 脉管体积变化影响

图11 、图12为同一长径比下改变脉管体积对微通道脉管制冷机冷热端温度的影响，即保持长径比不变，增大脉管体积（同比例放大）。

图11 同一长径比下改变脉管体积对微通道脉管制冷机冷端温度的影响Fig.11 Effect of changing the vessel volume on the cold end temperature of microchannel tube cooler at the same aspect ration

从图11、图12中可以看出，随着脉管体积的增大，脉管冷端温度逐步降低，热端温度逐步升高，同时升高与降低的幅度在减小，说明当压力比（12∶1）和长径比一定时，脉管体积存在一最佳值使得冷热端温差最大。综合微通道脉管长度变化与宽度变化的共同影响，在保持微通道脉管长径比不变的条件下适当提高脉管体积可以提高其制冷性能。这表明在微通道脉管的设计过程中适当地放大脉管体积可以提高脉管制冷机的工作性能。

图12 同一长径比下改变脉管体积对微通道脉管制冷机热端温度影响Fig.12 Effect of changing the vessel volume on the hot end temperature of microchannel tube cooler at the same aspect ration