方腔涡流运动对压力脉动噪声影响的数值分析

郝维，刘正先，陈丽英

(1.天津大学 机械工程学院，天津300072;2.天津大学 化工学院，天津300072)

流体流过方腔引起的脉动噪声是实际工程中经常遇到的一类问题，如风管道、涡轮机中的空腔以及水下航行体的空隙和排水孔等.对腔体内部流动机制的研究得到许多研究者的重视[1-4]，近年来，随着计算机技术的飞速进步带动计算流体力学(CFD)快速发展，基于计算流体力学与气动声学交叉学科的计算气动声学得到迅速发展.采用数值模拟的方法研究流体在固体腔中由流动形成的噪声，以及非定常湍流的涡流形成机理，为研究涡流噪声进而削弱和消除噪声提供了基础的数值依据.Powell等[5]认为低马赫条件下的等熵绝热流体，其流体动力场和辐射声场的基本且唯一的源是涡，可以通过研究流场中涡的运动分析其声压变化规律.

本文通过对方腔内流场的数值模拟，分析涡旋产生、发展的规律及与其相应关联位置处压力脉动变化之间的耦合性，采用快速傅里叶变换得到频域下的脉动压力级值，并与实验数据进行对比和验证.

1 数值分析方法

1.1 湍流模型

以方腔为研究对象建立计算模型，通过数值模拟求解湍流N-S方程得到腔体内的湍流流场和涡流运动信息.分别利用标准k-ε模型和大涡模拟(LES)实现对湍流脉动量的模化，标准k-ε模型用于确定流场的定常平均参数，为进一步应用LES方法模拟非定常流动条件下的涡旋产生、发展和耗散提供基础信息.

1.1.1 标准 k-ε模型

标准k-ε模型是目前应用较广泛的两方程湍流模型[6]，它对运动方程中的所有脉动量实施模化，湍流粘度通过求解湍动能k方程和湍动耗散率ε方程得到，k和ε的求解方程如下:

式中:t为时间，s;Xi为空间坐标(i=1，2，3);ρ为流体密度，kg/m3;μ 为流体的动力粘度，Pa·s;μt为湍动能粘度，Pa·s;Gk为由平均速度梯度引起的湍动能产生项;Gb为由浮力影响引起的湍动能产生项;YM为湍流脉动膨胀对耗散率的影响.其中，Gk和YM对不可压缩流动为 0;C1=1.44，C2=1.92，C3=0.09;湍流普朗特系数分别为:σk=1.0，σε=1.3.

1.1.2 大涡模拟

大涡模拟方法(LES)是介于直接数值模拟(DNS)与雷诺应力模型(RSM)之间的一种湍流数值模拟方法.LES假设湍流中的动量、质量、能量及其他物理量的输运主要受大尺度涡影响，其湍流运动通过N-S方程直接求解，而小尺度涡对湍流的影响则通过模化在N-S方程中体现出来.

LES的运动方程通过对N-S方程在波数空间进行滤波得到.过滤原则是去掉比过滤宽度小的涡旋，大涡控制方程如下:式中:ui表示坐标维度下平均速度(i=1，2，3);τij为亚格子尺度应力，

方程式(3)～(4)与雷诺平均方程相似，不同的是方程的变量是过滤后的物理量，而非时间平均量.式中的湍流应力采用Smagorinsky[7]的基本亚格子应力(SGS)模型，该模型在LES方法中应用较广，且取得了有效的结果[8-9].

亚格子应力具有下列形式:

1.2 压力脉动分析方法

1.2.1 压力级的傅里叶变换

数值模拟得到的关键位置点的噪声分析数据来源于确定时间段内的压力脉动信息，该脉动信息是时域下的流动信号，而湍流噪声是频域信息，故采用快速傅里叶变换(FFT)算法[10]实现LES结果从时域到频域空间的转换.

对有限长序列x(n)进行离散傅里叶变换(DFT)表示为

