基于转捩模型及声比拟方法的高精度圆柱分离涡/涡致噪声模拟∗

王光学 王圣业 葛明明 邓小刚

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

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

(2017年12月17日收到；2018年8月10日收到修改稿)

1 引 言

圆柱尾迹涡及涡致噪声问题在航空航天、风工程等实际工程中非常具有代表性.如何准确预测涡街特性、涡街产声大小一直是计算流体力学(CFD)及计算气动声学(CAA)的热点.尤其是在亚临界雷诺数范围内(1×103—2×105),圆柱边界层为层流而尾迹已转变为湍流涡,对CFD中的湍流模型是个巨大的挑战[1].近年来,随着计算机技术的飞速发展,一些简单的学术问题通过采用大涡模拟(LES)方法已经能得到很好的结果[2−5].然而,对于如大型客机起落架、大型通风管道等雷诺数较高并且特征尺度较大的工程外形,边界层内最大湍流尺度相对于几何尺度也会变得很小,使计算花费巨增.因此,基于雷诺平均Navier-Stokes(RANS)方程的湍流模型仍然是工业应用CFD/CAA的中坚力量.

目前工程应用中,转捩主要依靠经验或半经验公式确定.其中,Menter和Langtry[6]提出的基于当地关联的γ-Reθ转捩模型,凭借其与现代CFD方法良好的兼容性,在工程领域得到了广泛的应用.然而,当转捩后出现大规模非定常分离时,γ-Reθ模型仍受到传统RANS方法的局限,在分离区会高估涡黏性而无法准确预测气动力.另一方面,传统的尺度模拟方法,如分离涡模拟(DES)[7]、尺度自适应模拟(SAS)等[8],基于全湍流模型,对于高雷诺数的大分离问题能很好地求解,而对于中等雷诺数(1×104—3×105)时边界层附近存在转捩的问题却无法准确预测.

2006年,Langtry等[9]便针对该问题指出:“众所周知RANS湍流模型在大规模分离区无法准确模拟；而混合RANS/LES方法,如DES或SAS,在准确捕捉分离流动的物理特性方面可能是必要的.”其后,诸多学者为实现该目的,在结合γ-Reθ模型和尺度求解方法方面做了许多工作.2011年,Sørensen等[10]提出了一种结合DES和γ-Reθ模型的分离涡模拟方法,在雷诺数10—1×106大范围内的圆柱绕流问题中均能得到与实验符合的平均阻力结果.2013年,You和Kwon[11]发展出了基于γ-Reθ模型的SAS方法,也能在圆柱绕流问题中很好地预测平均气动力.2014年,Qiao等[12]通过转捩平板算例重点关注了DES方法产生的模化应力损耗(MSD)问题,表明该问题会阻碍边界层转捩的发生,而采用延迟分离函数可使该现象得到改善.2017年,Hodara和Smith[1]又通过结合γ-Reθ模型, 将Sánchez-Rocha和Menon[13]的Hybrid LDKM方法(一种加权函数型混合RANS/LES方法)扩展至转捩分离流动的模拟中,在NACA 63-415机翼和圆柱绕流的气动力预测方面取得了良好的效果.然而,上述学者对模拟结果的讨论主要集中于宏观气动力,对速度梯度、压力脉动以及声压级等更为精细的流场变化关注较少.同时,梯度和脉动量的准确模拟更依赖于低耗散的高精度数值格式,上述工作也未能涉及.

本文采用CFD和声比拟相结合的方法,在近场采用Tran-DDES方法结合七阶WCNS-E8T7格式对非定常湍流流场进行高保真求解；在远场使用基于Ffowcs Williams-Hawkings(FW-H)方程的积分法模拟观测区噪声.算例考虑亚临界雷诺数(约4×104)的圆柱分离/噪声问题,及圆柱尾迹中放置翼型后的流动干扰问题.

2 湍流及噪声模型

2.1 SST-γ-Reθ转捩模型

γ-Reθ转捩模型通过经验关联函数控制边界层内间歇因子的生成,再通过间歇因子控制湍流模型中湍流的生成.其无量纲守恒形式的输运方程为:

其中源项及系数具体参见文献[6].该模型包含三个关键的经验关联函数:转捩动量厚度雷诺数Reθt、转捩区长度Flength和临界动量厚度雷诺数Reθc,并对使用的CFD软件平台的差异较为敏感.本文作者在前期工作中已依据T3系列低速平板实验对经验关联函数进行了标定,并在低速问题中开展了应用[14].

