基于混合雷诺平均/高精度隐式大涡模拟方法的高升力体气动噪声模拟*

葛明明 王圣业 王光学2) 邓小刚†

1) (国防科技大学空天科学学院,长沙 410073)

2) (中山大学物理学院,广州 510275)

发展了基于七阶精度混合型耗散紧致格式(HDCS)的混合雷诺平均(RANS)/高精度隐式大涡模拟(HILES)模型(HRILES),并结合Ffowcs Williams-Hawkings (FWH)声比拟方法对30P30N 高升力体气动噪声问题进行了模拟.首先对雷诺数 Red=4.3×104 的单圆柱绕流算例开展验证,并与传统的延迟分离涡模拟(DDES) 模型进行对比.结果表明HRILES模型具有对亚临界态尾迹区转捩流动的模拟能力,平均阻力系数与阻力均方根值和实验结果符合更好,结合FWH 声比拟方法得到了合理的远场声压级(SPL)的功率谱密度(PSD)分布.然后将其应用于30P30N 高升力体气动噪声算例模拟,结果表明HRILES模型准确预测缝翼凹腔剪切层各站位的平均速度、涡量和湍动能分布,壁面脉动压力谱分布与实验符合较好,近、远场噪声频谱准确预测了缝翼低频窄带噪声,并得到了合理的噪声辐射指向性分布.

1 引 言

随着航空运输业的快速发展,航空器噪声问题受到广泛关注.欧盟早在2000年就在European Visions 2020中提出降低航空器噪声50%的目标计划[1].美国随后也在AST (advanced subsonic transport)中提出安静飞行计划.伴随着大涵道比发动机降噪技术的发展,机体噪声的研究越来越受到重视,而在机体噪声中高升力体部件贡献了起降期间的大部分声能.一般认为高升力体前缘缝翼的噪声辐射是主要声源[2],缝翼噪声主要部分是宽频噪声,同时宽频中会叠加多个窄带峰值(NBP),此类高雷诺数气动噪声问题对数值模拟能力提出较高要求,需要发展兼顾计算精度与计算效率的计算模型与算法[3-5].

根据NASA 2030报告,大涡模拟方法将是航空工业下一代主流的湍流模拟方法[6].传统的大涡模拟(LES)模型存在一些尚未被解决的问题,如亚格子项存在被截断误差掩盖的可能性、空间滤波中的不确定度、在高雷诺数下的壁面流动问题中难以建立亚格子模型等[7].隐式大涡模拟作为一种简单易用方法,被应用于各类问题计算,其中包括各向同性耗散湍流[8]、 自由剪切流[9]、 壁面流动[10],超声速流动[11-13]、流动主动控制[14]等.高精度隐式大涡模拟(HILES)是基于HDCS-E8T7格式发展的一种隐式大涡模拟方法,利用格式的截段误差代替显式的亚格子模型,同时满足几何守恒律[15],已经成功地应用于各类复杂外形流动与噪声预测[16-18],但对于高雷诺数壁面流动模拟,计算资源的需求仍然较大.

目前湍流模型实际应用仍然以非定常雷诺平均方程(URANS)、混合RANS/LES模型以及壁面建模的大涡模拟(WMLES)为主[19].混合RANS/LES类方法的发展通常分为分区与不分区两种思路,分区方法人为将计算域设置为不同的区域采用不同的模型计算,流场信息通过交界面进行交换.不分区方法对全域统一求解,利用混合函数或者改变源项实现不同模型之间的切换.前者的关键在于处理交界面边界条件,而后者在于交界面不确定性带来的灰区问题.不分区方法中最著名的是Spalart等[20]发展的分离涡模拟(DES)模型,经典的DES模型引入网格尺度来修改Spalart-Allmaras模型(SA模型)破坏项中壁面距离,从而实现RANS向LES的转换,在此基础上改进发展出延迟分离涡模拟(DDES)、改进的延迟分离涡模拟(IDDES)等方法[21,22].Nichols[23]提出了构造混合模型的另一种思路,根据局部网格尺寸和湍流长度尺度从RANS平滑过渡到LES,称之为多尺度模型.随后将其应用于三维圆柱和武器填埋仓的流动模拟中,结果表明该混合方法在湍流尺度与网格尺度的比率接近2时,实现了较为合理的模拟结果[24].