对复数序列x(n)中的一个k值，按上式计算X(k)需要N次复数乘法和(N-1)次复数加法.对所有N个k值，共需要N2次复数乘法和N(N-1)次复数加法.利用旋转因子的对称性)和周期性，将长序列大点数的DFT分解为小点数的DFT，利用多个小点数DFT的计算代替整体的DFT计算，达到降低运算量，明显提高DFT运算速度，缩短运算时间的目的[11].

1.2.2 压力级的倍频程处理方法

与本数值模拟结果进行对比分析的实验数据是由国内某大型水槽实验室通过对相同方腔模型的水动噪声测试提供，实验针对方腔模型的关键位置点进行噪声监测，并对监测数据进行1/3倍频程处理.

数值模拟模型除与水槽实验模型具有相同的几何尺度、流场边界条件外，还应保持相同的数据处理方法.实验采用1/3倍频程数据处理方法，总压力级等效为

式中:i为倍频程带宽内3个1/3倍频程带宽信号，Lpi为倍频程带宽内第i个1/3倍频程压力级测量值.

完成一个标准的1/3倍频程分析，需要分别进行32次测量[12]，因此实验测量得到32个频率及对应的脉动压力级.为便于数值模拟结果与实验数据的比较，对数值模拟数据经傅里叶变换后再进行相同方法的1/3倍频程频处理，使其转化为与实验测量频率同值下的噪声压力级数据.

2 计算模型

方腔计算模型如图1所示，计算区域总尺度为4.3 ×0.9 ×0.8(长 × 深 × 宽，单位:m);其中方腔尺度为0.3 ×0.3 ×0.2;流动介质为水;计算中入口给定均匀流速为7 m/s;对应雷诺数为2.1×106;出口按照不可压缩流动条件满足质量守恒;下边界为不渗透固体壁面，满足无滑移流动条件，上边界设速度与压力梯度均为0.

图1 方腔模型结构及边界条件Fig.1 The cavity model and the boundary conditions

图2为腔体在X-Y与Y-Z流面上的网格分布.在网格生成过程中重点对涡旋生成、发展区域进行加密，而对主流区域采用逐渐过渡方法减少网格总数.根据模型对壁面处网格的要求，壁面第1层网格的无量纲垂直距离y+保持在10左右，整体计算域的网格总数为186万，其中方腔内部网格数为94.5万.非定常流场的计算时间步长为0.000 1 s，共计算10 000个时间步长，即1.0 s间流场的运动过程.

图3为方腔关键点位置，前缘点p1用于考察左侧来流经过腔体前的流动参数变化，腔底p2点用于监测流体在腔体内部流动对底壁面的影响，后缘点p3用于分析来流撞击台阶时流动参数变化.首先分析3点处的脉动压力值，再进行声压分析.

图2 腔体网格分布Fig.2 The grid distribution for the cavity hole

3 数值模拟结果分析

图3 监测点位置Fig.3 The location of the reference points

3.1 流场分析

图4为方腔中间面在t=1.0 s时刻速度矢量图.可看到，图中一个明显的顺时针大尺度涡，同时在腔体下侧两边角处分别有一个逆时针角涡.腔体左侧前缘台阶的存在，使流动在此处出现分离，并在下游形成涡流.不同时刻的涡流运动分析发现(图5)，该涡流在向后运动过程中逐渐发展和合并，涡的尺寸不断增大，最终在腔体右侧后缘台阶处破裂，一部分能量沿水平壁面向右传递，并在后缘诱发产生新的涡旋;另一部分能量则沿垂直壁面向下传递，成为腔内中心大涡形成的主要因素.在流场中，图4的中心涡和所有角涡保持相对稳定的位置和尺度，而前缘和后缘处的涡则进行生成、运动、发展和耗散变化，呈明显的周期变化规律.

图4 腔体t=1.0 s时刻的速度矢量Fig.4 The velocity vector at the end of t=1.0 s

图5为3个不同时刻方腔涡流随时间的形成和发展过程，其中不等间隔时间t1、t2、t3分别取为方腔前缘台阶涡生成、后缘台阶涡到达和破裂时刻.在图5(a)中可明显观察到涡A后方连续的涡旋存在，且涡的尺度不断增大;涡C在t2时刻接近后缘处并具备了分裂趋势，在t3时刻可明显看到分裂后的2个涡分别向水平和垂直方向移动;涡B则为后缘台阶上新的诱导涡.

