钝体绕流现象广泛存在于自然界和工程应用中．一方面，研究者们通过分析作用在结构上非稳态力的特性达到控制涡激振动(VIV)的目的[1]；另一方面，非稳态力也会产生流激噪声，对声音的产生和传播进行研究，从而控制有害噪声，有着重要的学术和工程意义[2]．此前，较多研究者基于声类比采用有限元和边界元方法求解波动方程进行声辐射研究[3-7]，但是该方法不考虑流场对声音的影响，且采用等效声源，较难对声音的产生机理进行细致的分析，而基于声传播方程的计算声学混合方法在这方面展现了较大的优势．较早对钝体绕流噪声进行研究的学者Strouhal发现圆柱的辐射噪声频率和均匀来流速度存在一定联系，提出著名的斯特劳哈尔数(Sr)．近年来，较多研究者采用数值方法对结构绕流噪声进行分析．文献[8]求解可压缩纳维-斯托克斯方程对均匀流中二维圆柱发声机理进行详细分析并讨论多普勒效应对声传播的影响．部分研究者提出不同的声传播模型且用圆柱绕流噪声问题进行验证，如文献[9]提出的改进黏声分离方法及文献[10]推导的声扰动方程(APEs)等．也有部分学者对方柱绕流噪声进行研究．文献[11]采用混合方法讨论了较高马赫数下(Ma=0.3)方柱的声传播特性．文献[12]讨论了长宽比为7的矩形柱在0°迎角下的声学特性．文献[13]讨论了0°迎角下并联双方柱间距对流场和声学特性的影响．以上研究对钝体绕流声学特性分析有着重要的参考价值，但是矩形柱和圆柱又有着一些不同，主要是迎角和长宽比对其流噪声的影响，之前的研究侧重0°迎角和固定长宽比情况，在变迎角和变长宽比方面的声学研究还较少且缺乏系统性．基于上述背景，本研究提出一种混合声学计算方法对矩形柱绕流噪声进行分析．首先采用不可压缩直接数值模拟对非定常流场进行求解，再以流场结果作为声源输入声扰动方程获得其声学特性．流场和声学的固体边界条件都采用浸入边界方法处理．利用该混合方法开展了低雷诺数下矩形柱的流场和声学特性研究，讨论了方柱的迎角和矩形柱的长宽比对流噪声的影响规律，为后续矩形柱的发声机理和降噪控制研究提供参考．1 数值方法1.1　控制方程流场的控制方程是黏性不可压缩纳维-斯托克斯(NS)方程，其矢量形式为∇•U=0；(1)∂U/∂t+(U•∇)U=-∇P/ρ0+v0∇2U，(2)式中：ρ0，U，P分别为不可压缩流体的密度、速度矢量和压力；t为时间；v0为流体运动黏度．由于黏性对声传播影响较小，因此声学计算时假定流体是无黏的．将压力、速度和密度分解成时均量和扰动量两部分，可以得到声扰动方程[10]为∂u'/∂t+∇(U¯•u')+∇(p'/ρ¯0)=Smom；(3)∂p'/∂t+c¯02∇•(ρ¯0u'+U¯p'/c¯02)=c¯02Scont，(4)式中：c¯0为时均声速；u'和p'分别为速度扰动和声压；U¯和ρ¯0分别为流场的时均速度和时均密度．式(3)和(4)的等号左边描述了声波在非均匀流场中的传播和反射过程，等号右边为声源项．对于低马赫数问题，声源项为Smom=∇P'/ρ0[10]，其中P'为不可压缩压力扰动．1.2　离散格式流场的求解采用文献[14]提出的浸入边界-有限差分程序．该方法基于笛卡尔同位网格，采用有限差分方法直接离散不可压缩NS方程．空间离散采用二阶中心差分格式，时间上采用二阶精度的分步法进行推进．该方法被广泛应用于求解具有复杂几何和动边界的流动问题．声场空间离散采用文献[15]提出的7点4阶色散关系保持格式(DRP)，该格式经过优化具有低色散和低耗散的特点，被广泛应用于气动声学的研究．时间离散采用文献[16]提出的优化龙格库塔格式，该格式也满足低耗散低色散特性．高阶差分会产生高频振荡，为了保持数值稳定性，须增加高阶滤波．本研究采用11点10阶滤波[17-18]，在边界处降阶处理．详细的声学方法验证可参考文献[19]．1.3　边界条件本研究采用的边界条件有两种，一种是固体壁面边界条件，另外一种是计算域边界条件．流场计算域边界条件为：入口和侧面都采用速度入口边界，而出口的速度梯度为零，所有压力边界条件采用诺伊曼(Neumann)边界条件，即压力梯度为零．对于声学计算须设置无反射边界条件．本研究采用文献[15]推导的辐射边界条件，在计算域边界外部构建三个虚拟网格并采用单侧DRP格式进行离散求解．流场计算的固体边界须满足无滑移边界条件．在声学计算时假设流体是无黏性的，满足无穿透边界条件．本研究的流场和声学固体边界条件都是基于单个虚拟网格的浸入边界法对边界进行计算．首先对整个计算域采用均匀或者非均匀笛卡尔网格离散．依据结构边界和流体的位置将网格分为3类：在固体外部的是流体网格；在内部的是固体网格；第三类是虚拟网格，虚拟网格本质在固体内部，但是至少有一个相邻的流体网格．在搜寻完所有虚拟网格以后，以结构表面为镜面找到虚拟网格在流体中的镜像点．镜像点单元的物理量须通过周围的网格插值得到，然后虚拟网格的值通过各自的边界条件计算得到．固体边界条件处理方式示意图如图1所示，详细的浸入边界方法介绍和验证可参考文献[14]．10.13245/j.hust.211017.F001图1固体边界条件示意图2 计算模型和参数流场和声学均采用无量纲形式进行计算．计算模型如图2所示，矩形柱高度为D，长度为W，计算参考长度为D，参考时间T=D/U0(U0为来流速度)，参考声压为ρ0c02，参考涡量为U0/D，雷诺数(U0D/ν0)为150，来流马赫数(U0/c0)为0.2．本研究的流场和声场采用相同的笛卡尔网格，须指出的是：在实际应用中，因为流场和声学计算对计算域和网格要求不同，所以流场和声学网格和区域可以不同，此时需要网格插值处理．10.13245/j.hust.211017.F002图2计算模型示意图划分网格时须要同时考虑流场和声场对网格的要求．对于流体计算，须要对固体周围和尾部加密捕捉边界层和尾涡．对于声学模拟，采用DRP格式，须要满足一个波长内至少有6~8个网格[15]．考虑到计算成本，本研究使用非均匀笛卡尔网格，选取3组不同网格进行无关性验证，最小网格尺寸dx分别为0.03D，0.015D和0.01D．当沿着矩形柱向计算域边界网格尺寸逐渐增大到1.6D时，仍然将均匀网格用于声波的模拟．此时每个波长内网格数约为20个，满足声学模拟的要求．整个计算域尺寸为0≤(x，y)≤200D，其中x和y分别为横轴坐标和纵轴坐标，采用较大的计算域用于声传播模拟并减少边界反射．流场的无量纲时间步长为0.01D/U0，声学的无量纲时间步长为0.005D/c0，声学模拟在每一时刻的声源输入通过与该时刻相邻的两个流场时刻的压力扰动线性插值计算得到．3 计算结果和分析3.1　数值模型验证为了验证本文数值模型，首先对方柱在0°迎角下进行计算，3组不同尺寸网格下方柱的升力系数和阻力系数时间曲线如图3(a)所示，方柱的升力系数(Cl)和阻力系数(Cd)曲线都随时间周期性波动，随着网格细化差异逐渐减小，表明网格是收敛的．此外，本研究的平均阻力系数为1.401，和文献[20]的计算结果(1.406)及文献[13]的计算结果(1.415)接近．升力系数的均方根为0.280，和文献[21-23]结果较为符合(0.162~0.292)．流场的涡量分布如图3(b)所示，可以看到方柱尾部上下交替脱落的卡门涡街．对升力系数曲线进行快速傅里叶变换(FFT)，得到斯特劳哈尔数为0.150，和文献[13]的仿真结果(Sr=0.151)符合较好．综上结果表明：本文流场计算方法和参数是有效的，可以作为声学计算的输入条件．10.13245/j.hust.211017.F003图3方柱流场计算结果图4为方柱声学计算结果．图4(a)展示了声压云图分布，可以看到方柱绕流产生的噪声以上下对称的形式向计算域四周传播，具有明显的偶极子声源分布形式．在方柱后侧区域有一带状以交替脉动形式向外传播的压力扰动，该扰动不以声速传播，而是沿着流场下游逐渐耗散，即拟声现象．此外，由于多普勒效应的影响，因此最大声压方向并不垂直于来流而是偏向上游(图中虚线所示)[8]．为了定量验证本文结果，计算远场半径r=75D处的声压指向性，并和文献[13]的计算结果对比如图4(b)所示，指向性呈现“8”字形分布再次验证了偶极子声源的本质，本文计算结果和文献[13]的直接数值模拟结果符合度较好，定量验证了混合模型的有效性，可用于后续结构的流噪声特性预报分析．为了减小计算量，最终本文的最小网格尺寸选为0.015D．10.13245/j.hust.211017.F004图4方柱声学计算结果3.2　不同迎角下方柱声学特性在实际应用中，流体经常从不同方向流过结构，影响流场和声学特性．本研究以方柱模型为研究对象，分析迎角对其流噪声的影响规律，考虑到方柱结构的对称性，只计算方柱在4种旋转角(0°，15°，30°和45°)下的声传播特性，计算模型如图5所示．10.13245/j.hust.211017.F005图50°和15°迎角方柱模型示意图图6显示了不同迎角下方柱在r=75D处的指向性分布．结果表明：不同迎角下方柱的指向性分布仍然保持着偶极子声源的特性，但是声传播有着较大的变化，说明迎角对声学特性有着重要的影响，主要包括对声压幅值和传播方向的影响．迎角显著增强了辐射声压级，随着迎角的增大，方柱上半部分的辐射声压先增大后减小，在30°迎角达到最大，而下半部分的声压随着迎角的增大逐渐增大，在45°达到最大．在15°和30°情况下，由于结构上下的不对称性，导致声压幅值也不对称，因此方柱上半部分的幅值较下半部分大．在45°迎角上下部分仍然呈现出对称的指向性．此外，迎角对声传播方向也有着较大的影响，在0°情况下，传播方向主要受到多普勒效应的影响，偏向来流上游传播，在15°和30°情况下，声传播角受到迎角和多普勒效应的双重影响，而迎角的存在导致了上下的不对称，在45°情况下，由于结构的上下对称，因此传播方向和0°迎角类似．10.13245/j.hust.211017.F006图6不同迎角下指向性分布3.3　不同长宽比下矩形柱声学特性研究矩形柱长宽比AR(W/D)对声学特性的影响，选取3组不同长宽比的模型(2:1，3:1和4:1)进行计算分析，半径r=75D处的声场指向性如图7所示，方柱作为参考对比．由图可知：矩形柱的声场指向性仍然保持着偶极子声源的特性，但是矩形柱的长宽比对辐射声压有着重要的影响，声压幅值都比方柱模型较小．随着长宽比的增大，矩形柱声压级逐渐减小．另一方面，矩形柱的声传播角变化较小，和方柱类似，多普勒效应的影响使传播方向更加偏上游．对于本文的马赫数(0.2)，最大声压所在位置(声传播角)约为105°，和文献[8]计算的单圆柱传播角为101.5°较为接近．10.13245/j.hust.211017.F007图7不同长宽比矩形柱指向性分布为研究远场辐射噪声特性，选取长宽比为4:1的矩形柱，记录观测点(100D，75D)声压变化曲线如图8所示，结果表明声压达到稳定并随时间做周期性变化．对声压时间历程曲线和升力时间曲线分别进行FFT变换，得到矩形柱的功率谱密度曲线如图9所示，图中：wSD为功率谱密度，代表单位频带内的功率(均方值)．结果表明：声压曲线的频域主频率为0.133，等于尾涡脱落的频率(Sr)，比方柱(0.150)略小，表明尾涡脱落和升力变化是主要的噪声源．10.13245/j.hust.211017.F008图8观测点处声压变化曲线10.13245/j.hust.211017.F009图9声压和升力功率谱密度曲线不同长宽比下矩形柱最大升力时刻的尾涡分布云图如图10所示，流体绕过矩形柱前端发生流动分离，在尾部产生交错排列的卡门涡街，但是涡街形成的位置不同．在方柱情况下，涡形成的位置约为方柱中心后4D(见图3(b))，随着矩形柱长宽比的增加，涡街形成的位置逐渐偏向下游，当AR=4时，形成点约在10D处，且尾涡的强度逐渐下降．10.13245/j.hust.211017.F010图10不同长宽比矩形柱涡量云图不同长宽比矩形柱的升力变化曲线如图11所示，由图可知：尾涡受到矩形柱的影响导致升力系数逐渐下降，波动逐渐减小，并导致了辐射声压的下降．10.13245/j.hust.211017.F011图11不同长宽比矩形柱升力变化曲线4 结论本研究采用浸入边界-声扰动方法对矩形柱声学特性进行分析，主要结论如下．a. 方柱计算模型的流场和声学计算结果与文献符合度较好，说明本文混合计算方法能够有效求解矩形柱绕流噪声问题．b. 迎角对方柱辐射噪声有较大影响，在30°迎角下辐射声压较大；对称结构(0°和45°)的传播角受到多普勒效应的影响，非对称结构受到迎角和多普勒效应的双重影响．c. 矩形柱绕流噪声低于方柱，并随着长宽比增大逐渐减小，主要是因为细长矩形的卡门涡街形成位置向下游偏移，涡脱强度和升力系数下降．d. 矩形柱尾涡脱落和升力脉动变化是主要的噪声源．