基于分子动力学的板叠间传热传质微观机理的研究

Study on Microscopic Mechanism of Heat and Mass Transfer between Plates Based on Molecular Dynamics Simulation

摘要: 本文采用分子动力学计算方法从微观上对热声发动机板叠间的传热传质机理及其起振频率进行了模拟研究。在往复振荡过程中，由于温度梯度的存在，使得能量面发生了迁移，最终在冷热气体间产生密度梯度，直至密度梯度为零时完成一个周期。自激振荡伴随有能量耗散，故而温度波峰与波谷随着周期数的增加逐渐减小。随着温度比的升高，起振频率逐渐降低，这将有利于微通道板叠结构缩短谐振管的长度，实现结构的微型化。

Abstract: In this paper, the mechanism of heat and mass transfer between thermoacoustic engine plates and the starting frequency are studied by molecular dynamics method. In the process of reciprocating oscillation, due to the existence of temperature gradient, the energy surface is transferred, and fi-nally a density gradient is generated between hot and cold gases until the density gradient is zero. The self-excited oscillation is accompanied by energy dissipation, so the temperature peaks and troughs decrease with the increase of period number. With the increase of the temperature ratio, the starting frequency decreases gradually, which will help the microchannel stack structure shorten the length of the resonance tube and realize the microminiaturization of the engine.

1. 引言

近年来，由于长期使用氟利昂带来的环境问题受到了越来越多的关注，加之传统的机械制冷装置因为震动，密封要求高等对制冷机的长期稳定运行有一定的影响，制冷装置的发展遇到了瓶颈。为此，人们开始寻找新的制冷装置，以期改变这种局面。1988年Swift最早对热声制冷机进行了报道，自此，热声制冷技术引起了人们的广泛关注。B. Yu等人 [1] 测试了由行波热声发动机驱动的热声制冷机，对不同直径惯性管的制冷机进行了研究。实验发现此方法可以提高了制冷机的COP和发动机的输出效率。杨卓等人 [2] 设计了一台气液双作用行波热声发动机上使用的行波制冷机，在300 K的环境温度，250 K的制冷温度下进行理论研究，COP达到了2.74，相对卡诺效率接近60%，效率高。汪建新与张彤 [3] 设计了一种高效热声制冷机谐振管，采用Fluent软件进行模拟，发现纵向谐振频率与声源驱动器频率相近，振动增强，热声制冷机制冷效率提高。何秋石等人 [4] 建立了热声制冷微循环模型，算例模拟结果表明，在热声制冷微循环中，制冷量、制冷率以及制冷机性能系数都得到了优化解，对热声制冷机的改进及优化起到了一定的理论指导作用。汪建新等人 [5] 发现谐振管截面的几何尺寸是热声制冷机比较重要的参数，他们通过模拟得出相同激励条件下，采用圆锥形谐振管活塞震动幅值最大，声强较高。李德玉等人 [6] 通过研究发现，板叠长度和中心位置对热声转换效率影响比较大，并以空气为工质，在常温常压环境下获得了冷热端37℃的温差。蒋智杰等人 [7] 根据热声制冷机参数振荡特点，建立了不可逆热声制冷循环模型，提出了目标函数，并在目标函数取最大值时，使设备达到了最优的效果。A. C. Alock等人 [8] 通过研究发现，通过调节装置的几何形状可以改变其性能，并且在需要较少的输出量时，与传统蒸汽压缩式制冷机相比，热声制冷机效率更高。汪建新等人 [9] 基于热声效应理论和有限元法对热声制冷装置谐振腔内声压幅值和气柱活塞的振动特性进行了分析，得到了谐振腔内声强越高，热声制冷机效果越好的结论。基于以上研究，本文分别建立了室温与低温、高温与室温的微通道板叠热声转化模型，为使热声结构进一步微型化、高效化提供理论依据。

2. 热声转化

板叠作为热声制冷机的核心部件，由回热器两端的加热器与冷端换热器基于温度梯度，基于起振温度，使气体微团在板叠间进行往复振荡，从而在板叠两端形成温度梯度。图1是板叠内气团振荡移动的示意图。板叠间距的增加提高了热传递效率但增大了熵损耗，因此减小板叠间距，提高自激振荡频率，使板叠结构进一步微型化且保持较小的熵损耗与热声转化效率。本文基于分子动力学计算方法，将复杂的板叠与气团间对流换热与气团振荡过程分解分析，估算板叠间的高、低温气团的温差对自激振荡过程的频率、压力波动的影响，对板叠间无传热温差的温度场内的传热传质微观机理进行分析，并为微型板叠结构的设计提供理论基础。

