波瓣混合器的大涡模拟

许锦锦，边一帆，赵马杰，叶桃红

（中国科学技术大学热科学和能源工程系，安徽合肥230027）

波瓣混合器是一种增强同向混合流动的装置，在航空发动机的排气系统中，中心射流与次流的混合效率对航空发动机的喷射噪声有重要影响，因而波瓣混合器最初用于降低排气噪声（见图1）[1-4]，随后研究发现波瓣混合器能够增加发动机推力[5]、减少燃油消耗[6-7]，降低红外辐射[8]，此外波瓣混合器也应用于增强燃料与空气在燃烧室的混合作用，用以提高燃烧效率和减少污染物的形成[9]，因而波瓣混合器在航空发动机中有着广泛的应用。所以开展对波瓣混合器增强混合机制的研究具有重要意义。

图1 波瓣混合器

许多研究人员利用实验对波瓣混合器的流场结构进行了捕捉，用以分析其增强混合作用的机制。Paterson[10]最早利用实验的方法对波瓣混合器的混合机制进行了研究，研究发现波瓣混合器内存在大尺度的流向涡、正交涡以及小尺度的涡流结构，但对其中的混合过程并不十分清楚。Werle[11]等通过对二维波瓣结构的实验测量分析，揭示了流向涡的形成过程，图2所示为波瓣混合器内存在的几种涡流结构。McCormick[12]等通过实验的方法指出由于K-H不稳定性形成的正交涡对加强混合具有同等重要的作用，同时发现正交涡和流向涡的相互作用产生了更小尺度的涡流结构。Glauser[13]等利用热线测速的方法对二维波瓣混合器进行了研究，研究发现流向涡导致了正交涡的变形，最终致使其破碎，从而形成了小尺度湍流结构。Hu[14]751等利用实验分析了正交涡和流向涡随着流向距离的衰减情况，发现它们具有同样的衰减趋势，进一步证明了正交涡与流向涡之间的相互作用。Belovich[15]总结了波瓣混合器的混合过程主要受到三个因素控制：①由于波瓣混合器的几何外形而诱导产生的流向涡；②自由射流的剪切层中K-H不稳定性而形成的正交涡；③由于波瓣混合器特殊的几何外形结构增大了射流与次流的接触面积，加强了剪切作用。但是，实验方法对波瓣混合器流场的测量受到诸多因素的限制，对波瓣混合器的研究具有局限性，因而开展对波瓣混合器的数值模拟具有重要作用。

图2 涡流结构示意图

目前对波瓣混合器的数值模拟主要采用雷诺平均方法。Koutmos[16]利用k-ε双方程模型模拟计算，并与实验测量结果进行对比，发现可以很好地再现实验中的结果。Cooper[17]等则采用不同的湍流模型模拟了波瓣混合器内的流场结构，验证了不同的湍流模型对实验结果的吻合情况，研究发现k-εRealizable湍流模型相对于其他双方程模型，对波瓣混合器内最大流向涡和正交涡的值吻合的最好。而更多的数值模拟研究集中在优化波瓣混合器的几何外形，进而提高混合效率。O’Sullivan[18]和Du[19]研究了不同波瓣张角对混合效率的影响，Wright[20]研究了不同波瓣数目以及波瓣面积对波瓣混合器增强同向混合效率的影响，Lei[21]研究发现流向涡的破碎过程对于波瓣混合器提高混合效率具有重要作用，因此设计出新的波瓣混合器几何结构用以增强混合效率。

但是，基于雷诺平均方法对波瓣混合器的数值模拟研究只能得到混合器内大尺度的流场特性，而缺少了小尺度的结果和波瓣混合器内的瞬时流场特性，尽管Ooba[22]等采用大涡模拟方法（LES）对波瓣混合器的流场进行了数值模拟，但也只给出了数值的平均结果，对波瓣混合器内小尺度的湍流结构变化缺乏研究，因此波瓣混合器增强同向混合作用的混合机制仍然有待完善。本文采用大涡数值模拟计算，研究了波瓣混合器流场中涡流结构，分析了不同涡流结构的发展过程。

1 数学物理模型及数值方法

本研究中波瓣混合器的具体尺寸参照Hu[14]759等人的实验数据，在建立三维模型时，圆管部分的特征直接通过对圆环进行拉伸得到，而波瓣处的外形特征，由于其复杂的波瓣外形则需要首先通过二维尺寸添加波瓣出口处的外张角和内张角的辅助线，最后通过各个辅助线的位置创建外形曲面，如图3所示，其主要参数如下：D为波瓣混合器内径40mm；H为波瓣高度15mm；d为波瓣宽度6mm；θout-in为波瓣外壁内张角22°；θout-out为波瓣外壁外张角28°；θin-in为波瓣内壁内张角17°；θin-out为波瓣内壁外张角9°。

图3 波瓣混合器三维结构

