采用管束式热声核心的新型驻波热声发动机的高保真数值模拟

《International Journal of Heat and Mass Transfer》：High-fidelity numerical simulation of novel standing-wave thermoacoustic engines with a tube-bundle thermoacoustic core

　　**摘要**
热声发动机（TAEs）能够在无运动部件的情况下实现能量转换，具有结构简单、可靠性高和长期运行稳定性等优点。然而，传统的热声核心（TACs）通常采用分离的堆栈和热交换器，这限制了其紧凑性。本研究探讨了一种管束式TAC，将热交换器（HHX）、堆栈和冷交换器（CHX）集成到一个统一的结构中，作为一种更为紧凑的替代方案。开发并验证了一种轴对称非连续Galerkin（DG）模型，用于求解相关的非线性流动、热量传递和能量转换过程。研究了管半径与粘性穿透深度之比（??/????）在2.73–5.46范围内的影响。最大压力幅度出现在??/????=4.37–4.91时。尽管瞬时轴向速度剖面的空间分布相似，但最大速度值出现在??/????=3.28时，此时压力幅度并非最高。时间平均速度幅度场和涡度场进一步揭示了声流强度和膨胀诱导涡旋的发展变化。较小的管半径将高温区域限制在靠近壁面的地方，而较大的管半径则增强了径向热传导，从而提高了堆栈温度并使得轴向温度梯度更加均匀。时间平均壁面热流密度分布在各研究案例中仅有微小差异，但一致表明热能存储集中在堆栈的上游部分。正的局部能量密度主要集中在HHX-堆栈接口附近，而负值则出现在靠近壁面的粘性边界层以及TAC两端的膨胀区域内，这是由于粘性和涡流的耗散作用。最高循环功转换率出现在??/????=4.37时，而??/????=4.37–4.91的范围在热声功产生、局部耗散和热输入之间提供了良好的平衡。这些结果证明了所提出的集成TAC的可行性，并为紧凑型热声发动机的设计提供了指导。

**1. 引言**
热声发动机（TAEs）是一种利用热声效应通过可压缩流体的自激振荡将热能转换为声能的能量转换装置[1]。典型的TAE由一个声学谐振器和一个热声核心（TAC）组成。TAC包括一个热交换器（HHX）、一个冷交换器（CHX）以及位于它们之间的堆栈，以建立温度梯度[2]，而谐振器则充当声波传播的波导。由于没有运动部件且使用环保的工作流体，TAEs表现出高可靠性和耐用性，在能量转换[3]、[4]、废热回收[5]、[6]、[7]、制冷[8]、[9]和航空航天应用[11]方面具有巨大潜力。堆栈和热交换器是热声装置的关键组件，通常独立设计。在之前的研究中，提出了多种堆栈配置并进行了研究，也采用了不同的热交换器结构。平行板堆栈是最常见的配置，并已在大量研究中得到验证。Arnott等人[12]从理论上比较了几种孔隙几何形状，表明在短堆栈近似下，平行板堆栈能够产生有利的热量和功流动。Tijani等人[13]基于线性热声理论开发了一种热声制冷器的设计程序，并选择了平行板堆栈，因为其具有有利的热声性能和实际制造性。Mao等人[14]使用粒子图像测速（PIV）实验研究了平行板堆栈端部的振荡流动，并揭示了不同雷诺数下的不同涡流脱落模式。Aben等人[15]也使用PIV研究了平行板堆栈端部的涡流结构，并考察了板端形状和孔隙率对流动模式的影响。Alamir[16]实验研究了平行板堆栈的几何参数对驻波热声制冷器共振频率和性能的影响。Piccolo等人[17]构建了一个带有平行板堆栈的实验性驻波TAE，并在振荡流动条件下研究了具有圆形孔隙和矩形孔隙几何形状的热交换器的热性能。Piccolo和Jaworski[18]进一步比较了带鳍片的管式热交换器和圆形孔隙热交换器对驻波TAE性能的影响，发现圆形孔隙热交换器使发动机性能下降了约23%。除了平行板堆栈外，其他多种堆栈配置也引起了研究人员的关注。Lin等人[19]对具有轴对称环形堆栈的驻波热声-压电发动机进行了高保真全可压缩模拟，展示了堆栈几何形状对热声增长率和极限环行为的重要影响。Wantha[20]研究了不同堆栈配置（包括螺旋形、圆形孔隙和针阵列堆栈）对热声制冷系统性能的影响，并报告针阵列堆栈实现了更低的冷端温度。Chen等人[21]进行了三维数值模拟，研究了圆形孔隙和针阵列堆栈结构中水力半径与热穿透深度之间的最佳比例，发现在最佳几何参数下，针阵列配置的热声性能优于圆形孔隙配置。Wang和Hu[22]提出了一种使用印刷电路板热交换器（PCHE）和超临界CO2的驻波TAE，其中堆栈区域和热交换器均采用PCHE结构。他们的实验结果表明这种配置可以实现高压力幅度和低起始温度差。他们[23]进一步研究了PCHE通道尺寸对发动机性能的影响。