本文利用Nichols[23]提出的多尺度模型构造思想,结合剪切应力输运(SST)模型和HILES方法,希望实现对于高雷诺数流动的绝大部分区域使用HILES精确求解大尺度流动,仅边界层内层利用SST 模型充当壁面模型.本文第2节主要介绍混合节点半节点型耗散紧致格式(HDCS-E8T7)与HRILES 方法,第3节模拟亚临界态的圆柱流动以验证HRILES方法,第4节采用HRILES和IDDES模型针对高升力体噪声问题进行模拟,并将结果与实验数据进行对比分析.

2 数值方法

2.1 高精度数值格式

本文所有的算例都采用七阶混合型耗散紧致差分格式HDCS-E8T7,该格式采用显式八阶差分格式对对流项进行空间离散,通过在中心插值的基础上增加奇数阶导数来调高数值耗散,抑制数值色散,具体推导过程见文献[25].对流项差分格式表示为

其中 α 为可调节耗散参数,根据色散关系保持格式优化取值为0.3085.黏性项采用八阶中心差分格式加线性插格式进行计算.时间推进采用二阶隐式双时间步方法以节约计算资源.HDCS 格式已经应用于阵列圆柱、串联圆柱-翼型、喷流噪声[16-18]、三角翼[26]等典型算例中,为复杂外形气动声学问题模拟奠了定基础.

2.2 HRILES模型

由于RANS方程与LES方程的相似性,可以通过过渡函数实现不同模拟方法的切换.Piomelli等[27]提出了基于SST的多尺度模型,对Red=8×106的圆柱绕流进行了模拟,认为多尺度模型相比于DES模型的优点主要在于能够基于当地网格尺度平滑地从RANS方法过渡到LES方法.本文采用Morgan优化的过渡函数：

其中

(5)式中的 LT和 LG分别代表湍流长度尺度和网格长度尺度,原始文献中定义如下：

本文计算中网格长度尺度参照Spalart等[21]在IDDES中的建议设置为

可以计算出湍流黏性系数 νt为

对于隐式大涡模拟,可认为亚格子黏性由数值格式隐式提供,因此湍流黏性系数可以简单表示为νt=νtRANSfd.实际计算中,RANS模型只在边界层内层启动,大部分边界层流动采用隐式大涡模拟,根据Wang等[26]的综述,这类混合方法也可以看作为WMLES方法的一种.混合模型构造思路是想通过尺度判据将边界层内层采用RANS 模型求解,它最显著的优点是可以推广到任意的RANS模型和SGS模型的组合.SST作为广泛应用的工程湍流模型在高雷诺数复杂外形流动中表现良好,HILES方法具有高精度高效率的特点,HRILES模型可以结合二者的优势.

2.3 噪声计算方法

对于气动声学问题,远场噪声计算通常是采用声比拟方法对近场计算得到的声源面的时间序列数据进行积分,其中FWH方法最为常用.声源面可以采用壁面或者空间任意可穿透面,区别在于后者包含了可能存在的四极子体声源.由于本文算例属于低亚音速流动,四极子声源占比可以忽略不计,故采用固体壁面作为声源面,采用Francescantonio提出的第二KFWH方程,积分求解远场噪声[28]：

其中

(11)式中的下标ret代表积分项对应的是推迟时间的数据,S代表积分声源面.方程(10)等式右边前两项代表厚度噪声(单极子声源),后两项代表载荷噪声(偶极子声源).对于二维构型的算例,由于网格量的限制,往往展向长度小于实验长度,这会影响远场噪声结果,需要对声学结果进行校正.本文采用紧致声源校正常用的Kato公式[29]：

其中SPLexp为校正后声压级,SPLsim为数值积分得到的待校正声压级,Lsim,Lc和 Lexp分别代表数值模拟展向长度、流动展向相干长度、实验展向长度.

3 单圆柱绕流

3.1 计算设置