γ-Reθ转捩模型通常需要耦合k-ω剪切应力输运(SST)模型构成四方程模型进行求解.无量纲SST模型输运方程的守恒形式为:

其中源项及系数具体参见文献[6].γ-Reθ模型和SST模型的结合主要是通过间歇函数γ来修正k方程的生成项和破坏项:

其中F1,orig是原SST模型中的混合函数；γfi是γ-Reθ模型经过分离诱导转捩修正后的间歇函数；γD,min和F3参见文献[6],ω方程不做修正.

2.2 Tran-DDES方法

分离涡模拟方法由Spalart等[15]提出,将LES方法和SA模型结合,并通过比较当地网格尺度与壁面距离实现两种方法的自动切换.其后,Strelets[16]借鉴Spalart的思想,通过比较当地网格尺度与湍流长度尺度,将DES方法引入SST模型.然而,早期DES方法的限制器在复杂网格上处理RANS和LES的转换过程中过于生硬,会造成模化应力损耗等问题.2006年,Spalart等[17]借鉴Menter的SST模型构造思想,采用“延迟LES函数”改善了原版本的MSD问题,发展出了DDES方法.由于γ-Reθ模型是耦合SST模型求解的,而DDES方法仅是针对k方程耗散项进行的修正,因此转捩-分离涡模拟(Tran-DDES)方法是直接对(5)式的耗散项修正为:

这里DDES限制器为

转捩限制器和分离涡模拟限制器均作用于k方程耗散项上,两者之间可能会产生干扰,需要引起关注.要避免DDES判断在转捩区开启,保证在层流-转捩区中经验关联函数不被DDES修正干扰,本文采用了如下保护函数:

而fd是原DDES方法中的延迟函数,

这里,ν和νt分别为运动黏性系数和涡黏系数,κ为冯.卡门常数0.41；dw为壁面距离.

2.3 Ffowcs Williams-Hawkings声比拟方法

远场噪声采用FW-H方程对近场计算得到的非定常流场参数进行积分得到.声源面可采用固体边界面或可穿透声源面两种.本文采用固体边界面,并且假设声源面包含了所有声源信息时的情况下,时域求解的FW-H方程[18]为

其中声源面用f(xi,t)=0表示,r表示观测点距声源点的距离,Li和H分别代表载荷噪声(偶极子声源)和厚度噪声(单级子声源).tadv代表推迟时间,即观测点接收到声波的时刻.

由于计算网格量限制,计算模型的展向长度往往与实验模型存在差别.因此需要根据声源信息的相干性对结果进行校正.较为常用的校正公式为Kato公式[19]:

其中,(Spp(f))exp为校正后声压；(Spp(f))sim为数值积分得到的待校正声压；Lsim,LC和Lexp分别代表计算模型展向长度、流动展向相关长度和实验模型展向长度.

3 高精度数值格式

加权紧致非线性格式(WCNS)首先由Deng和Zhang[20]提出.之后,一系列WCNS格式被发展[21,22]并广泛用于流动问题中[23,24].本文采用的WCNS-E8T7格式是由Liu等[22]通过在经典的5阶显式加权插值模板基础上添加紧致项构造而来.该格式在保持较短的插值模板长度的基础上,实现了精度及色散、耗散特性的显著提高.同时,在遇到激波或数值不稳定情况时,通过设计加权系数使紧致项关闭,7阶紧致插值格式又可退回至5阶显式插值,从而提高格式在较复杂外形求解时的鲁棒性.

具体的,7阶加权紧致插值公式为

权ωk由下式给出,

这里,

其中β0,β1和β2是光滑因子,和5阶显式加权插值相同.β5被设计用来控制紧致项,具体参考文献[22].

Wang等[25]在涡输运、三角翼等典型涡流动算例中,对WCNS-E8T7格式进了验证.结果表明该格式能够在较粗糙的网格上达到甚至超过传统二阶格式在极密网格下对涡结构的分辨率,即计算效率更高.Seo等[2],Boudet等[3]以及Jiang等[5]在模拟圆柱或圆柱-翼型算例时均采用了高阶精度的空间离散格式,以更好地避免格式色散、耗散误差对流场/声场的影响.因此为更方便地对比湍流模型的差异,本文空间离散同样基于高阶精度格式.为本文研究的转捩后涡街结构的保持、涡致噪声的准确模拟奠定了基础.另外需要指出的是,在计算网格导数及雅克比时,本文采用了对称守恒网格导数方法(SCMM)[26],以保证WCNS格式在曲线网格求解时的数值精度.