3.2 压力脉动数值分析

图5 腔体不同时刻涡量Fig.5 The vorticity distribution in the cavity hole on three time points

图6为方腔3个位置点即p1、p2和p3在时间步长为4 000～10 000时段的压力脉动强度分布.可以看出，三者的压力值均显示出周期或准周期波动特征.就波动幅值而言，p3点最大，p2点最小，p1点的周期性效果最不明显;从平均值看，p1点基本在零值附近，p2和p3点均为负值，p3点的负压值最大.分析认为各点的压力脉动特征与其位置相关:p3点位于腔体后缘台阶处，前缘处流体产生的涡旋经生长和合并，在此处破裂，发生较剧烈的变化，同时沿后缘水平壁面再次形成涡流，由此导致压力波动剧烈且幅值较大;p2点位于方腔底部中心，由于腔体的中心涡尺寸虽大但较稳定，因此p2点的压力波动幅值很小;方腔进口端的p1点，受水流入腔时边界条件变化影响，突然从稳定流态演变为剧烈的后台阶涡旋流动，并受中心涡的制约，在台阶下方的垂直壁面附近诱发出逆向小涡，由此形成了虽呈周期性波动但幅值不稳定的压力变化规律.

以p3点的压力脉动为例，提取其波峰值和波谷值对应时刻的涡量位置如图7所示，可以明显看到，图6中的波峰a点对应图7(a)中涡流到达后缘处对壁面形成的挤压状态，而波谷b点正是图7(b)涡旋沿水平和垂直方向分裂为2个涡时的对应状态.

图 6 p1、p2、p3 位置点压力脉动Fig.6 The pressure pulsation for p1，p2and p3

图7 p3点压力峰、谷时刻涡量Fig.7 The vorticity at peak and valley point of pressure for p3

4 涡流压力级分析

4.1 傅里叶压力级

数值模拟得到流场的压力信息通过以下处理得到脉动声压值:以MATLAB程序为平台，确定压力点的采样点数为6 000(由数值模拟数据分析，确定0.4～1.0 s时间段，压力脉动呈现准周期且稳定的波动规律)，采样频率为10 000，对所得压力信号进行傅里叶转换并提取幅值，再进行频率转换得到频谱图，即频率-脉动压力级图.

图8分别为 p1、p2、p3位置点的频率-脉动压力级.可看到，随着频率值的增大，压力级呈下降趋势，同时出现若干峰值点.把大于10 Hz的第1高峰值做为主频值，则 p1、p2、p3点的主频值分别为42、26.67、37.14 Hz.各点的次高和其余主频值见表 1.由压力级峰值对应的频率可以发现，p1的频率普遍大于p2和p3;p3的脉动压力级最大，即涡旋撞击台阶产生的压力脉动对声压值贡献最大.

另外，还考察了p1位置处涡旋随时间的变化过程，发现涡旋在向右发展过程中，不断与周围小涡旋合并、生长，这导致p3的脉动频率减小.

图8 监测点的频率-脉动压力级对比Fig.8 The frequency-SPL diagrams of the reference points

表1 参考点频率与脉动压力级峰值对应表Table 1 The correspondence between frequency and SPL peak value for the key points

4.2 压力级的数值与实验数据比较

实验测量脉动压力数据和数值模拟计算数据均运用前节的噪声处理方法得到关键位置的脉动压力级.图9分别为p1、p2、p3这3个位置点压力级的数值与实验数据比较，其中横坐标为各频率对应的序号Ni，频率对应值见表2.

从图9可以看出，数值与实验结果的总体趋势一致，均随参考点频率值的增大，脉动压力级呈下降趋势.3个压力点的数值与实验数据在40 Hz之后呈现很好的符合度，在40 Hz之前则一致地表现出模拟值大于实验值.

图9 压力级数值与实验数据比较Fig.9 Comparison between simulation and experiment