采用HRILES模型对单圆柱算例对进行模拟,基于圆柱直径d的雷诺数为 Red=43000 ,来流马赫数0.21.该条件下流动处于亚临界态,包括边界层分离,尾迹区转捩以及大尺度涡结构等复杂现象[30],对复杂流动结构的模拟精度直接影响阻力计算结果.计算采用O型网格进行,其中周向和法向分布为 181×181 ,展向均布31个点,壁面法向第一层网格无量纲高度 h/d=1.0×10-5.时间推进采用双时间步法,真实时间步长 dt=1×10-6s.计算100 个对流周期得到近似拟序结构,又计算了200个周期用于统计平均与声源面信息收集.

3.2 流动结果

表1中给出了流动结果的统计数据,并与经典SST-DDES模型的结果和Szepessy和Bearman[31]的实验结果进行了对比.从阻力系数来看,相比SST-DDES方法的结果,HRILES方法得到的平均阻力系数 CD,ave和阻力均方根值CD,rms都更加接近实验值.升力系数的斯特劳哈尔数Sr表征圆柱尾部的涡脱落频率,本文计算得到的结果与Seo等[32]的LES 计算结果更加接近,采用DDES模型会高估壁面的涡脱落频率.θsap代表分离起始位置角度,θ 计算以来流方向为起点,顺时针为正,两种模型计算得到的分离位置接近.

图1展示了HRILES模型计算得到的壁面平均压力系数 Cp分布结果,与实验结果和SST-DDES模型的计算结果进行了对比.HRILES模型计算出的整个背风区压力系数分布都要比SST-DDES模型结果更接近实验值.从图2给出的统计平均得到的流场流线分布来看,HRILES 模型与SSTDDES方法的结果差别较大,前者能捕捉到背风区尾迹的中等尺度的分离泡,而SST-DDES计算结果只得到了大尺度尾部回流区.

表1 单圆柱算例流动参数统计结果Table 1.Statistical results of aerodynamic coefficients for the single cylinder.

图1 圆柱表面平均压力系数分布Fig.1.Mean wall pressure coefficient distribution of the rod.

图2 流线分布 (a) HRILES; (b) SST-DDESFig.2.Distribution of streamlines：(a) HRILES; (b) SSTDDES.

3.3 噪声结果

Jacob等[33]在德国宇航中心(DLR)的声学风洞进行了相似条件下的实验研究,测量得到圆柱正上方,距离 r=185d 处观测点的声压级功率谱密度分布.本文声学结果计算采用FHW方程对壁面瞬态流动数据积分得到,并利用Kato公式修正展向长度差异的影响.每个非定常时间步采集壁面数据,共得到约3万个采样点,这些数据被分为7个窗,数据重叠率 50%,最终功率谱密度的分布通过平均7个窗的结果得到.图3对比了两种模型得到的对应实验远场观测点的声压级功率谱分布,r代表观测点至圆柱中心的距离.结果表明HRILES方法能够准确捕捉尖频噪声,并且主频幅值相较SST-DDES模型更接近实验结果.

图3 远场 θ=90° ,r=180d 观测点声压级功率谱密度Fig.3.Farfield acoustic result of the rod：PSD at (θ=90°,r=180d).

4 高升力体

4.1 计算设置

算例几何外形与机体噪声会议中给出的30P30N缝翼噪声标准算例一致,气动弦长C为18 in (1 in=2.54 cm).缝翼与襟翼的长度分别为0.15C与 0.3C.计算采用等效飞行条件,襟翼缝翼张开角度均为 30°,来流马赫数0.17,基于弦长的雷诺数 ReC=1.7×106,攻角 5.5°.远场边界采用特征边界条件以消除边界声波反射影响,展向边界采用周期性边界条件.

计算网格见图4,采用日本宇航中心(JAXA)提供的多块对接结构网格,计算域沿壁面法向方向延伸100个弦长.JAXA关于网格具体介绍参考文献[34].本文为验证方法,壁面法向网格雷诺数(y+)最大值约1.5.二维网格总量70000,展向长度设置为2 in,均布61 个网格点,空间网格总量为4200000.