4 计算结果

4.1 单圆柱转捩分离

单圆柱算例考虑ReD=4.3×104雷诺数情况,此时圆柱壁面附近的层流边界层开始分离,剪切层随后发生转捩.流场结果参考Szepessy和Bearman[27]在ReD=4.3×104,圆柱展长Lz=2.7D的实验结果和Seo等[2]在ReD=4.6×104,圆柱展长Lz=3.0D的不可压缩LES结果.

本算例采用7阶紧致格式WCNS-E8T7,同时添加准3阶格式MUSCL(κ= −1/3)[28]作为参考.湍流模型选择Tran-DDES并和传统SSTDDES进行对比.计算网格采用O型拓扑,包含180×180×30网格单元,与文献[2]中网格设置保持一致.时间推进采用双时间步方法,100个左右对流周期(T=D/U)后,拟周期的涡街流动建立.此时,统计平均过程开始并维持200T.

图1给出了120T后升阻力系数随时间的变化过程,对比图2中Seo等的LES结果,Tran-DDES方法得到升阻力系数波动幅值与LES更为接近.表1给出了气动力的统计结果,升力系数随卡门涡街脱落周期振荡,实验测得的斯特劳哈尔数为St=0.19,WCNS结合Tran-DDES与之符合很好,而结合传统DES则表现稍差.另一方面,观察MUSCL结合Tran-DDES结果,其明显偏离了实验值.这是由低阶格式较大的色散误差造成的,并且会对后续声学结果的处理造成影响.CL,rms和CD,rms分别为升、阻力系数的均方根值,反映拟周期振荡的振幅大小.MUSCL和WCNS结合Tran-DDES均能与实验符合,比SST-DDES表现得更好.

图1 单圆柱算例升阻力系数随时间变化Fig.1.Time variations of the drag and lift cofiicients for the single cylinder.

图2 Seo等的单圆柱LES结果[2]Fig.2.The LES results obtained by Seo et al.using LES[2].

本文中采用的网格较密,在平均流场特性方面,WCNS与MUSCL结果差异不大.下面主要对比湍流模型之间的差异.图3对比了WCNS-E8T7格式得到的表面压力分布,右图给出了Seo等的LES结果作为参考.Tran-DDES在背压区包括吸力峰部分均与实验值符合得很好,而SST-DDES的在背风区负压值则偏小.观察图4中平均流场的涡强分布,采用全湍流模型计算时,圆柱两侧的湍流剪切层抵抗失稳能力更强.即SST-DDES预测的回流区长度更长,在背风面附近流动的动能也更大,因此静压比Tran-DDES得到的要小.虽然在圆柱表面,湍流边界层的摩擦阻力大于层流状态,但湍流尾迹引起的压差阻力减小占主导,因而造成总阻力偏小.

表1 单圆柱算例气动力统计结果Table 1.Statistical results of aerodynamic cofiicients for the single cylinder.

图3 单圆柱表面压力分布,右图为Seo等的LES结果[2]Fig.3.Mean pressure cofiicient around the cylinder,right:results by Seo et al.using LES[2].

图4 单圆柱平均流场的涡强分布Fig.4.Distributions of mean vorticity intensity for the single cylinder.

然后分析远场观测点噪声结果,参考Jacob等[29]在ReD=4.8×104,圆柱展长Lz=30D的实验数据.采用FW-H声比拟方法模拟远场观测点的噪声大小,观测点位于圆柱垂直于来流的轴线上R=185D处.由圆柱绕流声场的辐射指向特性可知,垂直于来流的轴线为偶极子轴,并且声压级最大.由于实验模型的展长为Lz=30D,而计算模型的展长仅为Lz=3D.因此本文采用了Kato公式,根据声源信息的相干性对结果进行校正.图5给出了该观测点处功率谱密度计算结果与实验的对比.随着卡门涡街的周期性脱落,声压级在能谱图上随频率(St数)形成了主峰、二次峰和三次峰等.主峰为圆柱壁面交替的正负压力脉动产生的偶极子声源,Jacob实验测得的主峰处St数接近0.2,相较Seo实验偏大,这与两次实验的来流雷诺数及模型差异有关.另一方面,由于FW-H积分公式中时间导数通常仅采用二阶差分格式,使观测点处声压级谱上的St数均略低于气动力统计得到的St数(表1),也部分导致了WCNS结合SST-DDES的主频更接近实验值.然而,对比主峰处声压幅值,WCNS结合Tran-DDES的结果优于SST-DDES.MUSCL结合Tran-DDES得到的主频和幅值均偏离实验稍远,这与表1中St数趋势相符合.

