双开口气波制冷机非定常流动三维数值模拟

2022-08-25吴响响王泽武许晓飞胡大鹏

流体机械 2022年7期

吴响响，王泽武，许晓飞，胡大鹏

（1.中铁第四勘察设计院集团有限公司，武汉 430061；2.大连理工大学化工机械与安全系，辽宁大连 116024）

0 引言

双开口气波制冷机（DOGWR）具有结构简单紧凑、参数可调节和制冷效率较高等优势，被广泛应用于航空发动机试验、天然气脱水和工业尾气回收等领域［1-2］。目前国内外学者对于泵等水力机械的非定常流动特性进行了大量研究［3-7］，但主要集中在不可压缩的液体，而对于非定常瞬态气体这类可压缩、具有能量转化和波系产生的流体研究较少。

王静娴等［8］建立了双开口气波制冷机实验平台，监测了气波机振荡管内气体脉动压力，对气波机内部波系和流动机理进行了研究。于洋等［9］搭建了气波冷凝装置，对气体热力学特性进行了详细研究，获得了影响气波机制冷性能和脱湿干气压力恢复性能的主要参数。代玉强［10］和刘培启等［11］通过实验研究了气波机振荡管制冷性能，对影响气波机制冷效率的膨胀比、转速和高温端口压力等参数进行了测试。这些实验研究在一定程度上为气波机的结构优化和操作提供了指导，但现有实验条件只能研究振荡管内压力振荡特性，而对于管内流体温度、速度等其他重要物理量无法捕捉，且实验难以直观体现振荡管内气体的非定常流动。

随着近些年数值计算的迅速发展，数值模拟可以直观体现复杂流场、多物理场特性，从而弥补传统实验的不足，因此数值计算作为一种科学的研究方法也被广泛应用于气波制冷机中［12-14］。SPALDING等［15］利用一维非定常流动理论，忽略了波转子的周向和径向，对振荡管进行了非定常流动分析，并在DANIEL等［16］研究下进一步实现了程序化。于洋等［9］将气波制冷机沿周向展开得到平面移动模型，采用二维周期性边界条件对气波机进行了流场、温度场的数值模拟，并优化了气波机冷凝波图。刘培启等［17-18］通过数值模拟和实验定量分析了高温口压力对气波机制冷和内部流场的影响，并对气波机最佳偏角设置进行了预测。然而目前使用的气波机一维、二维数值模拟尚未考虑波转子的旋转，忽略了高压气体的入射损失，因此不能真实反应管内流体边界层效应和漩涡等复杂三维现象，也不能真实模拟高压气体射入振荡管这一关键渐开渐闭过程。同时，一维、二维数值模拟获取的激励载荷也无法加载到振荡管上对气波机动态力学响应等问题进行科学研究。

为此，论文首先分析双开口气波制冷机的结构和原理，基于周期性条件开发气波机三维非定常流动数值模型，并对数值结果进行验证。然后分析振荡管内非定常流动流体的压力场和温度场等特性，重点关注高压喷嘴气体射入振荡管这一关键过程。最后提取振荡管内多个关键点的激励载荷，从而为进一步研究气波机压力交换制冷机理、振荡管流致振动问题和延寿提供技术支持。

1 结构和原理

双开口气波制冷机是一种利用激波、膨胀波和压缩波之间波系运动和能量交换实现气体制冷的特殊设备［8，11］，其常见结构如图 1 所示［9］。

图1 双开口气波制冷机结构Fig.1 Structure of the DOGWR

气波制冷机的核心部件波转子由内、外圆筒和插片组成，从而在内、外筒之间的整个圆周上形成多个通道（振荡管）。同时，带有螺栓孔的波转子与旋转轴之间通过螺栓连接，进而实现波转子随转轴同步转动。波转子两端分别设置有静止的高压进气端口、高温出气端口、低温出气端口和中压进气端口，且高压进气端口与高温出气端口中心存在一定的偏角，偏角和气体进出端口的大小以及波转子的旋转速度可以调节，从而改变气波机的工况，实现气波机的最高制冷效率。另外，进、出气端口固定端板与波转子两端面保持有一定的间隙，防止波转子运转时与端板接触挤压而发生卡死现象。

双开口气波制冷机正常运转后，高压喷嘴持续向波转子通道中射入高压气体，并对管内原有气体进行压缩，使得原有气体温度升高，射入的新鲜高压气体膨胀降压，温度降低。由于振荡管周期性的与进、出气端口逐渐接通和闭合，因此振荡管内产生非定常流动流体。本文双开口气波机的工作介质为空气，其数值模型尺寸与文献［8］实验装置结构尺寸参数相同，具体参数见表1。

表1 双开口气波制冷机结构参数Tab.1 Structural parameters of the DOGWR

2 三维数值模型开发

2.1 模型简化