图4 30P30N计算网格Fig.4.Mesh of 30P30N airfoil.

4.2 流动结果

图5展示了IDDES模型和HRILES模型计算得到的瞬态Q判据等值面(QC/U∞=5000),U∞为来流速度.Q等值面通过瞬态无量纲展向涡量 ωz进行着色,参考量为 U∞/C.结果清晰展示了缝翼凹腔涡流的尺度范围和剪切层发展轨迹.缝翼下侧凹腔内外壁面流动在过尖缘后混合,形成剪切流动,并迅速转捩为湍流.湍流结构在空间充分发展后,最终在凹腔上壁面再附,再附点是数值模拟的重要指标.流动再附后部分流动向下游发展,部分沿上壁面回流进入凹腔形成涡流.IDDES和HRILES 模型模拟出的剪切层是相似的,但是HRILES模型能分辨出更精细的涡结构,特别是在缝翼凹腔内部.

图6展示了x-y截面平均展向涡量 ωz,ave云图,计算结果沿展向进行了平均.两种模型计算结果分布相似,IDDES模型得到的剪切层涡量峰值比HIRLES 的更大.Pascioni[35]的实验测量了凹腔剪切层七个不同站位的流动参数分布,在图6(b)中分别对应 L1- L7.图6(b)中 P1- P6代表高升力体表面非定常压力脉动数据测点.图7的x-y截面平均流线分布展示了高升力体流动中典型的流动特征,即缝翼和主翼凹腔中的两个回流区.本文计算状态条件下背风面均未出现流动分离现象.

图5 QC/U∞=5000 等值面 (a) IDDES; (b) HRILESFig.5.The isosurfaces of the Q-criterion (QC/U∞=5000)：(a) IDDES; (b) HRILES.

图6 平均展向涡量云图 (a) IDDES; (b) HRILESFig.6.Contours of meanmean spanwise vorticity：(a) IDDES; (b) HRILES.

图7 平均流线分布 (a) 缝翼; (b) 襟翼Fig.7.Distribution of streamlines：(a) Slat; (b) flap.

图8给出了翼型表面平均压力系数的分布结果与JAXA的实验结果的对比[36],统计了0.05-0.10 s时间段的非定常数据,并沿展向进行了平均.HRILES模型结果大部分位置与JAXA的实验结果符合得很好,缝翼上表面负压略微偏低.IDDES模型的结果缝翼上表面明显低估了背风区吸力效应.这与前人的模拟结果规律相似,Terracol等[34]采用DDES和分区DES(ZDES)模拟的结果也低估了背风面吸气效应,他指出可能是由于缝翼上表面边界层在当前雷诺数条件下是层流,使用全RANS 模型模拟引入的误差; Gao等[37]采用通量重构(FR)方法进行的模拟认为是其计算几何模型没有考虑缝翼尾缘厚度造成的; Zhang等[38]的模拟结果也表明缝翼上表面压力系数脉动的网格敏感性强于下表面.由于本文采用的也是JAXA的网格,几何构型相同,因此得到与其类似的结果应该是由于IDDES 模型在背风面边界层大部分区域也采用了RANS模拟.图9给出了缝翼表面压力系数脉动均方根的分布,并对比了实验给出的 P2- P6五个点的数据.两种模型的结果得到的再附点的位置(压力脉动峰值处)接近,但HRILES模型得到的均方根值比IDDES模型更接近实验值.

图8 壁面压力系数分布Fig.8.Distribution of wall pressure coefficient.

图9 缝翼表面压力系数脉动均方根分布Fig.9.RMS of the fluctuating pressure coefficient on the surface of the slat.

图10 各个站位的平均速度分布Fig.10.Mean velocity magnitudes along the seven lines across.

图11给出了缝翼凹腔剪切层对应站位的平均展向涡量 ωz,ave的结果.在剪切层起始位置,两种模型的计算结果都要高于实验数据,而随着剪切层向下游充分发展,计算结果与实验结果逐渐接近.近壁面剪切层的计算误差可能是由于模拟剪切层转捩现象滞后,剪切层未充分掺混.