图5 圆柱上方R=185D处声压功率谱密度结果对比Fig.5.Comparison of PSD of acoustic pressure at R=185D directly above the cylinder.

4.2 圆柱-翼型干扰流动

单圆柱算例中,主要研究了转捩模型对圆柱壁面附近非定常流场的影响.可以看到,是否考虑转捩对预测的涡街形态有很大影响.下面进一步研究圆柱绕流与下游翼型的相互干扰.如图6,圆柱的直径为0.01 m,基于圆柱直径的雷诺数和单圆柱算例近似,为ReD=4.8×104.下游翼型为NACA0012,弦长c=0.1 m.翼型与圆柱的距离为0.1 m.以翼型前缘为原点,远场噪声观测点Pfar的坐标为(0.05 m,1.85 m).

图6 圆柱-翼型实验设置及观测点位置Fig.6.Setup of the rod-airfoil experiment and location of measurements of interest.

本算例采用7阶紧致格式WCNS-E8T7,湍流模型选择Tran-DDES和SST-DDES进行对比.计算网格采用H-O型拓扑,圆柱表面仍为180个网格单元,翼型表面为320个网格单元.展向拉伸0.03 m,均布30个网格单元,总网格单元数约300万.文献[3]采用传统LES方法模拟,圆柱和翼型表面分别布置了200,350个网格单元,展向布置30个网格单元.但由于法向方向网格量稍高,当前网格的总单元量反而高于文献[3].计算采用双时间步方法,时间步长取1.5×10−6s,平均统计时长0.03 s.结果参考Jacob等[29]在里昂中央理工大学(Ecole Centrale de Lyon)的亚声速无回声风洞的实验结果,展向长度0.3 m.

首先分析气动结果. 图7展示了0.03 s时,Tran-DDES模型和SST-DDES模型结合高精度WCNS-E8T7得到的Q等值面分布.圆柱两端剪切层失稳形成卡门涡街,尾迹涡打到机翼并破裂.特别是涡街破裂后的较小尺度湍流结构也能被清晰地分辨.

以翼型顶点为原点,观察平均流场中四个典型剖面的速度分布(图8).x/c= −0.87和x/c=−0.255为圆柱尾迹中两个典型站位,平均速度在圆柱流向轴线处达到最小值.在x/c=−0.87站位,两种模型的结果相比实验值均存在巨大偏差.Agrawal和Sharma[30],Jiang等[5]学者采用LES方法预测的速度值在中线处也明显低于实验值.Agrawal对该站位的实验值提出了质疑,因为尾迹速度亏损随着远离圆柱而减弱,即在x/c=−0.87处的速度最小值应该小于x/c= −0.255处.在x/c= −0.255站位处,两种模型得到的速度型没有明显差异,并且速度最小值均比实验结果低.x/c=0.25位于翼型最大厚度处,圆柱尾迹流过翼型前缘到最大厚度处是加速的过程.观察该站位实验值,在y/c=0.3附近速度达到最大值.两种模型均预测出了该趋势,与实验值符合较好.x/c=2.0为翼型尾迹,两种方法也均能较好地预测出速度型,SST-DDES比Tran-DDES稍接近实验值,但差异并不大.

图7 圆柱-翼型干扰流动Q等值面分布,颜色由马赫数标识Fig.7.Iso-surface of the Q-criterion=0 for the rod-airfoil,colored by Mach number.

图8 圆柱-翼型平均流场中四个典型站位速度分布图Fig.8.Comparisons of mean velocity at four slice stations in the flow-field of rod-airfoil.

图9 圆柱-翼型算例圆柱附近平均流场的涡强分布Fig.9.Distributions of mean vorticity intensity for the single cylinder.

图10 圆柱-翼型流动中四个典型站位速度均方根值分布图Fig.10.Comparisons of root mean square value of velocity at four slice stations in the flow-field of rod-airfoil.

上述四个站位的平均速度分布,Tran-DDES与SST-DDES的差异并不明显.进一步观察两种方法得到的平均涡量分布(图9),并与单圆柱构型情况的数值结果(图4)进行对比.可以发现,全湍流模型无法模拟层流-转捩过程而导致剪切层失稳推迟的问题减弱甚至消除了.文献[22]指出翼型的加入也会反过来影响圆柱流场,包括圆柱表面气动力及涡街形态.由于翼型的影响,SST-DDES方法模拟的圆柱尾迹回流区被压缩,并且恰好更接近真实情况.这也是许多文献中[31,32]采用全湍流方法仍能在圆柱-翼型干扰流动中取得较好结果的原因.

