热声技术是一种比较高效的能量转换技术，其核心的热声效应能够实现热能和声能的相互转化。声波成为能量转化的媒介，避免了许多因为机械运动而产生的能量损耗[1]。
谐振管中声场的调节和回热器的制冷效果控制在驻波热声制冷机研究中占有极为重要的地位。Ceperley[2]最早将热声现象引入回热器， 并指出在热声现象中回热器起调相作用，即在回热器中调制出振荡流体工质的声波和热波之间的相位差，使热量和压力波动满足瑞利（Rayleigh）判据[3] 。此判据在热声学家Rott的《热声学》中给出，主要说明了热量波动与压力波动的向量夹角为锐角时供热会加强声振动，其为钝角时，则供热会减弱声振动[4]。随后Rott[5-7]结合《热声学》以及流体力学相关方程定性描述了热声现象。近年来，我国学者在此方向上也取得了许多成果。2017年，内蒙古科技大学对驻波热声制冷机板叠的温度分布进行了分析，描述了驻波热声制冷机板叠中高温低温交替呈现的温度梯度[8]。2020年，华中科技大学对驻波热声制冷机分层回热器的温度分布进行了模拟，分析了不同类型的回热器在不同压力和频率条件下的特性，得出了回热器的热流量与回热器的孔隙度密切相关的结论[9]。2023年，中国科学院低温工程学重点实验室对室温温区的驻波热声制冷机进行了研究，通过测试谐振腔和空腔组成的调相机构对整机系统的影响时发现，在环境温度为50 ℃的情况下，可获得整机5.8 kW的制冷量[10]。
相对于传统的单侧声源驻波热声制冷机，双声源驻波热声制冷机作为热声制冷机的一种新的结构形式，在声场调节和回热器温度控制上具有更加便捷的优势。2012年，中国科学院低温工程学重点实验室对双声源驱动的热声系统声场调制进行了相应研究，通过对双声源的调制，得到较合适的声场，对实际工作有一定的指导作用[11]。2017年，长春理工大学针对双侧激励源的相位差与谐振管长度的匹配问题进行了研究，发现双侧激励(若谐振腔长度为 170 mm，则受到的双侧激励为同相激励，若谐振腔长度为340 mm，则受到的双侧激励为反相激励)谐振腔内产生的平面驻波声场更加稳定、单色性好， 声压幅值是单侧激励的1.78倍（半波长谐振腔），甚至1.81倍（1个波长谐振腔），均接近于2倍关系[12]。
驻波热声制冷机中回热器两端产生的温差是评价其制冷效果的一个较为显著的指标，但已有的研究却很少提及谐振管两端声源间相位差对回热器两端温差的影响。本文主要基于热声学的基本原理，通过改变两端声源输入信号间的相位差使得驻波热声制冷机取得最佳的声场分布，进而探究此声场中回热器两端产生最大温差时所对应的回热器位置及所需的最小回热器长度，有效简化热声制冷机的设计过程。
1 驻波热声制冷机基本关系式
1.1 基本热声学原理
热声制冷机中回热器两端存在温差，所以回热器中气体为非等温流体。为了模拟热声制冷机谐振管中流体所处的声场及温度场情况，假设其满足线性热声理论的理想条件：（1）小幅声值扰动，二阶量忽略不计；（2）固体介质为理想刚体；（3）板叠长度远低于波长。热声系统回热器中的气体介质可被视为是由一系列往复运动的气体微团首尾相连组成，应满足质量守恒、能量守恒和动量守恒，上述守恒关系可以写为欧拉坐标系下的表达式：
[?p?t+?(ρu)=0] （1）
[ρ?p?t=ρh-?p+μ?2u] （2）
[ρCp?T?t+u?T-?p?t+u?p=k?2T] （3）
式中：[ρ]为密度，[u]为速度，t为时间[,μ]为动力黏度，h为体积力分布函数，[Cp]为定压比热容，k为热导率，[p]为压力， [?]为拉格朗日算子，[T]为温度。
气体状态方程为：
[ρ]=[ ρ]([p]，[T]) （4）
联立式（1）~式（4）可解出u、T、[ρ]、[p]。由于工质处于振动状态，故可将[u]、T、[ ρ]、[p]写成扰动加上平均值的形式：
[p=pm+Re[p1(x)eiωt]] （5）
[u=Re[u1(x)eiωt]] （6）
[T=Tm+Re[T1(x,y,z)eiωt]] （7）
[ρ=ρm+Re[ρ1(x,y,z)eiωt]] （8）
式中，[pm、Tm、ρm]分别表示压力、温度、密度的平均值，[p1、u1、T1、ρ1]分别表示压力、速度、温度、密度的一阶波动幅值，x、y、z为直角坐标，ω为角速度。
将式（5）~式（8）代入式（1）~式（4）得到基本热声线性方程组，表达式如下：
[iωρ1+ρmdu1dx=0] （9）
[iωρmu1=-du1dx+μ(?2u1?y2+?2u1?z2)] （10）
[ρmCp(iωT1+u1dTmdx)iωp1=k(?2T1?y2+?2T1?z2)]
（11）
[ρ1=ρmpmp1-ρmTmT1] （12）
1.2 压力对温度的影响
通过上述热声学关系的描述可知，气体微团在回热器中振荡时，压力和密度的波动会导致温度的波动。由热力学第二定律可知：
[Tds=de+pdv] （13）
式中：s为比熵，e为比内能，v为比体积。
对热力学状态方程求偏微商得体膨胀系数：
[β=1v(?v?T)p=-1ρ(?ρ?T)p] （14）
将式（14）代入式（13），得：
[Tds=CpdT-Tvβdp] （15）
[ds=CpdTT-(?v?T)pdp] （16）
由式（15）和式（16）可知：
[ds=(?s?T)pdT+(?s?p)Tdp] （17）
将式（16）代入式（17）有：
[Tdsdt=CpdTdt-βTρdpdt] （18）
通过线性化可得以下方程：
[s1=CpTmT1-βρmp1] （19）
[T1Tm=βp1ρmCp+s1Cp] （20）
其中，[s1]为熵的一阶波动量。
2 模拟设置
2.1 双声源的设置
本次模拟在谐振管两端选用端口输出模式边界条件，用于激励和吸收进入或离开声学模型中的波导结构（如小管道或通道）的声波。
多物理场仿真软件COMSOL中，给定的端口条件支持一种特定的传播模式。通常，只有平面波模式在需要热黏性表示的小结构中传播。与简单阻抗条件或能完美吸收反射声波的完全匹配层配置相比，端口条件为波导提供了优越的无反射或辐射条件。因此本模拟选择两个平面波端口作为激励系统的声源。
[p0=l∈DAineiφ(Slj+δlj)Pl] （21）
式中：[p0]为合成声压；l为端口总数；j为模式总数；在软件定义的边界条件D上，[Slj]为散射参数；Pl为第l个端口的模式形状（定义端口的物理量组成的物理量数组参数）；[Ain]为入射声源振幅；φ为相位；[δlj]为模式对应系数。
不同散射参数对应不同的端口，端口扫描功能会根据不同的数字端口名称（即不同的l值和j值对应不同的给定散射模块）选择端口类型。根据端口扫描功能的选择，端口模式会显示不同的声源辐射形式。本次模拟选取平面波辐射声源对应的（0,0）模式，即[Slj+δlj]=1的定频率的平面波辐射。
2.2 模拟基本原理
通过COMSOL模拟双声源的驻波热声制冷机，初始条件为两边声源的振幅、频率均相同。来自左、右两端的声波的波动方程如式（22）和式（23）所示。
[pa=Acos(ωt-φa)] （22）
[pb=Acos(ωt-φb)] （23）
式中：A为单一声源输出压力幅值，[φa]、[φb]为各声波信号的相位。
根据声波的叠加定理，谐振管中任意位置处的合成声压为：
[p0=pa+pb=pccos(ωt-φ)] （24）
其中，
[p2c=p2a+p2b+2papbcos(φb-φa)] （25）
[φ=arctanpasinφa+pbsinφbpacosφa+pbcosφb] （26）
综上可知，该位置处合成声波依然是一个相同频率的声振动，但是合成的声压幅值并不等于两列声波声压的振幅之和，而与两列声波的相位差值有关。在COMSOL中，通过调节两边声源的相位使驻波热声制冷机获得效果最佳的声场，然后代入热声学方程组求解其他物理量。
3 驻波热声制冷机模拟
3.1 模拟模型的建立
根据上述模拟分析以及热声基本关系可知，典型的热声热机由谐振管提供所需的声场，在不同波长的谐振管中，1/4波长最适用于微型热声制冷机，因为其损耗最小，且占空间小[13]。
本次模拟中，选择工作频率为250 Hz，工质为氮气，通过氮气中声速可知1/4波长为333 mm，回热器材质为不锈钢，换热器材质为铜，以上材质均可在COMSOL内置材料模块中直接选用。回热器由38根长34 mm、厚0.5 mm的板叠组成。通常，换热器的长度与气团的振荡幅值[n]有关。
[n=uω=pmωcρmsinn] （27）
其中，氮气的时均密度[ρm]和氮气中声速c可查表得到，谐振管内静压[pm]取3个标准大气压，[ω]为2πf，f为模拟选取的工作频率250 Hz。由式（27）可计算得到气团振幅为0.89 mm。文献[14]指出最佳冷端换热器长度为2n，热端换热器长度为4n，为了方便计算，取冷端换热器长度2 mm，热端换热器长度为4 mm，厚度均为0.5 mm，冷端换热器和热端换热器分别紧贴在回热器板叠的两端。
谐振管两端分别具有可控声源，即声源1和声源2，管内静压取3个标准大气压，管内热端换热器温度为333 K，管径为24 mm。建立简易驻波热声机的二维（two dimension，2D）模型和回热器内部构造模型，如图1所示。
<G:\武汉工程大学\2025\第3期\肖邦彦-1-1.tif><G:\武汉工程大学\2025\第3期\肖邦彦-1-2.tif>[（a）][声源1][板叠][谐振管][声源2][谐振管][（b）][回热器组件][ 热端换热器][ 冷端换热器]
图1 (a)2D驻波热声制冷机模型；
(b)回热器组件内部构造
Fig. 1 (a) 2D model for standing wave thermoacoustic
refrigerator; (b) Internal structure of the
regenerator component
3.2 驻波热声制冷机声场调节
在空谐振管中，通过调节两端声源的相位差，可获得不同的声场分布。参考文献[15]的取值，两边端口均设置为5 kPa振幅的平面辐射波。图2为相位差分别为0°、90°、180°时声场波形的模拟结果，其中[ε]为压力点距离左端声源的位移。
由图2可知，两端声源相位差为90°时，可在1/4波长的谐振管内产生最接近1/4波长的叠加波形，且此时的压力幅值达到最高，热声器件获得的声功也最高，可产生最优的声场条件。
<G:\武汉工程大学\2025\第3期\肖邦彦-2.tif>[0 50 100 150 200 250 300 350
ε / mm][12
9
6
3
0
-3
-6][p / kPa][相位差0°
相位差90°
相位差180°]
图2 不同相位差情况下驻波热声制冷机内的声场分布
Fig. 2 Sound field distribution in a standing wave
thermoacoustic refrigerator with different phase differences
3.3 回热器的位置选取
将两端声源信号间相位差设置为90°的条件下，在长度为333 mm的谐振管内取6个等间距点，探究板叠的最佳制冷位置。当回热器处于谐振管内不同位置时，回热器两端的温降情况如图3所示。其中,d表示板叠中心到谐振管左端的距离，ΔT表示板叠两端产生的温降。
<G:\武汉工程大学\2025\第3期\肖邦彦-3.tif>[50 100 150 200 250 300 350
d / mm][2.4
2.0
1.6
1.2
0.8
0.4
0.0][ΔT / K]
图3 不同位置处回热器两端的温降
Fig. 3 Temperature drop at both ends of the regenerator at different positions
由图3可知，当板叠中心到左侧声源的距离从303 mm变化至253 mm时，其两端温差有所下降，由于板叠处于d=303 mm位置时已非常接近右侧声源，所以其可以粗略视为右侧最高温差位置；在板叠从d=253 mm处左移至d=103 mm处期间，板叠两端温差持续增加；此后，随着板叠继续左移，其两端温差显著下降，所以板叠两端温差拐点出现在d=103 mm左右的区域。
考虑到冷、热端换热器和板叠自身的长度，可在d=（103[±]10）mm的范围内，继续探究能在板叠两端产生最大温降的板叠位置点，如图4所示。其中，α为气体工质在板叠通道中的位置。在板叠右端获得的低温反映了回热器的制冷效果。
<G:\武汉工程大学\2025\第3期\肖邦彦-4.tif>[0 5 10 15 20 25 30 35
α / mm][333.5
333.0
332.5
332.0
331.5
331.0
330.5][T / K][d=113 mm
d=103 mm
d=93 mm]
图4 回热器内部温度分布
Fig. 4 Temperature distribution inside the regenerator
由图4可知，d=113 mm时，回热器两端温差为1.81 K；d =103 mm时，回热器两端温差升为2.03 K； d =93 mm时，回热器两端温差又回落为1.23 K。可见，在给定工况条件下，当板叠的中心位置距离左侧声源103 mm时，板叠两端间产生的温降最大，即热声制冷机能获得最佳的制冷效果。
4 结 论
本文借助COMSOL模拟了两端声源间相位差对驻波热声制冷机谐振管内声场分布的影响，继而模拟了回热器组件置于管内不同位置时的温降效果，得到以下结论：
（1）波形的调节取决于两端声源间的驱动相位差，根据模拟结果可知，当两端相位差为90°时，热声制冷机谐振管内可获得最佳的1/4波长波形，与前人的实验结果相符。
（2）当回热器组件置于距离谐振管左端103 mm左右位置时，可在其两端获得最大温差。
（3）初始选取的板叠长度为34 mm，但模拟结果显示：设定工况条件下，距离热端换热器20 mm以外的回热器内温度变化极小。故该制冷机板叠长度取20 mm即可。
综上，利用COMSOL对谐振管内的声场和温度场模拟，确定回热器在管中的位置及回热器长度，显著简化了设计过程。