为了研究热声装置，广泛采用了数值模拟方法，并在之前的研究中提出了多种数值模型。Rott的线性热声理论[24]通过线性化基本守恒方程为分析热声现象奠定了重要理论基础。这一框架后来被Swift[25]、[26]扩展，使得在低幅度工作条件下可以预测和解释声学特性和热声装置的行为。Blanc-Benon等人[27]结合PIV测量和低马赫数模拟，可视化了堆栈边缘附近的非稳态速度和涡度场，并揭示了板厚对涡流结构的影响。Besnoin和Kino[28]使用基于涡度的低马赫数方案进行了非稳态数值模拟，研究了驱动比、板厚和热交换器几何形状对热声制冷器流动行为和性能的影响。Ma等人[29]、[30]使用混合低马赫数框架数值研究了驻波热声发动机，分析了阻抗负载对启动性能的影响，并进一步研究了活跃单元中的周期性流动动力学，强调了临界温度比、涡流脱落和热交换器锐角在热声发动机性能中的重要作用。Scalo等人[31]对行波热声热发动机进行了三维全可压缩模拟，展示了声流和非线性热声传输在极限环行为中的重要作用。Bouramdane等人[32]使用COMSOL Multiphysics进行了数值模拟，分析了堆栈参数、工作气体和边界条件对TAE性能的影响。Ja’fari等人[33]使用OpenFOAM数值研究了TAE中的声流。随后，Ja’fari和Jaworski[34]提出了“分布式流体泵”类比，以解释共振器沿线的非均匀平均压力分布。Baccoli等人[35]将基于DeltaEC的线性设计与基于OpenFOAM的非线性CFD模拟相结合，用于废热驱动的热声发电机，其中环形装置被简化为轴对称模型，再生器和热交换器使用多孔介质源项进行建模，包括热交换器的平行板多孔介质模型。Yao等人[36]应用基于MFEM框架的不连续Galerkin（DG）方法研究了四分之一波长TAE中的非线性效应，并进一步研究了不同堆栈参数对发动机性能的影响[37]。在上述研究中，两侧的热交换器在数值模型中被简化为与堆栈具有相同厚度的平板。Nowak等人[38]使用ANSYS CFX开发了一个简化的分析模型和数值模型来模拟TAE的操作。Chen等人[39]进一步使用ANSYS Fluent进行了大涡模拟，研究了四分之一波长驻波TAE。在这两项研究中，堆栈壁和热交换器都施加了简化的热边界条件，忽略了热交换器。Blanc等人[40]使用高保真软件PyFR模拟了高幅度驻波热声发动机中的非线性效应，在堆栈热边界施加了线性温度剖面。Slimene等人[41]使用格子Boltzmann方法研究了驻波TAE中的对流热传递和Rayleigh流动，其中热交换器和堆栈被简化为多孔介质。Yang和Hu[42]提出了一种CFD-ANN-DRL框架，用于优化超临界CO2 TAE中TAC的通道形状。为了简化数值模型，PCHE通道被近似为圆形截面通道，将模型简化为二维轴对称问题。

尽管对堆栈配置和热交换器设计进行了大量研究，但这些组件在传统热声系统中通常是独立优化的，导致其紧凑性有限。此外，数值研究常常采用简化的热交换器表示方法，这限制了准确捕捉TAC内部流动和热量传递过程的能力。为了解决这些限制，本研究提出了一种新颖的管束式TAC，将热交换器、冷交换器和堆栈集成到一个紧凑的结构中。这种集成设计简化了制造和组装过程，同时能够更真实地表示热交换器内部的流动和热量传递过程，而这些过程在以前的模型中常常被简化。考虑到管束结构的轴对称性，系统采用二维轴对称公式进行建模以降低成本。全可压缩Navier-Stokes方程使用基于MFEM框架的高保真DG方法求解[43]，并通过水平Rijke管的模拟对该模型进行了验证。系统研究了管半径对热声性能的影响，重点关注流动结构、热量传递和能量转换之间的耦合关系。

本文的其余部分安排如下：第2节介绍数值方法，包括控制方程和轴对称DG弱形式的推导。第3节通过水平Rijke管的模拟验证了数值模型。第4节对不同管半径的TAE性能进行了比较分析。最后，第5节总结了主要结论。

