斯特林制冷机间隙密封泄漏的分子动力学模拟

穆德富，祁影霞，孟祥麒

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

1816年，Stirling提出了一种由两个等温压缩和膨胀过程与两个等容回热过程组成的闭式热力循环，称为斯特林循环。1862年，Krik将斯特林循环的逆循环用于制冷并获得成功，该循环称为逆向斯特林循环，也称斯特林制冷循环［1］。早期的斯特林制冷机的可靠性和寿命问题限制了它的应用空间。传统的斯特林制冷机一般采用活塞环实现密封，这种密封存在磨损，因而工作寿命受到限制。由运动部件上的密封环摩擦带来的污染及密封环本身的磨损是影响制冷机使用寿命、增加泄漏的一个重要因素。直到牛津型斯特林制冷机的出现，制冷机的稳定性和寿命才得以大大提高。牛津型斯特林制冷机是20世纪80年代初，由牛津大学Davy首次在制冷机中采用线性驱动、间隙密封、摆线悬梁板弹簧支撑、工质泄漏污染等关键技术而研制成功的制冷机。

这类制冷机的压缩机活塞由直线电机驱动自由压缩活塞，采用线性驱动的制冷机在寿命上占优势，不但减少了运动部件，简化了结构，而且大幅度减少了干扰。非线性压缩机的金属部件之间使用润滑油以形成液膜用于润滑和冷却运动部件［2］，润滑油也起到一定的密封作用。而自由活塞斯特林制冷机采用板弹簧支撑的间隙密封，利用板弹簧的轴向刚度小、径向刚度大来保证活塞与气缸的间隙密封［3］。与传统的环密封相比，间隙密封轴孔两零件采用间隙配合，使两零件无接触。因此，摩擦损失减小，提高了制冷机的寿命。斯特林制冷机一般采用少量氦气作为制冷工质，整个循环无相变［4］。由于有间隙的存在，当密封两端压力不相等时会引起气体的泄漏，所以必须考虑由于气体泄漏而造成的冷量损失。

Bailey等［5］通过测量活塞移动位置来测量间隙密封的层流流动，以此估计在特定行程下的间隙密封损失；陈曦等［6］推导了活塞运动和交变压力波同时存在情况下，环形间隙的泄漏量和一个周期内的平均泄漏量的计算公式，指出泄漏由活塞振动和压差两部分组成；龚俊等［7］建立了层流工况下斯特林发动机气缸与活塞间隙密封的数学物理模型，推导了密封间隙的泄漏量，最后与理想条件下的泄漏量进行比较，得出实际工况和由位置偏心引起的密封间隙的泄漏量较理想状态下泄漏量大；马诗旻等［8］研究得出，压力波和间隙宽度对泄漏量的影响较大，随充气压力、压比以及间隙的增大，泄漏量增加；卢明［9］分析了几种形式的间隙密封的流动特性，并使用Fluent软件对间隙密封在稳态层流、不可压、定温、定黏度、内外壁面无相对滑动的条件下进行了数值模拟，模拟结果与理论推导非常接近。

过去对间隙泄漏的研究主要集中在稳定层流、不可压缩、恒定温度的条件下使用流体力学方法的数值模拟。本文提出一种通过分子动力学模拟的方法来研究间隙密封的泄漏机理，分析泄漏机理、压差对泄漏量的影响。

1 模拟的基本原理和方法

1.1 基本原理

分子动力学模拟是根据牛顿力学原理发展起来的计算方法。该方法最早由Alder［10］于1957年引入分子体系，其基本原理是通过牛顿经典力学计算物理系统中各个原子的运动轨迹，然后使用一定的统计方法计算出系统的力学、热力学、动力学等性质。在分子动力学中，首先将由N个粒子构成的系统抽象成N个相互作用的质点，每个质点具有坐标(通常在笛卡儿坐标系中)、质量、电荷及成键方式，按目标温度根据Boltzmann分布［11］随机指定各质点的初始速度，然后根据所选用的力场中相应的成键和非键能量表达形式对质点间的相互作用能以及每个质点所受的力进行计算。接着依据牛顿力学计算出各质点的加速度及速度，从而得到经一指定积分步长(通常为1 fs)后各质点新的坐标和速度，这样质点就移动了。经一定的积分步数后，质点就有了运动轨迹，同时设定一定的时间间隔对轨迹进行保存。最后对轨迹进行各种结构、能量、热力学、动力学、力学等分析，从而得到所需要的计算结果。

1.2 模拟方法

利用Materials Explorer软件模拟建立由铁原子晶体与氦原子组成的物理模型。物理模型如图1所示，图中：L为压缩腔长度；d为间隙宽度。首先，建立He制冷剂在间隙内部的周期性边界条件，在背压腔端和壁面1的右侧，设置一定容积的He，标记为He1，在壁面3的左侧空间，设置与He1同样压力的He2，壁面设置为一定空间规则排布的Fe原子。整个泄漏的模拟过程为He气体从背压腔通过宽度为d的间隙到达压缩腔，壁面1和壁面2均由Fe原子构成，其中长度、宽度、高度方向上分别设置800、5、2个晶胞。由于Fe原子的晶胞都是紧密排列的，每个晶胞为 2.860 6 Å（1 Å=10-10m），由铁原子晶胞组成壁面。该狭缝用来模拟斯特林制冷机中的间隙密封的宽度，将氦气固定在狭缝的另一端，作为气体源，通过改变狭缝一端的氦气压力来模拟密封间隙内压力变化对泄漏的影响。

