惯性管型脉管内气体振荡制冷机理的分子动力学模拟

Molecular Dynamics Simulation of Gas Oscillation Refrigeration Mechanism in Inertia Tube Pulse Tube

摘要: 运用分子动力学模拟的方法，研究了惯性管型脉管内振荡的热力学机理。结果表明，与基本型脉管相比，惯性管型脉管温度，压力和质量流的相位均延迟，质量流和压力波之间的相位差减小，提高了脉管的制冷效率。同时，热端附近的质量流振幅有所减小，而冷端附近质量流振幅有所增大。此外，惯性管内气体温度有所降低。仿真结果有助于了解惯性管型脉管内部振荡的热力学机理。

Abstract: In this paper, the molecular dynamics simulation method is used to study the thermodynamic mechanism of the oscillation in the inertial tube (IT). The results show that the temperature, pressure and massflow phase of the inertia tube are delayed, the phase difference between the massflow and pressure wave is reduced, and the cooling efficiency of the pulse tube is improved. At the same time, the amplitude of massflow near the hot end decreases, while that near the cold end increases. In addition, the gas temperature in the IT is reduced. The simulation results are helpful to understand the thermodynamic mechanism of the internal oscillation of the IT.

1. 简介

根据Radebaugh R.等 [1] 提出的焓流调相理论，脉管冷端压力波和质量流的相位差为0˚时，理论制冷量最大；相位差为90˚时，制冷效率最低。目前最常见的调相机构为惯性管，利用其内部阻力、惯性及容抗调节内部交变流动气体的压力和速度之间的相位关系，从而优化整机的相位关系。Luo E. C.等 [2] 通过热声理论计算了湍流状态与层流状态惯性管的调相性能，研究表明，用湍流模型来描述惯性管调相比较合适。ROCH等 [3] 采用电路类比模型对惯性管进行了理论分析，研究表明与小孔型相比，惯性管调相结构能够使整机效率提高1.6倍。Dai W.等 [4] [5] 对比了脉管制冷机纯惯性管和惯性管带气库两种情况下的整机的制冷性能，发现惯性管带气库具有更好的制冷性能。在此基础上，他们通过合理的选取惯性管尺寸，使之调相角度与惯性管带气库相近。胡剑英等 [6] 对惯性管的调相能力做了一系列的研究，并研究了超高频下的调相能力。SCHUNK等 [7] 提出了惯性管分段式部件模型，将惯性管划分为多段，每一段分别由阻性、容性和感性来表示。刘少帅 [8] 等利用热声软件DeltaEC对单段和双段惯性管加气库型调相部件的调相能力进行研究，研究表明第2段惯性管会带来压力波和体积流之间更大的相位变化，第2段惯性管的长度比第1段长时更容易达到较大的调相角度。Schunk L. O.等 [9] 通过理论分析并结合实验研究了惯性管的调相能力，并对千瓦级脉管制冷机应用时调相机构的选择进行了研究。

惯性管常用的研究方法包括实验，理论计算和CFD模拟，DeltaEC可以用来计算惯性管在交变流动下的电阻项、惯性项和电容项，并推导出了固定尺寸的进出口的相位变化。这些研究方法都是从宏观角度，采用节点分析法或者有限元网格法，基于质量守恒和能量守恒的求解方程，设定脉管的进口或出口参数，但对脉管内部热力学性质和流动过程描述不精确，与实际流动存在差异。本文运用分子动力学模拟的方法，对气体工质的振荡不作假设，减少CFD模拟中假设条件与严重耦合引起的误差，对脉管带有惯性管和气库的热力学性质和流动特性进行了分析，研究其相位变化机理。

2. 原理

2.1. 分子动力学模拟

随着计算机的迅速发展，分子动力学模拟(MD) [10] 越来越多地应用于各种体系的模拟。分子动力学主要依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。考虑到运动系统中含有N个分子或原子，系统的能量为系统中分子的动能与系统势能的总和，依照经典力学，系统中任一原子所受的力为势能的梯度：