**2. 数值方法**
**2.1. 控制方程**
为模拟TAE，求解了全可压缩Navier-Stokes方程和傅里叶热传导方程。这些方程以保守矢量形式表示为：
(1)
???????? + ???????(??) ? ???????(??, ???) = 0，
其中保守变量??和原始变量??定义如下：
(2)
?? = ?? ? ? ??????3??????3???????? ? ? ??,
?? = ?? ? ? ???????????????? ? ? ??，
总通量????包括对流通量??????和粘性通量??????，其定义如下：
(3)
?????? = ?? ? ? ????3??????3????3???? + ??3????3??????3(??3?? + ??)??? ? ? ??,
?????? = ?? ? ? ??0????3??????3????3?? ? ????????????3????, ??3?????????????? ? ? ??，
其中总能量??、粘性应力张量????3??和粘性热通量????定义如下：
(4)
?? = ?????1 + 1/2???3????3????,
(5)
????3?? = ??3(??????3?????? + ????????????) ? 2/3??????????????3????3,
(6)
???? = ???3??????????，
其中??、??、??表示空间坐标方向。
动态粘度?被认为依赖于温度，使用Sutherland定律计算：
(7)
? = ?????3??3?? / (????3??3??3 + ?????? + ????)，
其中?????3??3??是参考温度????3??3??下的参考动态粘度，????是Sutherland常数。根据Prandtl数的定义，热导率从动态粘度计算得出：
(8)
? = ?????????????，
其中????是Prandtl数，????是恒压比热容。因此，在当前模型中，动态粘度和热导率都随温度变化。

为了在轴对称假设下简化求解过程，方程被重写为圆柱坐标（??,??,??）形式，其中????、????和????分别表示径向、方位角和轴向速度分量。为了简洁起见，这里不再展开圆柱坐标下的方程完整表达式。

**2.2. 轴对称DG弱形式**
为了使用DG方法求解控制方程，推导了轴对称弱形式[44]。圆柱坐标系中的计算域??和选定的轴对称半截面???3??如图1所示。任何点??∈??都可以用坐标（??,??,??）表示，其中（??,??）∈???3??且??∈[0,2π]。弱形式是通过将控制方程乘以一个测试函数并在域??2???上积分得到的。我们用????2???????????表示与??的外边界重合的??2???的边界部分，即????2???????????=????2???∩????。考虑了一个具有特征尺寸?的??2???的网格???，其中每个元素??∈???支持?∈??阶的多项式逼近。轴对称DG格式的标量试验/测试函数空间???定义为：(9)???={??∈??2?(??);???2???∈??2????;???(??,??,??)=??2????(??,??),?(??,??,??)∈??}，其中(10)??2????={??2???∈??2?(??2???);??2???∣??∈???(??),???∈???}。通过对控制方程施加??/????=0，得到了轴对称格式。还应该注意的是，在轴对称假设下，径向和方位角速度分量在??=0时必须满足????=????=0。对于在轴对称计算域??上定义的任何函数??，以下两个重要关系成立：(11)∫????d??=2????∫??2???????d??, (12)∫??????d??=2????∫????2???????d??，其中d???和d???分别表示??2???子午线截面中的面积和线元素。因此，可以对??上的轴对称函数进行积分，将其表示为在2D子午线截面??2???上的积分，对于????上的积分也有类似的关系。弱形式是通过将方程(1)乘以相应的测试函数并在每个元素??上积分来逐元素得到的。由于试验函数和测试函数都是轴对称的，因此可以将关系(11)、(12)应用于得到的表达式：(13)∫???????????????d??=∫??????????d???∫???????????????d???，其中数值通量??????用于确定元素界面之间的正常通量并确保相邻元素之间的信息交换。可以使用几种近似的Riemann求解器来实现这一目的，例如Godunov、局部Lax–Friedrichs和Roe通量。在这项工作中，采用了局部Lax–Friedrichs通量，垂直方向的数值通量表示为：(14)??????=???+??+2+??2????(??????+)，其中上标“?”和“＋”分别表示元素界面处的内部和外部状态。参数??代表与方向通量雅可比矩阵相关的局部最大特征速度，可以通过以下方式近似为[45]：(15)??=max??∈[???,??+]?(√?????????+√???????)。在DG方法中，可以通过减小网格尺寸、提高有限元阶数或同时应用这两种策略来获得数值收敛性，这一过程被称为hp收敛[45]。先前的研究表明，有限元阶数高于2即可提供足够的精确度，而进一步提高阶数虽然不会带来显著的改进，但会显著增加计算成本[36]。因此，在本工作中有限元阶数固定为2，并使用Gmsh[49]生成了四种不同元素大小的网格分辨率，以评估它们对TAE封闭端压力幅值的影响。为了确保数值稳定性，Courant–Friedrichs–Lewy（CFL）数保持在0.233，相应的时间步长约为3.0×10^-9秒。图7展示了在不同网格分辨率下封闭端（R=0.25毫米）的压力幅度。结果表明，当网格元素总数达到27121个时，压力幅度几乎保持不变，这表明已经实现了网格独立性。图8显示了管半径R=2.5×10^-4米时自激压力振荡的时间演化及其对应的FFT频谱。如图8(a)所示，压力振荡从初始的小扰动逐渐增长，最终达到近乎周期性的状态，表明稳定的自激热声振荡已经建立。这一过程大约需要0.25秒。图8(b)展示了压力信号的频域特性，主导频率约为595.04赫兹。图9展示了R=0.25毫米情况下HHX和CHX周围的网格分布。由于计算模型中两端的热交换器都采用等温壁面边界条件，因此没有明确建模它们的物理区域，因此没有生成相应的网格区域。由于HHX附近存在较大的温度梯度（尤其是在初始阶段），会产生强烈的压力梯度，需要更高的分辨率才能实现数值收敛，因此该区域的网格更为细腻。相比之下，CHX和其余区域的温度梯度较小，因此可以使用较粗糙的网格。在堆叠的上边界附近，由于不需要直接解析复杂的流动结构，因此使用了相对较粗的网格。对于热声核心外的流体区域，当温度和压力梯度较弱时，网格密度也会降低。这种网格策略在所有后续模拟中都得到了采用，无论管径如何。