双开口气波制冷机三维流动过程数值计算对计算机性能要求极高，且旋转振荡管与进、出气端口渐开渐闭，因此，数值模拟的主要难点是合理简化模型，使现有计算机能够满足计算要求。开发的双开口气波制冷机三维模型如图2所示，主要考虑以下3个方面：

图2 双开口气波制冷机简化三维模型Fig.2 Simplified 3D model of the DOGWR

（1）振荡管是气波机关键也是薄弱零件，因此三维模型主要考虑振荡管、转毂和进、出气端口，而忽略气波机旋转轴、接管和外壳等结构。同时，以部分振荡管代替整机，并将进、出气端口轴向长度缩小，以节约计算资源。

（2）3根振荡管和管内流体模型以及滑移边界条件既可以真实模拟高压端口气体入射，又便于脉动载荷传递和考虑振荡管之间的相互影响。同时，带有螺栓孔的转毂模型有利于准确设置振荡管在高频激励载荷下动态响应的边界条件，从而最大程度的模拟气波机实际运行工况。

（3）中间振荡管初始位置与气波机的高压进气端口十分靠近，便于在较短的计算时间内获取振荡管与高压端口渐开渐闭这一关键过程的管内脉动流体流场，从而节约计算资源。

2.2 网格划分

气波制冷机流体域和固体域网格一起划分会导致三维模型网格数量倍增，且算法复杂。因此，为了提高计算效率，进一步简化数值计算模型，只对气波机三维流体模型进行网格划分和流场计算分析。为了更好地捕捉振荡管内脉动流体的流场，尤其是管内产生的各种激波和膨胀波，采用“多区域网格划分”方法，使用“六面体”和“六面体主导”的网格类型。为了消除网格大小对计算结果的影响，进行了多次细化网格试算，最终确定了整体结构网格大小为1 mm×1 mm，同时为了准确模拟流固交界处的流体流动，对振荡管内壁接触的流体表面设置6层膨胀层和一定的生长率，使得边界层网格向流体内部区域逐渐变大，从而与流体内部区域网格完好衔接，便于流体表面区域与内部区域间进行信息传递和计算。最后得到结构化流体网格模型的单元数量和单元节点分别为644 400，707 323个，气波机三维流动数值计算网格模型如图3所示。

图3 气波机三维流动数值计算网格模型Fig.3 3D finite element mesh model of the DOGWR

2.3 边界条件和算法设置

（1）边界条件。

双开口气波制冷机三维流体流动模型边界条件主要有5种：高压进气端口（HP inlet）和中压进气端口（MP inlet）的“压力进口”边界（Pressure inlet）；高温出气端口（HT outlet）和低温出气端口（LT outlet）的“压力出口”边界（Pressure outlet）；振荡管内流体左端面和高压喷嘴、低温气腔的接触面设置为一对“滑移边界”（Sliding interfaces），同时振荡管内流体右端面和高压喷嘴、中压气腔的接触面设置为另一对“滑移边界”；此外，振荡管、喷嘴和气腔内与设备接触的流体面均设置为“壁面”边界（Wall）；其他内部流体软件默认为“内部”边界（Interior）。气波机数值计算进、出端口边界条件见表2。

表2 进、出气端口边界条件Tab.2 Boundary conditions of the inlet and outlet ports

另外，振荡管内流体非定常流动的数值计算首先需要对计算流体域赋予初值，使流场特性较快到达气波机正常运转时工况。为保证计算的收敛性，喷嘴和气腔的初始化条件按照表2中气波机对应端口的边界条件给定，振荡管内气体初始压力和温度分别设置为0.5 kPa（表压）、299.8 K。

（2）算法设置。

基于FLUENT软件中的瞬态分析、密度基和滑移网格计算振荡管内脉动流体的流场特性。其中，密度基求解器适用于激波干扰、伴有燃烧等的高速可压缩流体和超音速流动流体［19-20］。振荡管内流体采用理想气体作为计算介质，其在低压、高温下与实际工作介质空气的偏离在工程允许误差范围内，已广泛应用于气波机数值模拟［8-9，17-18］。计算模型中近壁面采用“标准壁面函数”，具体控制方程采用Realizable k-ε湍流模型，该模型对于带有旋流的平面混合流、平面射流和圆形射流管内充分发展的湍流流动数值计算更加准确［17，21］，适合高压喷嘴射流工况数值计算。采用AUSM离散格式，其对间断处较为敏感，具有捕捉激波强的优势，已广泛应用于可压缩气体动力学领域［22］，可以准确获取管内多种波系和强间断。模型计算Courant数缩小设置为3，可以使计算收敛更好。同时，对振荡管内流体域设置2 340 r/min的旋转速度模拟管内流体随振荡管的转动。转动振荡管与固定进、出气端口间采用滑移网格技术，滑移网格边界采用Interface边界条件，最后得到双开口气波制冷机三维流体计算模型如图4所示。