在计算域的选取方面，对于波瓣混合器的内部计算域，直接采用几何抽取的方法，抽取计算域如图4所示，需要注意的是计算域只包含管道内部流动空间，对于管壁厚度不予考虑。对于波瓣混合器的外部计算域，为了消除边界对于数值计算结果的影响，采用10D的圆柱状，其流向距离为12.6D，最终得到的波瓣混合器的计算域如图4所示。

图4 波瓣混合器计算域

在网格划分方面，考虑边界层的影响，在波瓣混合器壁面处采用边界层网格，第一层网格高度为0.01，网格层数为15层。由于波瓣混合器出口位置射流与周围流体之间存在很强的剪切作用，因而对波瓣混合器出口截面位置进行局部加密，加密长度为波瓣混合器出口截面沿下游方向6D，最终划分的网格数量为38 867 454，如图5所示。

图5 非结构化网格划分

初始和边界条件方面，本文根据Hu[14]760等人实验测得波瓣混合器出口处的速度U0=20 m/s，通过波瓣混合器几何外形上出口与入口面积之比为1.36，流体的入口速度为14.5 m/s。流体介质为空气，密度为1.225 kg/m3，动力黏度为1.789 4 Pa·s。波瓣混合器的外部流场区域采用压力出口边界条件，静压为1个大气压；波瓣混合器壁面采用无滑移边界条件。

数值求解的方程包括不可压缩的连续方程和动量方程。在大涡模拟中，采用密度加权过滤（Favre平均）后得到的控制方程为

式中：“-”代表网格尺度过滤后的物理量，“～”表示为Favre过滤；上标sgs表示亚格子量（subgrid scale）和分别为过滤后的速度分量和压力，v为运动黏度，指小尺度脉动与可解尺度之间的动量输运，本文采用动态Smagorinsky模型对亚格子应力进行封闭，表达式为

式中：δij为克罗内克算子为应变率张量是亚格子涡黏系数为亚格子湍动能；CD和C1是通过局部动态模型得到的模型系数；Δ为过滤尺度。

本文基于流体力学计算软件OpenFOAM开展波瓣混合器的大涡模拟研究，物理模型的流体进口为常温常压的空气，马赫数较小，忽略其可压缩性，求解质量守恒方程和动量守恒方程，采用动态Smagorinsky涡粘模型对亚格子应力进行模化。采用Piso算法求解压力速度的耦合，采用有限体积法对控制方程进行离散，对流项和扩散项使用二阶中心差分格式，而时间导数项使用二阶隐式差分。

2 结果讨论

定义沿X、Y、Z方向上各截面上的无量纲涡量ωx、ωy、ωz分别为

式中：U0为波瓣混合器出口中心线速度，U0=20 m/s，u，v，w均为瞬时速度分量。

2.1 流向涡分布

流向涡的产生对波瓣混合器提高混合作用具有关键作用，流向涡的形成是由于波瓣混合器的特殊波瓣结构。

图6为不同截面上波峰处ωz分布，在波瓣混合器的圆管与波瓣喷管的连接截面上，部分射流流体开始进入到波瓣喷管当中，进入到波瓣处的流体分别向两侧流动，形成了两个旋涡，随着流动距离的增加，由于波瓣张角的存在，使得进入到波瓣内的流体越来越多，从而导致了流向涡越来越大并沿着波瓣截面分布。而后随着距离波瓣混合器出口截面距离进一步减少，流向涡进一步增大，并始终分布在波瓣两侧。可以预见波瓣张角、高度以及宽度直接影响着进入到波瓣内的流体质量，从而影响着流向涡的大小以及其分布形状。随着流向距离的进一步发展，当射流气体通过混合器截面以后，在波瓣混合器的出口截面，尽管外部流体的速度为0，但由于流向涡对周围流体的卷吸作用使得在波瓣混合器的外壁面存在较小的速度，并与波瓣混合器外壁面相互作用，存在较小的涡量。此时，内侧涵道产生沿径向方向向外的二次流动，外侧流体则有向内的二次流，内外两侧涵道诱导产生二次流发展成为流向涡，最终形成了一对反向旋转的涡对。同时，由于在波瓣壁面边界层的存在，在流向涡内靠近波瓣壁面处形成了狭长的沿着波瓣壁面分布与流向涡方向相反的涡量。值得注意的是，在流向涡的形成过程中，在波瓣内部径向向外的压力梯度的作用下，边界层沿波瓣表面向波峰方向卷起，形成了一对反向旋转的涡对，称为马蹄涡，这部分涡量并没有像流向涡一样随着距离的增加逐渐增大，而是始终存在于波瓣顶部，并在波瓣出口截面以后迅速消散，而由于波瓣壁面产生的这部分涡量也在下游位置迅速消散。这两部分涡量的值以及分布相对流向涡较小，且耗散速率较快，因而对提高混合效率的影响较少。