为了进一步评估近壁分辨率，所采用的网格尺寸与粘性和热穿透深度进行了比较。粘性和热穿透深度分别计算为λv=√2μ/πη和λk=√2κ/πη，其中μ是运动粘度，κ是热扩散率，πη=2πf是角频率。在300–600 K的温度范围内，粘性和热穿透深度分别变化在9.16×10^-5米到1.66×10^-4米和1.09×10^-4米到1.97×10^-4米之间。在整个TAC区域内，径向网格间距保持在约1.67×10^-5米。由于本研究采用了二阶有限元近似并且验证了网格独立性，因此目前的网格分辨率被认为足以捕捉粘性和热边界层效应。因此，在所有后续模拟中采用了相同的近壁网格分辨率标准。由于计算域大小的变化，元素总数也随之变化。

下载：下载高分辨率图像（697KB）
下载：下载全尺寸图像

图9. 不同区域计算网格的示意图：(a) HHX周围；(b) CHX周围。黄色区域代表固体域，绿色区域代表流体域。

4.4. 结果与讨论
为了表征管径的影响，引入了管径与粘性穿透深度之比R/λv。因此，对于研究的管径范围2.5×10^-4米到5.0×10^-4米（间隔为5.0×10^-5米），相应的R/λv值为2.73、3.28、3.82、4.37、4.91和5.46。图10显示了不同R/λv下TAEs的最大压力幅度变化。图表表明，最大压力幅度最初随着R/λv的增加而增加。具体来说，当R/λv从2.73增加到3.28时，最大压力幅度增加了约66.38%。当R/λv进一步增加到3.82和4.37时，压力幅度继续增加，尽管增幅逐渐减小。当R/λv增加到4.91时，压力幅度与R/λv=4.37的情况相比变化不大，这表明最大压力幅度出现在R/λv介于4.37和4.91之间。继续将R/λv增加到5.46时，封闭端的压力幅度有所下降。为了进一步理解这一行为，以下部分从流动、传热和能量转换的角度分析了由R/λv表征的相对管径的影响。

下载：下载高分辨率图像（180KB）
下载：下载全尺寸图像

图10. 不同R/λv下TAEs封闭端最大压力幅度的比较。

4.4.1. 流动特性分析
图11展示了在不同R/λv下，圆形管两端径向速度分布的比较，相位角度分别为0°和180°。在一个振荡周期内，0°和180°的相位角度分别对应于流体速度向左右方向达到最大值的时刻。由于模拟中考虑了不同的管径，径向坐标使用特征长度进行归一化。归一化值0对应于管中心线，而值1代表管壁。在z=0.0241米处，向左流动（φ=0°）时的最大速度幅度高于向右流动（φ=180°）时的速度幅度；而在z=0.0459米处，观察到相反的趋势。此外，当流动方向对应于管的汇流侧时，速度最大值从管中心线向近壁区域移动。这种速度剖面的相位依赖性重新分布与振荡管流中的环形效应有关。在振荡运动过程中，轴向压力梯度的周期性反转导致惯性主导的核心流动和壁附近的粘性斯托克斯层产生不同的响应。核心流动受粘性影响较小，主要受流体惯性驱动；而近壁区域则受到粘性扩散的强烈影响。这两个区域之间的相位差可能在粘性边界层的外缘产生速度超调，导致近壁处的速度达到最大值，而壁处的速度由于无滑移条件而保持为零。这种现象被称为环形效应，Richardson和Tyler在振荡管流中通过实验观察到[50]，并且在具有近壁速度增强的热声振荡流中也有所报道[51]。流动面积的局部变化及其相关的涡旋结构可能会进一步改变速度分布，从而在热声核心的两端产生不同的相位依赖性剖面。尽管不同R/λv值的速度分布呈现相似的趋势，但将R/λv从2.73增加到3.28会导致两端峰值速度显著增加。当R/λv进一步增加到5.46时，速度分布的幅度逐渐减小。图11(b)展示了z=0.0459米处0°和180°相位角度的径向速度分布比较。远离管壁时，速度分布逐渐向右端移动；而在靠近壁的地方，它保持圆形管中Poiseuille流动的特征。在CHX侧观察到一个显著的差异，其中R/λv=2.73时的速度分布与其他R/λv时的差异更大。对于R/λv在3.82到5.46范围内的情况，两端在归一化径向范围0到0.2内的速度明显高于R/λv=3.28时的速度。相比之下，远离管壁时，所有R/λv的速度幅度排序与z=0.0241米时观察到的顺序一致。此外，值得注意的是，尽管R/λv=3.28时管内的速度达到最大值，但相应的压力幅度并不是最大的。