$F_i=-{\nabla }_{i}U=-\left(i\frac{\partial }{\partial x_i}+j\frac{\partial }{\partial y_i}+k\frac{\partial }{\partial z_i}\right)U$ (1)

由牛顿运动定律可的原子的加速度为：

$a_i=\frac{F_i}{m_i}=\frac{\partial v_i}{\partial t}=\frac{{\partial }^2{r}_i}{\partial t^2}$ (2)

a为原子加速度，F为所受力，m为其质量，v和r为其速度和位置。而且

$F_i={\displaystyle \underset{j\ne i}{\sum }{F}_{ij}}=-{\displaystyle \underset{j\ne i}{\sum }{\nabla }_{ij}{U}_{ij}}$ (3)

从原理上讲，只要准确知道每个原子的对势大小，就可以用分子动力学来构造一个系统的精确图象。对于过度金属和半导体，通常描述方法为嵌入原子势(EAM/FS) [11] 的，本文使用的EAM/FS势函数来自Mendelev和Han [12] 的工作，Fe-Fe之间的相互作用势为：

$U_{ij}=G_i{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\rho }_{h,i}+\frac{1}{2}}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\displaystyle \underset{i\ne i}{\sum }{\varphi }_{ij}\left(r_{ij}\right)}}$ (4)

${\rho }_{h,i}={\displaystyle \underset{j=1,j\ne i}{\overset{N}{\sum }}{\rho }_j\left(r_{ij}\right)}$ (5)

${\rho }_i\left(r_{ij}\right)={\displaystyle \underset{k=1}{\overset{2}{\sum }}{A}_k{\left(R_k-r_{ij}\right)}^3}H\left(R_k-r_{ij}\right)$ (6)

${\varphi }_{ij}\left(r_{ij}\right)=\frac{1}{4\pi {\epsilon }_0}\frac{Z\left(i\right)Z\left(j\right)}{r_{ij}}$ (7)

U为原子对势，G为嵌入能，表示原子i嵌入其他邻近原子在该处产生的电荷密度的能量， $\rho$ 为电荷密度， $\varphi$ 为原子间相互斥能。

对于He-He和He-Fe之间的相互作用，采用lennard-jone势模型 [13] 进行描述：

$U_{ij}=4\epsilon \left[{\left(\frac{\sigma }{r_{ij}}\right)}^{12}-{\left(\frac{\sigma }{r_{ij}}\right)}^6\right]$ (8)

势能参数如表1所示。

Table 1. L-J potential energy model parameters

表1. L-J势能模型参数

系统气体温度计算为：

$T_{cal}=\frac{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{m}_i\left[{\left(v_{i,x}-v_{ci,x}\right)}^2+{\left(v_{i,y}-v_{ci,y}\right)}^2+{\left(v_{i,z}-v_{ci,z}\right)}^2\right]}}{3N{k}_B}$ (9)

其中m为原子质量，v为原子速度，N为系统中气体原子个数，k_B为玻尔兹曼常数，v_c为气体质量中心速度，计算式为：

$v_{ci}=\frac{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{v}_i}}{N}$ (10)

由理想气体状态方程计算压力：

$p=\frac{NRT}{V}$ (11)

运用LAMIMPS程序进行仿真，在三个方向上设置周期性边界条件，截断半径取盒长的一半为7.15 Å，时间步长设为0.4 fs。首先，将模型在NVT系综下运行20 ps，获得稳定的系统。然后再在NVE系综下运行，每50,000步记录一次He原子的速度和位置。

2.2. 分子动力学模型

基本型脉管和惯性管型脉管如图1所示，气体为He原子，脉管段为一段高压气体和一段的组成，左边为高压时为膨胀过程，反之为压缩过程，高低压力比为0.6 MPa/0.1 MPa。气库压力和惯性管压力均为0.35 MPa，边界设置为虚拟的反弹壁面，即He原子将在边界上反弹。He原子与Fe原子之间没有能量交换。