在不同截面处大涡数值模拟得到的流向涡与实验结果的对比如图7所示。从图7中可知，大涡数值模拟很好地反映了流向涡沿流向方向的变化。在波瓣混合器的出口位置，内侧涵道产生沿径向方向向外的二次流动，外侧流体则有向内的二次流，内外两侧流体相互作用，使得流向涡发展成为一对反向旋转的涡对，以波瓣波峰线为对称轴，两侧的涡量几何相等，但符号相反。相对于实验的测量值，大涡数值模拟得到的流向涡更为狭长，分布在波瓣混合器的波瓣两侧。随着流动距离的增加，速度不断下降，流向涡的值并不是变得更小，而是整个涡结构开始破碎，由最初的一对反向旋转的涡对，逐渐分解，形成多个尺度更小的涡结构，在Z/D=1和1.5截面上，大涡数值模拟以及实验测量都再现了流向涡的破碎过程。随着流动距离的进一步增大，可以看到在Z/D=2和3截面上，已经破碎为小尺度结构的流向涡开始逐渐耗散，直至消亡。因而流向涡在波瓣混合器流场中的过程可以分为产生、发展、破碎以及消亡四个阶段。相比于Hu[14]761等人基于雷诺平均方法的计算结果，大涡数值模拟更好地反映了流向涡变化，尤其在流向涡破碎为小尺度的涡流阶段。

图6 不同截面上波峰处ωz分布

2.2 正交涡分布

图7 不同截面上流向涡分布

在内外两侧涵道中，射流与外部流体存在较大的速度差，从而形成较强的剪切层，导致在围绕波瓣混合器的出口截面上形成了正交涡量。图8显示了正交涡的瞬时特征，大涡数值模拟较好地反映了正交涡的分布。相比于流向涡，正交涡的形成是由于剪切层的作用，因而不同于流向涡，正交涡形成在波瓣混合器的出口截面位置，其形状特征与波瓣混合器截面的几何结构一致。同流向涡一样，正交涡也存在着产生、发展、破碎以及消亡四个阶段。在Z/D=0.25截面上，正交涡连接成环状结构，而在波峰与波谷的连接处，正交涡出现了扭曲变形，而随着流向距离的增加，在Z/D=1处截面上，正交涡开始破碎，破碎的位置发生在波峰与波谷相连的扭曲变形处，破碎的正交涡分成两部分，一部分在沿着波峰分布呈现月牙状，另一部分沿着波谷分布形成花冠结构，随着流向距离的进一步增大，这两部分逐渐破碎，正交涡也不再连接成环状，各自分别形成更小尺度的涡量，最终耗散消亡。采用RANS的计算方法尽管捕捉到了正交涡最初的破碎过程，但是并没有反映出之后的衰减情况。值得注意的是，Hu[14]762等人对不同截面上的最大流向涡以及正交涡进行了研究，发现它们具有相同的衰减规律，这表明两种涡结构存在着相互作用。

图8 不同截面上正交涡分布

2.3 速度分布

图9（a）为沿波瓣混合器波峰方向横断截面的速度分布。由于波瓣混合器的出口截面相对于圆截面的面积减少，因而在波瓣喷管处，流体速度将增加。在出口位置的流向方向，速度呈现“三叉”结构，波峰出口处速度与中心射流速度形成三个明显的“尖叉”，这使得剪切层不仅存在于外部流体与射流流体之间，在中心射流与波瓣波峰流体之间的区域也存在剪切层，增强了混合作用。随着流向距离的增加，速度逐渐降低，尤其在波峰处的流体速度降低的最为明显，这是由于这部分流体的速度梯度较大，受到剪切层的影响较大，耗散效应更为明显。在Z/D=2的位置处，从波峰处射流出的流体速度已经基本耗散，这与流向涡和正交涡的消亡阶段保持一致。而在沿波瓣混合器波谷方向横断截面上如图9（b）所示，此时波瓣混合器的外形是向内凹进，因而不会出现“三叉”结构，中心射流仅受到外部流体的剪切作用。

图9 速度分布（t=0.467 009 56）

3 结论

本文基于大涡模拟方法（LES）开展了波瓣混合器涡结构的数值模拟，在波瓣混合器流场中主要存在大尺度的流向涡、正交涡以及小尺度的马蹄涡，通过对比PIV实验结果，大涡数值模拟很好地再现了波瓣混合器内不同涡结构的发展过程。流向涡产生于波瓣混合器的波瓣处，流体进入到波瓣内并产生旋涡，随着进入到波瓣内的流体越来越多，旋涡逐渐增大并形成流向涡，在波瓣混合器出口处，射流流体卷吸次流流体，使得流向涡进一步发展增大，流向涡随着流向距离的增加，破碎成更小尺度的涡量而并不是逐渐衰减变弱。相比于流向涡，正交涡的形成是由于射流流体与次流流体的剪切作用，因而正交涡形成与波瓣混合器出口截面，并呈现波瓣状，同时在波瓣混合器流场中存在多个剪切作用，一是中心射流与周围流速较低的流体之间的剪切作用，二是次流流体与射流之间的剪切作用。