图12比较了不同R/λv下管中心线处的轴向速度分布，相位角度分别为0°和180°。尽管速度幅值随管径变化，但在两个相位角度下的总体分布模式保持一致。如图12(a)所示，在0°相位角度下，流体向左振荡，轴向速度达到最大值。当流动从右侧进入圆形管时，速度迅速增加，达到峰值后逐渐减小。经过HHX后，速度再次急剧增加。值得注意的是，对于不同的管径，中间区域的速度衰减速率大致相同。在180°相位角度下，流体向右振荡，轴向速度再次达到最大值，如图12(b)所示。对于所有管径，中心线速度最初迅速增加，并在HHX的右端附近达到最大值，随后逐渐减小。更远的右侧位置的速度幅度更大，速度衰减逐渐停止。随后，速度开始向CHX的右端逐渐增加。不同管径的速度增加速率几乎相同。

下载：下载高分辨率图像（554KB）
下载：下载全尺寸图像

图11. 不同R/λv下管两端径向方向的轴向速度剖面比较，相位角度分别为0°和180°，其中(a) z=0.0241米，(b) z=0.0459米。虚线表示φ=0°，实线表示φ=180°。负号表示流向左侧。

图13展示了整个TAC区域内不同R/λv下的时间平均速度幅度场。时间平均速度场是通过在一个振荡周期内对181个瞬时速度场进行平均得到的。对于所有研究的管径，管内的流动显示出具有相反流动方向的双层结构。近壁的平均流动向右移动，而核心流动向左移动，表明存在由振荡流动引起的声流。此外，在TAC两端突然的区域变化附近也可以观察到类似的双层流动结构。这些区域由于振荡流动与局部几何变化的相互作用而表现出复杂的流动模式，其中产生了涡旋结构，并进一步影响了流动分布。如图13所示，可以清楚地识别出声流强度随R/λv的变化。对于R/λv=2.73，近壁的平均流动相对较弱，主要局限于HHX侧，表明声流的发展受到限制。当R/λv增加到3.28时，近壁流动区域沿轴向显著扩展，表明流动强度增强。随着R/λv进一步增加到3.82，声流进一步增强，核心平均流动速度达到相对较高的水平。对于中等管径R/λv=4.37和4.91，流动仍然明显，尽管高速度平均流动区域相比R/λv=3.82时较为分散。然而，当R/λv增加到5.46时，总体时间平均速度幅度减小，表明较大管径处的声流减弱。

图14展示了不同R/λv下TAE的堆叠和CHX周围的涡度场比较，相位角度为0°。对于所有R/λv值，壁涡度层在右端截面扩展区域附近发生扭曲，并随后被抛入下游区域，表明流动分离是由流动面积的突然扩大引起的。分离的剪切层在管出口下游产生局部涡旋结构，导致形成回流区。这一回流区改变了有效的流动路径，并挤压了主要的振荡流动，从而影响了右端附近的轴向速度分布。这解释了图12(a)中右端观察到的速度突变现象。随着 ??/???? 从 2.73 增加到 3.28，膨胀区域附近的涡旋结构变得更加明显，并且向下游延伸得更远。当 ??/???? 进一步增加到 3.82 和 4.37 时，涡度分布依然强烈，并且在空间上发展得更加明显。对于更大的 ??/???? 值 4.91 和 5.46，下游的涡旋结构变得相对较弱且更加分散。这些结果表明，右端膨胀区域附近的流动分离和涡旋发展受到 ??/???? 的强烈影响。此外，强壁涡度层和分离的涡旋结构有助于解释图 11(a) 中在 ??=0° 时显示的轴向速度分布。最大轴向速度主要出现在 CHX 区域，其位置随着不同 ??/???? 下形成的涡旋结构而变化。

图 12. 不同 ??/???? 下沿管中心线的轴向速度分布，不同的相位角：(a) ??=0°，(b) ??=180°。负号表示流向左侧的方向。红色、灰色和蓝色区域分别代表 HHX、堆叠区和 CHX 区域。

图 15 比较了不同 ??/???? 下 TAE 堆叠区和 CHX 区周围的涡度场，在相位角为 180° 时的情况。在这个相位角下，振荡流从左向右。分离的涡旋结构主要集中在管出口和下游膨胀区域。与 ??=0° 的情况不同，??=180° 时的涡度分布不会在 CHX 区导致明显的局部阻塞。因此，图 12(b) 中的轴向速度在 CHX 附近没有突然变化，而是逐渐向管子的右端增加。随着管径的增加而管间距保持不变，管间间隙的有效流动面积变大。涡度场显示，管出口附近的分离涡旋结构逐渐向径向扩散到膨胀区域的上部和下部，而它们的轴向位移相对受限。