Figure 1. Schematic diagram of air mass oscillation movement in the stack

图1. 板叠内气团振荡移动示意图

3. 分子动力学理论与模型建立

分子动力学模拟方法基于牛顿经典理论，从微观角度对温度场内热–声能量转换机理进行研究与分析，仅考虑两侧板叠对附近气团充分传热时刻冷、热气团的单周期混合振荡过程。图2是运行前的模型。模拟盒子的尺寸均为5.72 nm * 200 nm * 2291.2 nm (x * y * z)，初始压力为5 bar，模拟盒子左右两侧分别被区分为高温区域与低温区域，高温区域的z轴坐标范围为0 nm~1145.1 nm，低温区域范围的z轴左边坐标范围为1145.1 nm~2290.6 nm，板叠厚度维度采用周期性边界条件，以仿真板叠长度无限长，ylo、yhi、zlo、zhi维度均采用固定边界条件并设置反弹壁面。表1为运用理想气体状态方程计算的在预定温度、压力的原子个数，公式如下：

$n=\frac{pv}{RT}$ (1)

Table 1. Number and distribution of model atoms

表1. 模型原子数与分布

Figure 2. Molecular conformation before model operation

图2. 模型运行前分子构象

首先对高温与低温He气库采用无定型结构分别建立模型，之后进行模型能量最小化计算，并拼接形成完整的模拟盒子。表2是各模型的温度参数。先采用NVT系综对温度进行标定，采用Langevin控温方法调节颗粒间的摩擦阻力、黏性阻尼、颗粒随机碰撞力(耗散、波动定理)，将高温区域的气体平均温度设定为300 K，低温区域的气体平均温度设定为100 K。由温度差引导能量面的迁移，之后采用NVE系综对模型运行900万步，时间步长0.8 fs。模型采用LAMMPS (Large-scale Atomic Molecular Massively Parallel Simulation)进行并行计算，输出每个时刻的原子坐标与瞬时速度。之后将模拟盒子沿yhi维度每100 nm划分一个格子，采用Fortran编程统计每个格子内的温度、压力、速度、原子数、单位截面质量流率。

Table 2. Temperature parameters of each model

表2. 各模型温度参数

通过势能函数描述式赋予原子初始时刻位置，He原子之间采用LJ势能函数，并计算整个系统的势能。遵循牛顿经典力学，则模型中任意一原子i受到的力均为势能的梯度，并获得原子加速度，最后对牛顿运动定律方程进行时间积分，预测i原子在经过 $\tau$ 后的原子瞬时速度与位置。

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

$a_i=\frac{F_i}{m_i}$ (3)

其中 $r_i$ 为原子随时间迭代的瞬时位置； $v_i^0$ 为原子初始时刻速度； $v_i$ 为原子随时间迭代的瞬时速度。将所有原子置于模型中，初速度满足高斯分布，则系统中所有原子运动的动能总和应满足一下条件：

$K.E={\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{1}{2}{m}_i{\left(\left|v_i\right|\right)}^2}=\frac{i}{2}N{k}_B$ (4)

其中，N为总原子数；k_B为玻尔兹曼常数；T为热力学温度。由系统平均温度，赋予下一次循环的原子初始速度，并由该时刻原子的坐标球得到新的系统势能与动能，进行迭代计算，以下给出局部温度、压力的统计式：

$T_{cal}=\frac{{\displaystyle \sum m_i\left(v_{i,x}^2+v_{i,y}^2+v_{i,z}^2\right)}}{3N{k}_B}$ (5)

$P=\rho k_{B}T+\frac{1}{3v}\langle {\displaystyle \underset{i=1}{\overset{N}{\sum }}{F}_{ij}\left(r_{ij}\right)r_{ij}}\rangle$ (6)

4. 结果分析

4.1. 温度场与密度分布分析