图4 双开口气波制冷机三维流体计算模型Fig.4 3D fluid calculation model of the DOGWR

三维流体流动数值计算模型对应双开口气波制冷机正常运转时的较佳制冷转速是2 340 r/min，因此振荡管每旋转1°大约需要7.122 5×10-5s。为了更好的捕捉振荡管内非定常流动流体的激波、膨胀波、压缩波和接触面等强间断，设置计算时间步长为1×10-5s，计算稳定性较好。同时，设置数值计算的总时间步数为2 550步，则振荡管大约可以旋转整个周期，计算过程中所有物理量的残差基本小于1×10-3，计算收敛性较好。

3 结果验证与分析

3.1 数值模型验证

为了验证开发的双开口气波制冷机三维流动模型的准确性，利用现有文献实验数据与数值解进行对比分析。文献［8］在波转子的某一根振荡管上开设了6个测压孔，并在对应位置安装了高频动态压力传感器，从而对管内流体的脉动压力进行监测。

在气波机三维流动数值计算完成后，选取如5所示6个压力监测点，从而提取对应开孔处流体的压力实时数据。最后得到如图6所示的振荡管监测孔3，4的压力变化曲线。

图5 振荡管内流体压力监测孔Fig.5 Position of the pressure monitoring hole in the oscillation tube

从图6中可以看出，总体趋势上试验值与计算值基本吻合，也存在一定的误差，其产生的原因大致包括：数值计算中采用理想气体和Realizable k-ε湍流模型，与管内真实气体的物理特性和实际流动有一定差别；文献中试验气波机的高压端口入射气体到振荡管内有涡流产生，造成试验中压力传感器测量值有一定波动。虽然振荡管内脉动压力计算结果和试验测量值有一定的误差，但从整体上看，数值计算压力脉动值与试验测量值吻合度较高，计算模型可以较好地捕捉激波、压缩波和膨胀波的位置，因此本文开发的双开口气波制冷机三维流体流动数值模型可以用于研究气波机核心部件振荡管内的流场特性、入射损失、高频脉动载荷提取和瞬态热流固耦合力学响应等科学难题。

图6 振荡管监测孔3和4试验与数值压力时域对比Fig.6 Comparison curves of the pressure at openings 3 and 4 between experimental and numerical values

3.2 振荡管内流场分析

（1）压力场。

为了便于研究气波机振荡管内流体压力随振荡管旋转的波动特性，充分利用波转子振荡管周期性与进、出气端口接通和闭合，同时考虑振荡管之间的相互影响，选用3根振荡管模型的中间振荡管内流体为研究对象，提取其多个不同时刻的压力云图，并拼接成一个整体进行分析，得到中间振荡管与喷嘴接通和闭合过程中管内非定常流动流体的压力分布云图如图7所示。

图7 中间振荡管内脉动流体不同时刻静压分布云图Fig.7 Static pressure distribution nephogram of the fluid in the middle oscillating tube at different moments

从图7（a）中可以看出，当高压进气端口与波转子通道逐渐接通时，高压气体射入振荡管通道内，同时产生了激波，随着通道的旋转，激波向通道右侧快速运动，激波所到之处通道内气体压力上升。当中间振荡管旋转到与高压进气端口闭合时会产生膨胀波，膨胀波作用下气体压力降低。此外，当振荡管与高压、高温端口同时接通时，管内流体的压力普遍较大，同时波动剧烈；而当振荡管旋转至远离高压、高温端口时，虽然管内流体仍存在一定波动，但是波动幅度和压力都相对较小。从图7（b）中可以看出，当振荡管运动到低温排气阶段时，振荡管内流体压力普遍较小且相对稳定，压力在0.321～0.357 kPa之间波动，因此，可以将低温排气阶段简化为稳态进行数值分析。

为了便于研究振荡管与高压端口渐开渐闭复杂流动过程，单独取中间振荡管与高压喷嘴渐开、完全接通和渐闭3个时刻，并沿管内流体不同高度和不同轴向位置设置多个平面，得到振荡管与高压喷嘴不同时刻管内不同位置流体压力云图如图8所示。从图8（a）中可以看出，当振荡管与高压喷嘴部分接通时，射入管内的高压气体压缩管内原有空气，使得靠近高压端口的振荡管内流体压力快速升高，而远离高压端口的管内流体由于射入的高压气体还没有到达，因此压力变化较小。同时，由压力放大图可知，振荡管与高压喷嘴渐开时靠近高压端口的管内流体沿通道宽度（周向）方向有一定压力梯度，这是由于振荡管与高压喷嘴渐开导致射入管内高压气体在管内分布不均匀，但是远离高压端口的管内流体沿通道宽度方向的压力分布基本相同。从图8（b）中可以看出，当振荡管与高压喷嘴完全接通时，振荡管内流体压力较大，且沿着振荡管长度方向的压力出现明显波动。同时，振荡管与高压喷嘴完全接通时管内流体沿通道宽度方向的压力分布基本相同。从图8（c）中可以看出，当振荡管与高压喷嘴渐闭时，振荡管内流体压力较大，最大压力达到83.1 kPa，发生在振荡管中间。同时，振荡管与高压喷嘴渐闭时管内流体沿通道宽度方向的压力分布也基本相同。此外，从图8中可以看出，振荡管与高压端口渐开渐闭整个过程中管内不同高度截面的流体压力分布云图基本相同。