图12给出了缝翼凹腔剪切层对应站位的平均湍动能TKEave结果,整体上HRILES模型的结果明显更加接近实验数据.初始站位两种模型结果都偏低,这意味着计算出的剪切层相比实验存在转捩推迟现象.不同的是HRLES模型在 L2站位后与实验数据符合较好,而IDDES模型的湍动能始终偏低直至剪切层下游 L6站位,这也说明HRILES模型在空间剪切层发展的模拟上更具优势.

图13给出了壁面非定常压力脉动功率谱密度分布.P1点位于主翼下表面,处于流场的稳态区域,此处的压力谱可以近似表征噪声结果.实验数据中存在两个窄带噪声峰值,对应的频率分别是1330与1960 Hz,前人的研究认为这种低频窄带噪声来源于凹腔内部的涡-声耦合自激励反馈机制.两种模型的脉动压力谱分布相似,但IDDES模型计算出的尖频噪声比实验值偏高8 dB左右,而HRILES模型的结果与实验数据非常接近.高频宽带噪声则来源于缝翼后缘的涡脱落,其脱落频率不仅取决于尾缘厚度还与边界层位移厚度相关.P4点位于充分发展的湍流剪切层流动再附点附近,实验结果分布表明该点为宽频噪声.两种模型计算出的宽带噪声分布量级与实验值符合较好,但是IDDES模型模拟出了显著的窄带峰值,HRILES模型基本保持了宽带特性.在频率高于10 kHz的部分,P1和P4点的计算结果存在幅值显著衰减现象,这可能是由于空间网格密度不够造成的,但由于远场噪声声压级统计要求只对256 Hz-10 kHz频段进行宽带滤波,因此对远场噪声计算结果影响不大.

图11 各个站位的平均展向涡量分布Fig.11.Mean spanwise vorticity along the seven lines across.

图12 各个站位的平均湍动能分布Fig.12.Mean turbulent kinetic energy along the seven lines across.

图13 脉动压力功率谱密度分布 (a) P1 ; (b)P4Fig.13.Frequency spectra of pressure fluctuations：(a) P1 ; (b) P4.

4.3 噪声结果

图14展示了HRILES模型计算得到的x-y截面瞬态压力脉动云图,采用来流动量q进行无量纲.图中缝翼噪声的强度远大于其他区域,并有显著的指向性特征,计算结果清楚地模拟出了组件之间存在的声波反射、干涉现象.

图14 瞬态脉动压力云图Fig.14.Contours of pressure fluctuation.

Pascioni和Cattafesta[39]测量了位于r=1 m,θ=287.5°处的噪声谱,其中r以收起的缝翼前缘作为起点,θ 是以下游方向为起点逆时针为正.观测点噪声通过FW-H积分方法计算,声源面选择为壁面,声源信息在非定常流动模拟时采样,采样频率为10 kHz,共计8024个数据点.将采样数据分为7个窗,数据重叠率为50%,采用Hamming窗对每个窗进行傅里叶变换,最终将每个窗的变换结果进行平均得到最终的声压级谱,计算声压级参考压力值为 5×10-5Pa.由于实验几何模型展向长度为1 m,通过Kato公式对计算数据进行校正.从从图15的结果看,观测点获得的远场声压级谱也显示出窄带峰,这与 P1表面脉动压力谱规律一致.计算结果在1.5 -10.0 kHz范围内与实验值相比符合较好,尤其是HRILES幅值和频率都可以准确模拟,而IDDES的计算结果幅值偏大.

参考机体噪声会议基准算例要求,图16给出了位于 r=10C 处的声压级指向图.不同于图15实验结果的展长1 m,r=10C 处噪声计算要求展长统一为1 in,小于本文计算域展向长度(2 in),可以通过只统计z=0.5-1.5 in之间的壁面声源数据实现,对频段256 Hz 至10 kHz进行带宽滤波得到最终声压级,声压参考值为 5×10-5Pa.IDDES 模型的结果表现出更强的偶极子声源指向特性,第二第四象限的噪声强度更大.Choudhari和Lockard[40]给出了 θ=270°处声压级的参考范围为56.7-69.9 dB,平均值为61.5 dB.HRILES模拟得到的结果为61.5 dB,IDDES模拟得到的结果为66.6 dB,均在参考范围之内.