由图3可见，往复振荡过程伴随着传热传质过程。冷端换热器、加热器附近的气体温度伴随气体工质的往复振荡随时间发生了波动。当由温度梯度引导能量面发生迁移时，两侧的气体温度发生了快速的变化，高温气团大量向冷端气团扩散，冷端附近密度快速上升，产生密度梯度。当2000 ps时，两侧气体微团温差较小时会产生较大的密度梯度，进而由密度梯度引导能量面的迁移，靠近冷端的气团大量向热端移动，直至两侧密度相同，但同时在两侧再次产生了较大的温度梯度，完成一个周期的振荡，周而复始。一个完整的周期时长约为38,000 ps，由于未对两侧反弹壁面进行温度控制，第二个周期的温差波峰与密度波峰略低于第一个周期，自激振荡伴随着能量耗散，通过控制加热器的加热温度为系统泵入能量，维持自激振。

Figure 3. Temperature fluctuation of cold end and hot end of model III (heating temperature 600 K)

图3. 模型三冷端、热端的温度波动(加热温度600 K)

4.2. 压力波动与起振频率分析

回热器内板叠结构靠近冷端换热器一侧，由于气团温差驱动的往复振荡输出压力波，其中实线代表冷段换热器一侧的压力波动，虚线代表加热器一侧的压力波动。由图4所示，在冷端换热器侧和加热器侧均捕捉到了驻波，对于微通道板叠结构，当T_h/T_c为1.33时，冷端压比为1.437193；T_h/T_c为1.667时，冷端压比为1.81916；T_h/T_c为2时，冷端压比为2.18442。随着温度比的升高，冷端可提供的压比与之成正比。由此可见微通道板叠结构的起振温度比远小于常规尺寸的温度比。

Figure 4. Pressure fluctuation of cold and hot sides with time

图4. 冷、热两侧压力波动随时间的变化

同时，随着温度比的升高，冷端与热端均维持了均匀的周期时长，且相位差为90度，但起振频率与温度比成反比。由于热声驱动器的另一个重要参数时谐振管的长度，由声波波长、频率决定， $v=\lambda f$，其中v为当地声速， $\lambda$ 为声波波长，f为起振频率。谐振管管长通常为1/2波长或1/4波长。分子动力学计算模型可为起振频率、谐振管管长提供指导，以模型三为例，计算所得模型中周期时长为38,000 ps，按比例求得当板叠长度为10 cm时起振频率为603 Hz，所需的谐振管长度为0.829米。因此，微通道板叠结构有利于进一步缩短谐振管长度，使发动机微型化。

4.3. 压力波动与单位截面质量流率分析

图5(a)是模型一靠近zlo侧100 nm处压力波与质量流率随时间的变化曲线，图5(b)是模型一中间1100 nm处压力波于质量流率随时间的变化曲线，在这两个位置压力波都是超前于质量流率的；而图5(c)是靠近zhi侧90 nm处压力波与质量流率随时间的变化曲线，此时则是质量流率超前于压力波。可以看到，在三个位置，相同时间点的质量流率相位基本是没有变化的，只是波幅在变化，这主要是因为驻波声场气团获得的总能量变化时只会引起波幅的变化，并不会引起相位的变化。而压力波从zlo侧到zhi侧变化了180˚，主要是气体原子在整个模型盒子内做往复运动，所以在到达一端时会因为能量面的迁移而转变运动方向，就会出现模型两侧压力波相位相差180˚的情况。

4.4. 压力波与速度波分析

图6为压力波与速度波随时间的变化曲线，压力波与速度波在整个运动过程中相位差为90˚ [10]，因为沿着波传播的方向，气体微团压力波和速度波的波形始终是保持不变的，变化的只有因为能量迁移而带来的波幅的变化，所以在整个模拟中，压力波与速度波相位差保持在90˚附近。相位有微小的波动是因为气体微团在运动过程中会受到张力、粘性阻力等的影响。

(a) (b) (c)

Figure 5. Distribution of pressure and mass flow rate with time at zlo side, middle and zhi in model 1

图5. 模型一zlo、中间及zhi处压力与质量流率随时间分布

Figure 6. Pressure wave and velocity wave change with time

图6. 压力波与速度波随时间变化

5. 结论

本文基于分子动力学计算方法对热声发动机板叠间微通道内气团移动过程进行了仿真，于微观层面阐述了板叠间压力波动主要由温度梯度与密度梯度交替主导能量面迁移完成热声能量转换。进而讨论了加热温度对起振频率与压力波动的影响。随着加热温度的升高，压力波动也逐渐升高，振荡周期缩短，起振频率升高。对驻波声场而言，在整个循环中，气体微团的压力波与速度波相位差保持在90˚，改变的只有波幅。

NOTES

^*通讯作者。