不同温度下脉管内微气团自振荡热力学过程的分子动力学研究

Molecular Dynamics Study on the Thermodynamic Process of Self Oscillation of Micro Air Mass in the Pulse Tube at Different Temperatures

摘要: 采用分子动力学方法模拟了脉管中微气团在300~10 K温度范围内压比为12/6 bar时的自振荡过程，获得了不同温度下自振荡频率、冷、热端最大温差与冷端最低制冷温度。研究结果表明，随着初始温度的降低，自振荡周期时间呈幂指数增长，微气团冷、热端温差及最大温降也逐渐降低，因而，脉管的最佳运行频率是由自振荡周期而决定的。此外，温度和压力波自振荡强度随自振荡进行而衰减，且温度越高衰减速度越快，因而时均最低制冷温度与运行频率相关，考虑到沿脉管的温度梯度分布，单级脉管的冷、热端温差不宜过大，否则不能兼顾冷、热端的最佳运行频率。显然，级数越多，制冷效率越高，同时，高温级宜采用高频，低温级宜采用低频，以提高制冷效率。给出了脉管降温温差、最佳工作频率随制冷温度的变化曲线。

Abstract: The molecular dynamics method is used to simulate the self oscillation process of the micro air mass in the pulse tube when the pressure ratio is 12/6 bar in the temperature range of 300 - 10 k. The self oscillation frequency, the maximum temperature difference between the cold end and the hot end and the minimum cooling temperature at the cold end are obtained at different tempera-tures. The results show that with the decrease of the initial temperature, the time of the self oscil-lation period increases exponentially, and the temperature difference between the cold and hot ends of the micro air mass and the maximum temperature drop decrease gradually. Therefore, the optimal operation frequency of the pulse tube is determined by the self oscillation period. In addi-tion, the self oscillation intensity of temperature and pressure wave decays with the self oscillation, and the higher the temperature is, the faster the decay speed is. Therefore, the average minimum cooling temperature is related to the operating frequency. Considering the temperature gradient distribution along the pulse tube, the temperature difference between the cold end and the hot end of the single-stage pulse tube should not be too large. Otherwise, the optimal operating frequency of the cold end and the hot end cannot be considered. Obviously, the higher the number of stages, the higher the refrigeration efficiency. At the same time, high frequency should be used for high temperature stage and low frequency should be used for low temperature stage to improve the refrigeration efficiency. The curves of the temperature difference and the optimal working frequency of the pulse tube with the cooling temperature are given.

1. 引言

脉管制冷机的发展已有五六十年的历史，它因冷端没有运动部件，结构简单，可靠性高等特点，在抗电磁干扰、降低振动和长寿命方面有明显优势。近年来，尤其是在航空航天和军工等领域备受研究者的青睐。由于军工和航空航天领域对制冷机效率要求较高，如何在一定条件下，获得尽可能低的制冷温度就成为了一个难点。Zhihua Gan等人 [1] 研制的单级斯特林脉管制冷机获得57 K的空载温度。刘少帅等人 [2] 采用双段惯性管调相的方式获得了48.8 K的制冷温度。朱强龙等人 [3] 在250 W的输入功率下得到了24.1 K的无负荷制冷温度。为了使系统更加紧凑，提高制冷机的制冷性能，庞晓敏等人 [4] 研制了一台两级气耦合斯特林型脉管制冷机，通过添加双向进气调节，获得了13.4 K的无负荷制冷温度。自预冷脉管是一种可以显著降低多级脉管内热损失的方法，L. M. Qiu等人 [5] 在二级脉管制冷机中获得了15.87 K的空载温度。陈六彪等人 [6] 通过在制冷机第一级使用液氮预冷，最终在第二级获得了9.6 K的制冷温度。液氦温度下工作的脉管制冷机一直是科研人员努力开发的重点，Cao Qiang等人 [7] 采用预冷的方式，在预冷温度为13.3~22 K时，获得了3.6 K的空载温度。可见，不同单级和多级脉管所能达到的制冷温度各有不同。

为了研究沿脉管轴向温度梯度形成的内在机理，本文采用分子动力学的方法模拟了300 K~10 K温度范围内微通道脉管中气体微团的自振荡过程，以研究初始温度对微气团热声振荡过程的影响。

2. 理论分析

2.1. 物理模型

图1是带有惯性管的脉管制冷机原理示意图。直线压缩机提供的压力波使得气体在经过回热器、冷端换热器、脉管、热端换热器、进入惯性管、气库后，再逆向依次返回压缩机，完成一个循环。经历多次循环后，在回热器的蓄冷作用下，脉管冷端温度不断降低，直到在脉管冷端、热端之间形成稳定的温度梯度，冷端达到某一低温值。因此，脉管制冷机获得低温的关键部件是脉管。