图 13. 比较了不同 ??/???? 下 TAC 区域内的时间平均速度场。

图 14. 比较了相位角为 0° 时不同 ??/???? 下堆叠区和 CHX 区周围的涡度场。

图 15. 比较了相位角为 180° 时不同 ??/???? 下堆叠区和 CHX 区周围的涡度场。

4.4.2. 热传递特性分析

图 16 展示了相位角为 0° 时 HHX 周围的温度等值线对比。在温度等值线中，红色区域代表高温区域，蓝色区域代表低温区域。白色区域对应于从高温区域到低温区域的过渡区域。此时，换热器左侧的高温流体开始从右向左快速振荡，而热端持续向流体供热，形成明显的瞬时高温积聚。热膨胀流体像活塞一样运动，导致压力变化。如图所示，管中心的温度分布相对均匀，而靠近壁面的等温线更密集。这表明在近壁区域，由于热粘性边界层的作用，流体速度降低，热扩散占主导。相比之下，管中心的对流传输导致温度分布更加均匀。随着管径的变化，温度分布表现出不同的特点。对于 ??/????=2.73，较小的管径导致较大的热边界层比例，高温区域集中在换热器内部左侧。在这种情况下，强烈的粘性效应对流体振荡产生显著阻力。随着 ??/???? 增加到 3.28，壁温梯度和声流强度增强，而粘性效应相对减弱。因此，高温区域被推出 HHX 管并进一步延伸到下游区域。随着 ??/???? 进一步增加到 3.82–4.91，高温区域的径向分布更加均匀，等值线显示高温区域逐渐扩展。然而，当 ??/???? 达到 5.46 时，高温区域开始向管内的热端收缩。这归因于局部壁温梯度的进一步减小，从而减弱了声流强度及其对中心高温流体的影响。

图 17 展示了相位角为 180° 时 HHX 周围的温度等值线对比。此时，流体从 HHX 处远离振荡。对于较小的 ??/???? 值 2.73–3.28，高温区域主要限制在 HHX 管内部，向右端的扩散有限，表明壁面粘性效应对流体运动产生强烈阻力。随着 ??/???? 的增加，高温区域向右端更加明显地扩散，反映出壁面粘性影响的逐渐减弱。相反，对于较大的 ??/???? 值 4.91–5.46，中心区域的高温流体向右侧的影响减小，表明声流强度的减弱。通过结合相位角为 0° 和 180° 的温度等值线可以推断，当 ??/???? 处于 4.37–4.91 范围内时，声流和壁面粘性效应之间达到了相对最佳的平衡。这种平衡使热流体在管内进行更大的往复振荡，从而产生更高的 TAE 压力幅度。这也解释了为什么尽管在 ??/????=3.28 时速度幅度达到最大，但相应的压力幅度并不是最大的。

图 18 比较了 TAE 达到稳态时不同 ??/???? 值下堆叠壁温的轴向分布。温度升高主要集中在堆叠区域的左侧。随着 ??/???? 从 2.73 增加到 3.28，整体壁温相应上升。进一步将 ??/???? 增加到 3.82 时，温度分布只有轻微的变化。对于 ??/???? 值 4.37 和 4.91，则观察到更明显的温度升高，而在 ??/????=5.46 时，壁温下降，接近 ??/????=3.82 时的水平。温度等值线的分析表明，堆叠壁温分布与声流输送高温流体的有效性密切相关。更强的声流增强了堆叠内热流体的轴向传输，导致更高的壁温以及系统中储存的热能增加。

图 19 比较了不同 ??/???? 下 TAC 区域沿管壁的时间平均热流密度分布。所有情况都显示相似的整体分布。HHX 区域的时间平均壁热流密度为正，并逐渐向 ??≈0.0275 m 处的 HHX–堆叠接口增加。在该位置，振荡流体在一个声学周期内经历最大的温度变化，导致最大的壁热流密度。进入堆叠区域后，时间平均壁热流密度变为负值。这表明，在往复振荡过程中，热工作流体将热能传递给堆叠壁，部分能量暂时储存在固体材料中。在该区域下游，负热流密度的幅度迅速减小，在大约 ??=0.033 m 处接近零。因此，堆叠区域的后半部分对热能储存的贡献有限。

图 18. 比较了稳态下不同 ??/???? 下沿堆叠壁的轴向温度分布。

图 19. 比较了不同 ??/???? 下 TAC 区域沿管壁的时间平均热流密度分布。

4.4.3. 能量转换分析