图1 物理模型Fig.1 Physical model

用于Fe原子之间的势能函数是Johnson势能函数，He原子与He原子以及He原子与Fe原子之间采用UFF势能函数。模拟在NVT系综进行，温度保持在298 K，时间步长为0.8 fs，模拟以10 000步的速度执行。总的运行过程是以这个速度连续运行，将建好的物理模型上传到服务器中进行计算。运行一段时间后，统计从间隙中泄漏出的He原子数。

2 模拟结果及分析

2.1 He 原子运动轨迹分析

选取在压力p= 200 kPa，温度T= 298 K，d=2 000 Å，运行时间为 184 ps（1 ps = 1.0 × 10-12s）的He原子速度矢量图，如图2所示。从图中可看出，当d= 2 000 Å时，在靠近壁面处 He 原子的密度大于中间原子的密度，在壁面附近，边界层黏滞阻力以及He原子与Fe原子相互作用力影响He原子向压缩腔泄漏的进程，而中间的He原子由于受到的影响小，向压缩腔方向的泄漏进程变慢。从速度矢量上看，垂直于壁面的原子是不可能泄漏的，只有速度矢量与壁面呈一定夹角的原子，才有可能泄漏。

图2 He 原子速度矢量图Fig.2 Diagram of He actom velocity vector

2.2 不同压力下间隙泄漏量随时间的变化

在间隙宽度为2 000 Å，其他条件不变时，模拟了不同压力下的泄漏量，统计出在不同压力下泄漏分子数随时间的变化，结果如图3所示。

由图中可知，间隙内He原子分别在p= 200、400、600、800 kPa下开始泄漏的时刻不同，总体趋势均为泄漏量随压力的增大而增大。

图3 不同压力下泄漏分子数的变化Fig.3 Changes of the leakage molecular number under different pressures

图4为不同压力下泄漏速率的变化，其中时间段 0～8 ps以 1 代替、8～16 ps以 2 代替，依此类推作为横坐标数值。从图中可看出，密封间隙内He原子从背压腔到压缩腔的泄漏速率先递增后递减，最后稳定在一定值。随着压力增大，泄漏速率峰值增大，并且最大峰值出现的时刻也推迟。一方面，由于背压腔与压缩腔的压差增大，必然引起He原子向压缩腔运动的动能增大，泄漏速率增大。而随着时间推移，压缩腔内He原子数增加，其背压腔与压缩腔的压差逐渐减小，导致泄露速率逐渐减小。当两腔体之间压差稳定后，其泄露速率也相对平稳；另一方面，压差增大也会引起He原子之间的运动变得剧烈，相互之间的膨胀加剧，向压缩腔方向的运动也会受到影响，因此达到最大泄露速率的时间相对滞后。

图4 不同压力下泄漏速率的变化Fig.4 Changes of the leakage rate under different pressures

2.3 不同间隙宽度下泄漏量随时间的变化

在压力为200 kPa，其他条件不变时，模拟了在间隙宽度不同时的泄漏量。根据模拟结果，统计出当d分别为 2 000、3 000、4 000、5 000 Å时泄漏分子数随时间的变化，结果如图5所示。从图中可以看出：泄漏量随着时间的推移逐渐增大；不同间隙宽度的泄漏量不同，随着间隙宽度的增大，泄漏量增大。

图6为不同间隙宽度下泄漏速率的变化，其中时间段 0～8 ps以 1 代替、8～16 ps以 2 代替，依此类推作为横坐标数值。从图中可以看出，密封间隙内的He原子从背压腔到压缩腔的泄漏速率先递增后递减，最后趋于稳定。随着间隙宽度的增大，泄漏速率峰值也增大。这是由于随着间隙宽度的增大，泄漏的横截面变大，导致更多的He原子通过横截面进入间隙中，造成泄漏量增多。随着时间的推移，压缩腔内He原子数增加，压差减小，泄漏速率减小，并逐渐达到一个稳定的状态。

图5 不同间隙宽度下泄漏分子数的变化Fig.5 Changes of the leakage molecular number under different gap width

图6 不同间隙宽度下泄漏速率的变化Fig.6 Changes of the leakage rate under different gap width

3 结 论

（1）通过分析整个模拟过程发现：在壁面附近He原子的密度大于在中间位置的密度；速度垂直于壁面的原子是不可能泄漏的，只有速度矢量与壁面呈一定夹角的原子才有可能泄漏。

（2）泄漏量随着间隙内压力和间隙宽度的增大逐渐增大，所以应在保证斯特林制冷机正常运行下，选择合适的压力，尽量减小间隙宽度，减小泄漏量，提高制冷机的效率。