1-压缩机；2-回热器；3-冷端换热器；4-脉管；5-热端换热器；6-惯性管；7-气库。

Figure 1. Schematic diagram of pulse tube refrigerator with inertia tube

图1. 带有惯性管的脉管制冷机原理图

关于脉管的制冷机理，目前有表面泵热理论(图2) [8]、热力学非对称理论(图3) [9] 等。如图2所示，表面泵热理论认为在充气过程中，气体微团从1位置绝热压缩到 位置2，与壁面换热降温到T3，同时移动到位置3；然后放气过程开始，微气团绝热膨胀到位置4，经壁面吸热(制冷)回复到T1，同时位置回复到1，完成一个循环。如图3所示，热力学非对称理论认为，气体微团在流入脉管和流出脉管冷端时热力参数不对等，气体微团在离开脉管时温度低于进入脉管的温度，从而产生制冷。不同调相机构的脉管制冷机的制冷效应源于不同的热力学非对称效应的产生方式。

Figure 2. Pump heating process of gas micro cluster in pulse tube

图2. 脉管内气体微团的泵热过程

图4为复合脉管制冷原理示意图。脉管的制冷效应为沿着管壁的环形热黏层内气体的表面泵热和中心层柱状空间内气体的热力学非对称效应所产生的制冷效应之和。中心层产生的制冷量为：

${\stackrel{\dot{}}{Q}}_{cent}={\left(1-\frac{2\delta }{D}\right)}^{2}f\underset{j}{\overset{n}{{\displaystyle \sum }}}{Q}_j$ (1)

其中，δ为脉管中热黏层的厚度；D为脉管的内径；f为运行频率；Q_j为气体微元j在一个循环中的制冷量。

Figure 3. Thermodynamic process of gas cluster in the pulse tube

图3. 气体微团在脉管内的热力学过程

Figure 4. Schematic diagram of composite pulse tube refrigeration

图4. 复合脉管制冷原理示意图

$Q_j=\Delta m_j{c}_p\left(T_c-T_j\right)$ (2)

式中， $\Delta m_j$ 为在进气的半个周期中第j个时间间隔进入脉管冷端的气体微元质量； $c_p$ 为气体工质的比定压热容；T_c为制冷温度；T_j为 $\Delta m_j$ 在排气周期中气体微元离开脉管冷端时的温度。

在上述理论分析中，均提到了气体微团，但都未对微气团的热力学过程进行定量计算。本文采用分子动力学方法，通过对气体微团的压缩、膨胀振荡过程进行仿真，来获得定量的制冷温度的计算。

2.2. 分子动力学理论

分子动力学模拟(MD)是应用体系的力场及根据牛顿运动力学原理而发展起来的计算方法，是时下计算庞大复杂系统最广泛为人们所采用的方法。在含有很多个原子或分子的运动系统中，系统的能量为系统中分子或原子动能和系统总势能的总和。如式(3)：

$U=U_{VDW}+U_{\mathrm{int}}$ (3)

式中：U为原子间总势能；U_VDW为非键结范德瓦尔斯作用力；U_int为分子内部势能。

依照经典力学，系统中任意一个原子i所受的力是势能的梯度：

$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$ (4)

式中： $F_i$ 为i原子所受的力。

由牛顿运动定律得到i原子的加速度：

$a_i=\frac{F_t}{m_i}$ (5)

式中： $a_i$ 为i原子的加速度； $m_i$ 为i原子的质量。

要想预测i原子在经过时间t后的速度和位置，可以把牛顿运动定律方程式对时间进行积分：

$\frac{{\text{d}}^2}{\text{d}{t}^2}{r}_i=\frac{\text{d}}{\text{d}t}{v}_i=a_i$ (6)

$v_i=v_i^0+a_{i}t$ (7)

$r_i=r_i^0+v_i^0+\frac{1}{2}{a}_i{t}^2$ (8)

式中： $r_i$ 为i原子在时间t时的位置； $r_i^0$ 为i原子的初始位置； $v_i$ 为i原子在时间t时的速度； $v_i^0$ 为i原子的初始速度。t为时间。

本论文中，He-He原子之间采用UFF势能函数。

2.3. 分子动力学模型

由于脉管两端均有层流化原件，通常认为脉管内的流体为层流，因此，本文取图4中心处的细长微气团控制体为研究对象，如图5所示。由于中心处气团的周围均为气体，模型的上下、左右表面为气体弹性界面，第三个方向为周期边界条件，可模拟无限宽尺寸。为了观察在最佳频率下的制冷温度，采用自振荡的方式振荡，即给定初始压差，使气体自发振荡。