然后,对比四个典型剖面的速度均方根值分布(图10).在x/c=−0.87站位,计算结果仍与实验值存在较大偏差.在x/c=−0.255站位,Tran-DDES和SST-DDES均与实验符合很好,主要差异在y/c=0处速度均方根值是单个峰值还是两个峰值.来流流过圆柱后,分离涡在圆柱轴线两侧交替产生,因此在均方根图上会出现两个峰值的现象.然而实验也未捕捉到两个峰值的现象,这与该站位远离圆柱更接近翼型,并且受翼型流场干扰有关.在翼型最大厚度处的x/c=0.25站位,计算值均与实验符合得很好,仅在近壁面处Tran-DDES预测的速度均方根值优于SST-DDES.在翼型尾迹的x/c=2.0站位,实验预测出了明显的两个峰值,并且由于圆柱和翼型存在高低差,均方根分布并不对称.Tran-DDES成功预测出了两个峰值与实验趋势符合.风洞中来流通过实验模型湍流度会迅速升高,并且小尺度涡的速度脉动也会反映在均方根图上.但对于湍流模型而言,湍流度无法直观反映在流场中,并且小尺度涡会被平均(过滤)掉,因此在远离轴线后脉动值会趋近于0.

最后对比声学模拟结果,远场噪声采用FW-H方法模拟,并且同样通过Kato公式,根据声源信息的相干性对结果进行校正.图11给出了采用不同模型所得到的Pfar点声场频谱形状,并与实验结果进行了对比.Tran-DDES和SST-DDES得到的主频能较好地符合卡门涡街脱落的主频,并在整体上的形状也与实验符合较好.但功率谱密度的量级存在一定的偏差,相比实验值均偏大,部分原因可能是展向修正造成的.数值计算所得频谱高频振荡较为严重,这是由于数值计算采样时间不足采样信号不稳定造成的.后处理时需要对随机信号进行多次平均才能得到较为稳定的信号.这就需要延长采样时间,并在进行快速傅里叶变换时取多次平均可以改善高频振荡.

图12给出了以壁面为积分面在以翼型中点为圆心R=185D圆周上所得到的指向性结果与实验的对比,图中横坐标OASPL(overall sound pressure level)表示总声压级.可以看出两种计算模型预测的结果差异不大,并且均与实验得到的指向性规律一致,即表现为偶极子指向性,且偶极子长轴垂直于来流方向.说明圆柱脱落涡对翼型表面的周期性撞击所产生的非定常载荷是产生噪声的主要原因.

图11 圆柱-翼型算例Pfar点功率谱密度结果对比Fig.11.Comparison of PSD of acoustic pressure at Pfarfor the rod-airfoil.

图12 圆柱-翼型算例声压级指向分布图Fig.12.Directivity curves of OASPL in the flow-field of rod-airfoil.

5 结 论

本文基于七阶WCNS-E8T7格式,结合DDES方法和FW-H声比拟方法,对亚临界雷诺数下单圆柱、圆柱-翼型的分离涡/涡致噪声问题进行了数值模拟.针对亚临界雷诺数下圆柱尾迹中的转捩问题,发展了基于γ-Reθ模型的高精度Tran-DDES方法,并与传统基于全湍流SST模型的高精度DDES方法进行了对比.

在单圆柱算例中,由于无法估计圆柱壁面附近的层流-转捩过程,传统DDES方法延迟了圆柱两侧剪切层的失稳,致使平均流场中回流区增大,压差阻力偏小.同时,涡街产生的频率变高,使声压频谱整体向高频偏移.而耦合了γ-Reθ转捩模型的DDES方法能预测出与实验符合的结果.

在圆柱-翼型算例中,由于翼型对圆柱附近流场的影响,传统DDES方法对圆柱两侧剪切层失稳的推迟减弱甚至消除了.因此在平均速度分布方面,SST-DDES方法与Tran-DDES方法并无明显区别.但在脉动量预测以及脉动产生的噪声预测方面,Tran-DDES方法与实验符合得更好.

由于数值计算中无法承受加大展长所增加的巨大花费,本文采用了Kato公式拟合真实展长下的噪声.但在声压主频及功率谱密度量级方面,并不能完全与实验符合,这也是下一步需要继续改进的方面.

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