为了说明 TAE 内局部工作的空间分布，我们在之前的工作中记录了堆叠周围压力和比体积的时间变化 [37]。一个声学周期内压力-比体积曲线所围成的面积定义为局部能量密度 ??loc，作为 TAE 内不同空间位置的局部工作性能的衡量标准：(19)??loc?(??,??)=∫???(??,??,??)?d?[1???(??,??,??)]。图 20 展示了不同 ??/???? 值下 ??=0.02–0.0241 m 区内的局部能量密度等值线分布。该区域对应于圆管的左侧部分，在 ??=0.0241 m 处流截面面积发生突变。突然的面积变化引起了局部涡旋结构，导致额外的能量耗散。当 ??/????=2.73 时，等值线最右侧部分的局部能量密度为负，而在其余区域则大致为零。这表明该区域没有产生有效的工作，反而发生了局部能量耗散。随着 ??/???? 增加到 3.28，流动截面面积扩大，局部能量密度为负的区域向外移动。同时，管的中心区域逐渐较少受到涡旋结构的影响。等值线中心附近开始出现局部能量密度为正的区域，表明局部产生了有效工作。随着 ??/???? 进一步增加到 3.82，局部能量密度为正的区域扩大，正负局部能量密度的极端值都增加，表明局部能量交换的强度增强。然而，随着 ??/???? 继续增加，局部能量密度为正的区域开始缩小，其极端值逐渐减小。这意味着过大的流动面积削弱了局部能量转换强度，从而降低了有效工作输出。

图 21 展示了不同管径下 ??=0.0241?0.035 m 区内的局部能量密度等值线分布。对于所有管径，正能量密度区域大致位于相同的位置，并显示出明显的菱形图案。这个菱形区域的中心主要位于 HHX 和堆叠之间的界面附近，那里存在相对较大的温度梯度。这一特征主要归因于局部温度梯度和管壁附近的粘性边界层的共同作用。边界层内的强粘性耗散减弱了靠近壁面的能量转换能力，导致远离壁面的正能量密度区域更加明显。在该区域之外，大多数其他区域的能量密度接近零，表明有效工作生成主要集中在温度梯度较大的区域。当 ??/????=2.73 时，流体流动受到靠近壁面的粘性边界层的强烈影响。因此，峰值正能量密度相对较小，相应的区域范围也有限。将 ??/???? 增加到 3.28 会导致峰值正能量密度显著上升，正能量密度区域明显扩大。随着 ??/???? 进一步增加到 3.82，峰值继续增长，而区域的横向范围仅略有变化。然而，随着半径的进一步增加，峰值正能量密度逐渐减小。同时，远离壁面的区域中声流的驱动效应减弱，导致管中心附近的局部能量密度减小，正能量密度区域收缩。

图 20. 比较了不同 ??/???? 下 ??=0.02?0.0241 m 区内的局部能量密度场。

图 22 和图 23 分别展示了 TAC 下游部分和 CHX 下游膨胀区域内的局部能量密度分布。由于这些区域的温度梯度较小，局部能量密度的幅度通常较低。尽管如此，由于粘性边界层内的粘性耗散，管壁附近仍然可见负的局部能量密度。同样，CHX右端流面积突然扩大所产生的涡旋结构也会导致局部区域出现负能量密度。不同情况之间的比较表明，随着??/????的增加和流速的增大，负局部能量密度的幅度也在增加。根据方程（19）定义的局部能量密度，整个轴对称计算域内的正、负和净循环功可以计算为：(20)??+cyc=2????∫??????????max?(??loc,0)???d???d???, (21)???cyc=2????∫??????????min?(??loc,0)???d???d???, (22)??netcyc=??+cyc+???cyc，其中???????表示循环平均密度：(23)????????(??,??)=1??cyc?∫??cyc0???(??,??,??)????，这里??cyc是声学周期。从一个声学周期内从HHX壁传递到工作流体的热量可以计算为：(24)??cyc=∫??cyc0∫??HHX????d??d??，其中????>0表示从HHX壁到工作流体的热传递。对于轴对称模型，上述方程可以写为：(25)??cyc=2????∫??cyc0∫??HHX?????(??,??,??)???d??d????d???。然后定义循环功转换比为：(26)??cyc=??netcyc??cyc×100%。需要强调的是，??cyc表示整个声学周期内域积分的净循环功与从HHX传递的热量之间的比率。由于当前模型中没有包括外部声学负载，这个量不应被解释为TAE的整体热效率或可提取的声学功率效率。表3总结了不同??/????的正、负和净循环功，以及每个周期的热传递量和相应的循环功转换比。随着??/????从2.73增加到4.91，正循环功显著增加。这种增加部分归因于局部热声功生成的增加，部分归因于随着管半径的增加而包含在域积分中的流体体积的增大。与此同时，负循环功的幅度也增加了，反映了与粘性和涡旋耗散相关的更强局部声能消耗，以及积分流体体积的增加。由于当前的轴对称模型仅代表了TAC内的一个小特征流动通道，因此循环功和每个周期的热传递量的绝对值相对较小。此外，相对较低的循环功转换比可能与由粘性源区域引入的较强阻尼有关，这增加了系统中的声能耗散。然而，对所有调查的案例都应用了相同的粘性源处理方法。因此，虽然绝对值不应被解释为完整热声发动机的整体性能，但计算出的循环功转换比仍然适用于比较不同??/????之间的相对能量转换性能。净循环功在??/????=4.91时达到最大值。然而，最高的循环功转换比（约1.89%）出现在??/????=4.37时。当??/????从4.37增加到4.91时，净循环功的增加相对较小，而每个周期的热传递量继续增加，因此转换比略有下降。当??/????进一步增加到5.46时，净循环功明显减少，而热输入保持较高，导致循环功转换比进一步降低。这些结果表明，最佳??/????可能位于4.37和4.91之间，在这里热声功生成、局部耗散和热输入之间达到了一个有利的平衡。