Figure 5. Molecular dynamics simulation model

图5. 分子动力学仿真模型

模型尺寸取为850 nm × 50 nm × 1.4303 nm (长 × 高 × 宽)，分为两个相等的区域，左边为低压区，初始压力为0.6 Mpa，右边为高压区，初始压力为1.2 Mpa，压比为2，略大于斯特林型脉管制冷机的压比。模型两侧气体分子数目按照理想气体状态方程PV = NRT来确定。左侧气体发生压缩过程，在左顶端温度升到最高，为热端；右侧气体发生膨胀过程，在右底部温度降到最低，为冷端，中部气体温度不变，为初始温度。因此，本模型可以同时完成压缩与膨胀过程，获得压缩温度和膨胀温度。需注意的是，图2~图4是追踪研究对象的轨迹，因而研究对象不变，但位置发生变化。本模型采用控制体方法，观察的位置不变，研究对象则由于气体分子的流动而发生变化。

本文对300 K~10 K温度范围内的微气团的自振荡过程进行了研究。表1是不同温度下模型的原子数。

Table 1. Number of model atoms at different temperatures under the same initial pressure of 0.6 MPa/1.2 MPa

表1. 相同初始压力0.6 MPa/1.2 MPa下不同温度时的模型原子数

模型建好后，首先，在NVT (number, volume, temperature)正则系综下运行，对模型的初始温度进行标定，并使两部分模型内的气体分别混合均匀。然后在NVE (number, volume, energy)巨正则系综下运行。整个过程无外力输入，气体处于自振荡模式，且系统对外界绝热。

3. 模拟结果与分析

3.1. 模型验证

将图5沿轴向(x轴)方向，均分为100 nm的格子，对格子内的分子进行热力学统计计算，且所有参数均为模型横截面上的平均值。每80 ps取一次数据。

图6是初始温度为100 K时微气团内沿轴向各点压力随时间的变化，显然压力呈正弦波形式，说明气体发生振荡，但由于气体内部粘滞的作用，压力振荡波逐渐衰减，但其周期几乎不变。可以看出，沿着轴向始终存在压力梯度，仅在半周期处，压力梯度瞬时为零，但压力梯度随振荡进行而衰减。

Figure 6. Pressure wave change trend of micro air mass at initial temperature of 100 K

图6. 初始温度为100 K时微气团压力波变化趋势

图7为各点温度随振荡的变化曲线，可以看出，温度亦呈周期性振荡波动，热端温度类似正弦变化，冷端温度类似负正弦变化，冷端与热端温度基本上关于初始温度对称。与压力波类似，沿轴向始终存在温度梯度，仅在周期结束时，温度梯度瞬时为零。温度波亦随振荡进行而衰减，在第一振荡周期内冷、热端温差达到瞬时最大值，制冷温度达到瞬时最低值。

图8是温度波与压力波的相位对比图，可以看出，在自振荡模式下，温度波与压力波几乎是同相位的。从机械共振角度考虑，当驱动压力波与自振压力波同相时，声功效率最高。但温度波是在半周期时达到峰值，此时，制冷温度显然不是最低值。这也许是COP与最低制冷温度不能同时获得的原因。

理想气体等熵膨胀过程终了温度可由式(9)得到：

$T_2=T_1{\left(\frac{P_2}{P_1}\right)}^{\frac{K-1}{K}}$ (9)

式中： $T_2$ 为等熵膨胀终了温度； $T_1$ 为等熵膨胀初始温度； $P_2$ 为排气压力； $P_1$ 为吸气压力；K为绝热指数。

表2是在不同初始温度下分子动力学模拟所得到的瞬时最低制冷温度与式(7)计算得到的制冷温度的对比。可以看到，在相同压比下，模拟所得最低制冷温度均略高于等熵膨胀制冷温度，这是微气团存在轴向导热的缘故。这说明本模拟虽然是在微尺度下进行的，但所得模拟结果与宏观尺寸的规律一致，因而计算结果是合理的。

Figure 7. Variation trend of temperature wave of micro air mass when the initial temperature is 100 K

图7. 初始温度为100 K时微气团温度波变化趋势

Figure 8. Phase comparison of temperature wave and pressure wave

图8. 温度波与压力波的相位对比

Table 2. Comparison of molecular dynamics simulation refrigeration temperature and isentropic expansion refrigeration temperature at different temperatures

表2. 不同温度下分子动力学模拟制冷温度与等熵膨胀制冷温度对比

3.2. 结果分析

由于脉管制冷机在达到稳定运行状态时，沿脉管轴向会形成一个稳定的温度梯度，因此，处于从冷端到热端不同轴向位置的微气团可以用不同初始温度的微气团代表。