Figure 1. (a) Basic pulse tube; (b) IT pulse tube

图1. (a) 基本型脉管；(b) 惯性管型脉管(从左到右：气库、惯性管、脉管)

矩形通道为300 nm × 858 nm，惯性管为50 nm × 858 nm，全长3432 nm，He原子数124,423。

3. 仿真模拟结果与分析

为了进行热力学统计计算，沿x方向每隔100 nm对气体体积进行一维网格划分，原点位于脉管段的最左端，所有热力学性质是网格上的平均值，并在200 ps (50000个步骤)的时间周期内进行平均。

3.1. 惯性管型脉管模拟结果及分析

运行后惯性管型脉管如图2所示。

Figure 2. IT pulse tube after operation

图2. 运行后惯性管型脉管

图3示出了在脉管热端(100 nm)、中间端(900 nm)和冷端(1700 nm)处的温度随时间变化。可以看出，惯性管型脉管中的温度振幅低于基本型脉管中的温度振幅。这是因为在膨胀过程中，常温下的一些气体由于热膨胀流入脉管的热端，而在脉管中的热端部分气体通过膨胀流经惯性管进入气库，导致脉管中气体的温度降低一点。从能量守恒的观点出发，可以得出，使用该装置能产生的较低的冷却温度，是由于脉管热端通过惯性管将热附加到气库中的结果。因此，还可以得出结论，由于气库中的一些气体会进入惯性管，所以脉管和气库的温度会影响脉管的冷却温度。这与Kanao等 [14] 的实验结果一致，使用惯性管调相可以降低制冷温度。

Figure 3. Temperature changes with time

图3. 温度随时间变化

从图4可以看出，压力波的相位也因惯性管而延迟，而压力的振幅在后期略有降低，这意味着气体振荡因惯性管而衰减得更快。

Figure 4. Pressure changes with time

图4. 压力随时间变化

从图5可知，惯性管型脉管的质量流振幅比基本型脉管的质量流振幅小大(<900 nm)，比基本型脉管的质量流振幅小(>900 nm)。这是因为一些气体在脉管和惯性管之间振荡，使得靠近热端的质量流量降低，质量流量的降低导致回热器热负荷的降低。

Figure 5. Massflow changes with time

图5. 质量流量随时间变化

从图6可以清楚地看出，在热端(1700 nm)，质量流量和压力波之间的相位差最小，从脉管中端(900 nm)到冷端(100 nm)，相位差越来越大。也就是说，位置离热端越远，相移越弱。惯性管和气库的存在使得质量流和压力波之间的相位差减小。这与Fang L.等 [15] 研究结果一致，减小相位差可以提高制冷效果，与Mikulin E. I.等 [16] 实验结论相符。

Figure 6. Phase difference between massflow and pressure wave

图6. 质量流与压力波相位差

3.2. 惯性管型脉管特点

图7、图8和图9表示出了惯性管型脉管中各个位置的温度、压力和质量流量的时间变化。可以看出，气库中温度、压力和质量流量几乎是不变，惯性管中温度降低，其压力波振幅较小。沿脉管从热端到冷端的热力学性质变化与基本型具有相同的趋势。

Figure 7. Temperature changes with time

图7. 温度随时间变化

Figure 8. Pressure changes with time

图8. 压力随时间变化

Figure 9. Massflow changes with time

图9. 质量流量随时间变化

4. 结论

对绝热条件下惯性管型脉管内的气体振荡过程进行了分子动力学模拟。与基本型相比，惯性管型脉管内部温度更低一些，靠近冷端位置的质量流振幅有所增大，靠近热端位置的质量流增幅有所减小。惯性管能延迟温度、压力和质量流波的相位，使热端和冷端的质量流和压力波的位相差减小。这些结果与实验结果有很好的一致性。因此，从微观角度揭示了惯性管型脉管的作用机理。

基金项目

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

NOTES

^*通讯作者。