图8 不同工况下振荡管内流体压力云图Fig.8 Pressure nephogram of the fluid under different working conditions

（2）温度场。

中间振荡与喷嘴接通和闭合过程中管内非定常流动流体的温度分布如图9所示。

图9 中间振荡管内脉动流体不同时刻静温分布云图Fig.9 Static temperature distribution nephogram of the fluid in the middle oscillating tube at different moments

由图9可知，当高压进气端口与波转子通道逐渐接通时，高压气体射入振荡管通道内并产生了激波，随着通道的旋转，激波向通道右侧快速运动，激波所到之处气体温度上升形成一条高温带，最高温度达到372 K。当中间振荡管旋转到与高压进气端口闭合时会产生膨胀波，膨胀波作用下气体温度降低。同时，随着振荡管旋转产生的离心力以及管与喷嘴之间的渐开渐闭过程，管内流体温度在冷热交界面发生了扭曲，说明冷热气体在接触面强间断附近发生了混合，这也使接触面有一定宽度。此外，当振荡管与高压、高温端口同时接通时，管内流体整体温度较大，且靠近高温端口的管内温度最高，这与振荡管内高温气体从高温端口排出的设计原理一致。而当振荡管旋转至远离高压、高温端口时，管内流体温度都相对较小，且靠近低温端口的管内流体温度最低，实现了气波制冷机中低温气体从低温端口排出。

振荡管与高压喷嘴渐开、完全接通和渐闭时刻管内不同位置流体的温度分布云图如图10所示。

图10 不同工况下振荡管内流体温度云图Fig.10 Temperature nephogram of the fluid under different working conditions

从图10（a）中可以看出，当振荡管与高压喷嘴部分接通时，振荡管内靠近高压端口部分的流体温度急剧上升，而远离高压端口的流体温度变化相对较小。同时，由温度放大图可知，振荡管与高压喷嘴渐开时靠近高压端口的管内流体沿通道宽度方向有很大的温度梯度，这是由于射入管内的高压气体分布不均匀导致，且不均匀气体沿通道长度方向产生了漩涡，发生了冷热气体掺混，但是远离高压端口的管内流体沿通道宽度方向的温度分布基本相同。从图10（b）中可以看出，当振荡管与高压喷嘴完全接通时，振荡管内流体温度较大，且沿着振荡管长度方向的温度出现明显不连续，同时管内流体冷热温度交界面发生很大扭曲变形。另外，振荡管与高压喷嘴完全接通时管内流体沿通道宽度方向的温度分布基本相同。从图10（c）中可以看出，当振荡管与高压喷嘴渐闭时，振荡管内流体较高温度均分布在靠近高温出气端口的振荡管内，同时，振荡管与高压喷嘴渐闭时管内流体沿通道宽度方向有一定温度梯度。此外，从图11中可以看出，振荡管与高压端口渐开渐闭整个过程中振荡管内高度方向上的流体温度分布基本相同，与压力分布趋势一致。

3.3 激励载荷提取

研究高频激励载荷的精确获取方法，掌握振荡管内气体波动诱发的激励载荷特性，为气波制冷机结构设计优化和减振降噪等科学问题提供研究基础。基于本文开发的气波机三维流体数值计算结果，提取得到振荡管监测点No.1～No.6对应的压力、温度变化曲线如图11，12所示。从图11可以看出，振荡管内监测点No.1～No.6位置的流体压力在-15～85 kPa之间剧烈波动，其中压力上升最快的是No.1，最慢的是No.6。一个旋转周期内振荡管中同一区域气体压力的波动峰值随时间增加逐渐减小，且No.1～No.6中气体静压最后基本稳定到0附近微弱波动，管内流体压力恢复到初始状态，完成一个周期的载荷变化。从图12可以看出，振荡管内监测点No.1～No.6位置的流体温度在253～360 K之间剧烈波动，其中温度上升最快的是No.1，最慢的是No.6。No.1～No.6的温度时域图中温度最高的是No.5，最小的是No.1。No.1～No.6中最先稳定并恢复常温的是No.6，最后恢复常温的是No.1，且振荡管旋转一周后，管内流体温度逐渐恢复初始状态。