下载：下载高分辨率图像（480KB）
下载：下载全尺寸图像
图21. 不同??/????下??=0.0241?0.035 m范围内的局部能量密度场比较。
下载：下载高分辨率图像（539KB）
下载：下载全尺寸图像
图22. 不同??/????下??=0.035?0.0459 m范围内的局部能量密度场比较。
下载：下载高分辨率图像（513KB）
下载：下载全尺寸图像
图23. 不同??/????下??=0.0459?0.05 m范围内的局部能量密度场比较。
表3. 不同??/????的域积分循环功特性和转换比。
??/???? ??+cyc (J) ???cyc (J) ??netcyc (J) ??cyc (J) ??cyc (%)
2.73 5.35×1?0 ?7 ?5.9 ×1?0 ?8
4.76 4.82×1?0 ?7 4.8 ?5
0.99 3.28 1.81×1?0 ?6 ?3.0
1.51 1.58 ×1?0 ?6
9.52 3.82 2.55×1?0 ?6 ?5.5
1.99 1.78 ×1?0 ?6 1.12
4.37 2.88 1.89 ×1?0 ?6 ?5.9
2.28 2.28 ×1?0 ?6 1.20
4.91 2.92 1.83 ×1?0 ?6 ?6.16
2.30 1.83 ×1?0 ?6 1.26
5.46 2.63 ×1?0 ?6 ?6.29
2.00 1.27 ×1?0 ?6 1.58

5. 结论
本研究提出了一种新颖的管束热声核心（TAC）结构，用于热声发动机（TAE），具有紧凑和集成的配置。该核心由管束组成，管束之间由管板分隔，形成了包括热交换器（HHX）、冷交换器（CHX）和堆栈区域的统一结构。为了降低计算成本，问题被简化为二维轴对称形式，并开发了基于MFEM框架的轴对称不连续Galerkin（DG）数值模型。通过模拟水平Rijke管来验证数值模型的准确性。使用经过验证的模型，系统地研究了管半径对TAE性能的影响，分析了流动、热传递和能量转换特性。主要结论如下：
(1) 压力幅度受到??/????的强烈影响，其最大值出现在??/????=4.37–4.91，对应于管半径为??=0.4–0.45 mm。
(2) 不同??/????的瞬时轴向速度剖面具有相似的空间分布，但其峰值幅度出现在??/????=3.28，而不是在压力最优范围内。时间平均流动包括向右的近壁层和向左的核心层。声流在??/????=3.82之前增强，在??/????=4.37–4.91时保持相对较强，在??/????=5.46时减弱。下游扩管处的分离产生了涡旋结构，其强度和空间范围随??/????变化，并改变了局部轴向速度。
(3) 在较小的??/????下，高温区域局限于HHX管内；随着??/????增加到4.91，它向下游扩展得更远，堆栈壁温度升高。在??/????=5.46时，高温区域由于壁温梯度和声流减弱而收缩。对于所有情况，HHX中的时间平均壁热流是正的，并在HHX–堆栈接口附近达到最大值，然后在上游堆栈中变为负值，并在下游趋近于零，表明热能主要储存在堆栈的前半部分。
(4) 正的局部能量密度主要集中在HHX–堆栈接口附近，而负值出现在管壁附近和下游扩管区域，这是由于粘性和涡旋耗散。正循环功和净循环功增加到??/????=4.91，同时负循环功的幅度也在增加。最高的循环功转换比在??/????=4.37时获得。在??/????=5.46时，尽管热输入较高，净循环功和转换比都降低了。
这些结果验证了管束TAC的可行性，并表明通过平衡声流、堆栈中的热存储、热声功生成以及粘性和涡旋耗散，可以实现有利的操作范围。这些结果为设计紧凑高效的热声发动机提供了指导。

CRediT作者贡献声明
陈尧：撰写——原始草稿、软件、方法论、概念化。
李博涵：方法论、概念化。
高博：撰写——审阅与编辑。
刘继平：撰写——审阅与编辑、资金获取、概念化。
严俊杰：监督、概念化。

中国西安交通大学热流科学工程国家重点实验室，西安710049