图9是不同初始温度下微气团振荡交变流动时的瞬时温度变化曲线。可以看出从高温到低温的各微气团经历着相似的交变流动，每个微气团承担一定的温度梯度，连续不断的微气团构成了脉管的宏观轴向温度梯度，最终在冷端达到所需制冷温度，不同温度下的微气团其振荡周期不同。由于温度在随时间变化，通常以时均温度代表宏观温度。图10为不同温度下微气团冷、热端时均温度随振荡时间的变化。可以看出，各温度下达到最低制冷温度的时间不同，因而需要不同频率的外界驱动压力波。但是，对同一脉管，其驱动压力波频率是一定的。这就使得沿脉管长度的各微气团不能同时处于最佳周期，当满足了高温端，就不能满足低温端，只能采用一个适中值。还可以看出，从300 K到70 K，存在兼顾各温度下的最佳频率，50 K到10 K，也存在一个适中频率。这也许是脉管必须进行分级制冷的一个原因。

图11给出了不同温度下微气团具有最大温度梯度时的轴向瞬态温度分布。可以看出，在温度较低时，温度分布在热端和冷端具有较长的平缓线，而温度较高时，温度分布几乎是线性的。

(a)(b)

Figure 9. (a) Instantaneous temperature oscillation curve of micro air mass at different initial temperatures (100 K - 300 K); (b) Instantaneous temperature oscillation curve of micro air mass at different initial temperatures (10 K - 50 K)

图9. (a) 不同初始温度下微气团瞬时温度振荡变化曲线(100 K~300 K)；(b) 不同初始温度下微气团瞬时温度振荡变化曲线(10 K~50 K)

Figure 10. Time average temperature curve of micro air mass at different initial temperatures (10 K - 300 K)

图10. 不同初始温度下微气团时均温度变化曲线(10 K~300 K)

Figure 11. Instantaneous maximum axial temperature gradient distribution of micro air mass at different initial temperatures

图11. 不同初温度下微气团瞬时最大轴向温度梯度分布

图12表示了微气团冷、热端最大温差及最大温降随初始温度的变化曲线。可以看出，温差和温降都随着初始温度的下降而下降，尤其是在70 K以下，下降速度加快。当初始温度为10 K，压比为2时，最大降温仅为1.2 K。这解释了脉管制冷机降温速率随温度的下降而变缓慢的现象。

实际斯特林脉管制冷机的压比通常仅稍微大于1，因此微气团的冷、热端温差、温降会更小，且每个微气团不能同时达到最大时均温差，因而需要众多的微气团的连续接力才能达到整个脉管的冷、热端温差。目前，室温下经调相的单级脉管的最低制冷温度已经能达到30 K的数量级，二级脉管制冷机能达到10 K~20 K的制冷温度 [10] [11] [12]。

图13是微气团自振荡周期随温度的变化曲线。可以看出，自振荡周期随温度的下降呈指数上升。图14为根据图13换算的实际脉管运行最佳频率随温度的变化曲线。可以看到，随着温度的升高，最佳运行频率升高，但在50 K以下，随温度的下降而急剧下降。因此多级脉管制冷机在不同温级上，应输入不同频率的驱动压力波，低温级低频，高温级高频。需要注意的是，脉管制冷机的最佳运行频率除与最低制冷温度相关外，还与流动阻力相关。

Figure 12. Variation of maximum temperature difference and temperature drop with initial temperature

图12. 最大温差和温降随初始温度的变化

Figure 13. Temperature dependence curve of self oscillation period of micro air mass

图13. 微气团自振荡周期随温度变化曲线

Figure 14. Temperature dependence curve of self oscillation frequency of micro air mass

图14. 微气团自振荡频率随温度变化曲线

4. 结论

本文采用分子动力学方法研究了脉管内微气团在自振荡条件下的充放气的微观过程，探讨了不同温度下自振频率、冷热端温差及最大温降的变化。研究结果表明随着初始温度的降低，自振荡周期时间呈幂指数增长，微气团冷热端温差及最大温降也逐渐降低，这与等熵膨胀的宏观热力学规律一致。此外，温度和压力波自振荡强度随自振荡进行而衰减，且温度越高衰减速度越快，因而时均最低制冷温度与运行频率相关，高温下应采用高频，低温下应采用低频，以提高制冷效率。考虑到沿脉管的温度梯度分布，单级脉管的冷、热端温差不宜过大，否则不能兼顾冷、热端的最佳运行频率。显然，级数越多，制冷效率越高，同时，高温级宜采用高频，低温级宜采用低频。

基金项目

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

NOTES

^*通讯作者。