图17对比了分别采用缝翼、主翼、襟翼作为声源面计算得到的远场噪声结果.不同组件计算结果表征的是各自壁面的非定常压力脉动信息,其压力脉动来源于流动非定常特性,但同时包括对其他组件产生的声波的反射效应,可以作为研究声源强度分布的参考.整体上缝翼噪声在高升力体噪声中占主导地位,其偶极子辐射指向特性也更显著,其次为主翼噪声,指向性不明显.襟翼噪声占比最小,也呈偶极子辐射特性,噪声辐射在第一第三象限强度较大,声源主要来自尾缘的涡脱落.

图15 r=2.19C,θ=287.5° 观测点声压级功率谱Fig.15.Power spectra density of sound pressure level at r=2.19C,θ=287.5°.

图16 r=10C ,远场声压级指向图Fig.16.Directivity of SPL at r=10C.

图17 各部件远场(r=10C)声压级指向对比图 (a) IDDES; (b) HRILESFig.17.Directivity of components' SPL at r=10C ：(a) IDDES; (b) HRILES.

5 结 论

本文发展了高精度混合RANS/HILES方法.对亚临界态(Red=4.3×104)单圆柱绕流和30P30N高升力体基准算例进行了非定常模拟,结合声比拟方法得到了气动噪声结果,并将HRILES 模型与DDES,IDDES模型进行了对比.

从单圆柱模拟的流动结果看,HRILES模型能够模拟亚临界态尾迹区层流转捩现象,捕捉到背风面小尺度分离泡,因此对阻力系数的平均值和均方根值以及壁面压力分布的预测结果比SST-DDES模型更接近实验数据.远场噪声的功率谱密度分布的幅值与频率都与实验结果符合较好,证明了HRILES 模型具备预测复杂流动气动噪声问题的能力.

将HRILES模型应用于30P30N高升力体气动噪声问题,计算条件与机体噪声研讨会中规定的基准算例条件保持一致,并采用了JAXA提供的标准网格,同时采用IDDES模型模拟以对比模型差异的影响.从流动结果看,HRILES方法由于在分离区模型耗散更低,得到的缝翼凹腔中的涡结构更精细.对比剪切层各个站位的平均速度、展向涡量和湍动能分布,HRILES模型的结果与实验数据更加接近,IDDES 模型的结果存在显著的推迟转捩的现象,导致剪切层上、中区域展向涡量偏大,湍动能偏低.二者计算出的壁面非定常压力脉动的均方根值与实验数据相符,并得到了正确的频谱分布.从远场噪声结果看,噪声频谱分布特别是尖频噪声与实验值符合较好,这说明数值模拟捕捉到了缝翼凹腔自激励共振模态,但是IDDES模型的结果声压级幅值偏高.通过选取高升力体不同组件壁面作为声源面,对比了各自对远场噪声分布的贡献.结果都表明缝翼噪声贡献了大部分声能,缝翼和襟翼噪声都展现了显著的偶极子特性,主翼的噪声辐射指向特性不明显.

综上所述,本文通过算例应用与模型对比,确认了HRILES模型在典型气动声学问题中的适用性.HRILES作为广义RANS/LES模型的一种,能够在SST模型与HILES方法之间光滑过渡.对比HRILES和IDDES模型,二者在壁面区域都主要采用RANS模拟,且HRILES在大涡模拟的尺度判据计算上与IDDES 保持了一致,但后者在LES区域近似等效于Smagorinsky亚格子模型,这也是DES类方法的共性.HRILES方法具备了HILES方法对于含能尺度流动的高分辨率模拟能力,因此在混合层流动模拟结果优于传统的IDDES方法.同时在边界层采用SST模型降低隐式大涡模拟方法对壁面网格分布的要求,减轻计算资源需求,在中、高雷诺数流场模拟以及气动降噪应用研究中具有优势,未来将对模型在更高雷诺数的复杂外形应用中开展进一步研究.

感谢中山大学国家超级计算中心在计算资源方面提供的帮助.