表2 频率与Ni对照表Table 2 The correspondence between frequency and Ni

5 结论

1)通过对方腔内部涡流流场的非定常数值模拟和分析，确定腔体不同位置点处涡旋产生、发展和破裂运动具有脉动周期性，压力脉动的峰、谷特征与涡流状态具有良好的对应性.

2)采用快速傅里叶变换将数值计算的时域信息转换为压力脉动频域数据，通过1/3倍频程处理方法得到各频率下的压力级值.方腔不同的位置点具有不同的压力脉动频率和压力级特征，与其涡流运动密切相关.

3)与实验测量结果的对比表明，除40 Hz以下低频区外，两者具有很好的符合度.表明采用数值方法模拟方腔内的涡流流场，进而结合快速傅里叶变换和倍频程处理方法实施对方腔内涡流引起压力脉动的研究是可行的，可以为方腔类的水动力噪声和消声控制提供参考.

[1]耿冬寒，刘正先.大涡模拟-Lighthill等效声源法的空腔水动噪声预测[J].哈尔滨工程大学学报，2010，31(2):182-187.GENG Donghan，LIU Zhengxian.Predicting cavity hydrodynamic noise using a hybrid large eddy simulation-Lighthill＇s equivalent acoustic source method [J].Journal of Harbin Engineering University，2010，31(2):182-187.

[2]朱习剑，何祚镛.水洞中突出矩形腔的流激驻波振荡研究[J].哈尔滨工程大学学报，1993，14(3):41-52.ZHU Xijian，HE Zuoyong.Study of flow-induced standing wave resonance of rectangular cavity in water tunnel[J].Journal of Harbin Engineering University，1993，14(3):41-52.

[3]GHAR I M，ROSHKO A.The effect of flow oscillations on cavity drag[J].Journal of Fluid Mechanics，1987，177(10):510-530.

[4]WANG M，FREUND J B，LELE S K.Computational prediction of flow generated sound[J].Annual Review of Fluid Mechanics，2006，38(1):483-512.

[5]WANG Chunxu，ZHANG Tao，HOU Guoxiang.Noise prediction of submerged free jet based on theory of vortex sound[J].Journal of Ship Mechanics，2010，14(6):670-677.

[6]王福军.计算流体动力学分析-CFD软件原理与应用[M].北京:清华大学出版社，2004:116-123.WANG Fujun.Computational fluid dynamics analysis-principle and application of CFD software[M].Beijing:Tsinghua University Publisher，2004:116-123.

[7]张兆顺，崔桂香，许春晓.湍流大涡数值模拟的理论和应用[M].北京:清华大学出版社，2008:57-59.ZHANG Zhaoshun，CUI Guixiang，XU Chunxiao.The theory and application of large eddy simulation of turbulent flows[M].Beijing:Tsinghua University Publisher，2008:57-59.

[8]GERMANO M，P IOMELL I U，MO IN P，et al.Dynamic sub-grid-scale eddy viscosity model[J].Physics of Fluid A，1991，3(7):1760-1765.

[9]LILLY D K.A proposed modification of the Germano subgrid-scale closure model[J].Physics of Fluid，1992，4(3):633-635.

[10]张广超，宋文爱.基于递归法的FFT计算机仿真[J].国外电子测量技术，2008，27(6):9-11.ZHANG Guangchao，SONG Wenai.Computer simulating of FFT based on recursion[J].Foreign Electronic Measurement Technology，2008，27(6):9-11.

[11]谭子尤，张雅彬.离散傅里叶变换快速算法的研究与MATLAB算法实现[J].中国科技信息，2006(22):316-317，321.TAN Ziyou，ZHANG Yabin.Research of fast Fourier transformation and realization of MATLAB[J].China Science and Technology Information，2006(22):317-321.

[12]杨昌棋，秦树人，何辉.基于FFT的虚拟实时噪声倍频程分析仪[J].测控技术，2000，19(9):25-27.YANG Changqi，QIN Shuren，HE Hui.Virtual real-time noise octave analyzer based on FFT[J].Measurement＆Control Technology，2000，19(9):25-27.