2. 中国飞机强度研究所 航空声学与动强度航空科技重点实验室, 陕西 西安 710065
对于大多数声学问题, 流场计算和声场计算对离散格式的要求是大不相同的。特别是在需要关注于某些问题的声学特性时, 流场计算格式会引起声波的耗散和色散, 这对于声场模拟来说是不利的。一方面, CFD格式通常被设计为解决近场的问题, 因此扰动量在平均流中会衰减很快。另一方面, 声场模拟主要是个远场问题, 声源区域的大小相比于传播到观测点的距离而言是个小量。通常情况下, 声源区域外的压力扰动很小, 相比于流场脉动完全可以忽略不计。另外, 声波本质上是不稳定的运动, 其空间尺度、时间尺度都有较宽的分布范围, 声波波长一般大于流场特征长度。这些不同的特性(线性与非线性, 远场与近场, 时间依赖性与准静止, 大与小空间尺度)需要对计算方案作出不同的权衡。
由于计算流体力学和计算气动声学之间在很多应用上存在根本差异, 如果采用直接数值计算方法统一求解是很难处理两者的兼容性, 因此混合方法在处理噪声的预测问题上是一个比较明智的方式。目前, 对于大部分声学问题采用的都是混合计算气动声学方法, 声源生成一般采用CFD求解或统计声源的方法, 声波到远场的输运过程采用声类比输运或计算输运的方法。这种混合方式允许声源生成方法与声传播过程任意组合, 因此在实际应用中具备很好的灵活性。混合数值计算方法没有考虑声场与流场相互作用以及声波如何在流场中产生的问题, 但考虑到声波与流场能量的巨大差异, 大部分声学问题都可以完全忽略该问题。
针对气动噪声混合方法, 国内外主要集中在CFD耦合声比拟的方法。其中, Umesh Paliath等[1]采用高阶有限差分的大涡模拟(large eddy simulation, LES)模型耦合FWH(Ffowcs Williams & Hawkings)方程研究了发动机安装位置对气动噪声的影响。E Rebecca Busch等[2]采用了1.7亿网格量分别采用了CFD耦合声类比的混合方法和LES直接计算飞机降落状态下发动机反向旋转转子的气动噪声问题, 结果表明两者计算结果差异很小。A Fosso Pouangue[3]采用LES耦合FWH方程的混合方法研究了不同几何形状的发动机喷嘴对喷流噪声的影响。Y Colin等[4]研究了对转开放式转子安装位置对气动噪声的影响。L Wang等[5]采用IDDES(improved delayed detached eddy simulation)/FWH方法计算研究马赫数0.115和0.23情况下的起落架气动噪声, 计算与实验频谱曲线声压级峰值差异在2 dB以内。
对于CFD耦合声比拟的混合方法, 国内学者主要集中于工程应用, 对其理论研究相对较少。北京航天航空大学的孙晓峰等[6]基于时域无网格方法研究了平板尾迹辐射声场, 所得到的声场分布及其指向性均符合平板尾迹的偶极子声辐射特性。南京航空航天大学的龙双丽[7]通过DES/FW-H方法数值模拟及风洞试验研究了某型飞机主起落架噪声频谱特性、远场指向特性及其产生机制, 给出了各部件对总声压级的贡献, 这为后续起落架气动噪声的研究提供了参考依据。南京航空航天大学的招启军, 徐国华等[8]采用CFD/Kirchhoff耦合方法模拟了直升机旋翼高速脉冲气动噪声, 同时提出了一种谐波展开分析的方法, 简化Kirchhoff公式中被积函数的插值方法, 提高了插值效率和精度。
虽然CFD耦合声比拟的混合方法发展到今天已经相当成熟, 在工程中也得到广泛应用, 但是这种方法无法直观的反映各声源部件在传播过程中声波之间的相互作用。基于带源项的声传播方程数值求解的混合方法有效解决了该弊端, 并且能给出瞬时时刻的辐射声场。与CFD耦合声比拟的混合方法类似, 带源项的声传播方程数值求解的混合方法同样将计算域分为声源求解区域和声传播区域两部分。声源求解主要通过CFD方法提供非定常湍流速度或采用随机统计模型构造湍流速度。在求得声源项后, 通过声传播方程对声源进行传播。在对声源传播过程中, 为了满足声波能量级小、空间尺度分布广等特点, 数值离散往往需要采用低耗散、低色散的高阶格式。本文基于此, 研究了耦合RANS/LES模型与LEE(linearized euler equations)方程的气动噪声混合预测方法, 并利用该方法计算分析了M219空腔噪声问题。
1 线化欧拉方程声学问题常常具备非定常、能量级小的特点, 常规的CFD数值离散格式无法进行准确模拟, 其数值误差甚至可能会掩盖真实的声场。声场的另一特点是时间和空间尺度分布广, 声波波长往往大于流场特征尺度, 这使得声学问题的数值模拟相比于CFD而言, 其计算域更大。为了能准确模拟声波的这些特性, 若仅采用相容、稳定的高精度格式并不能满足计算要求, 数值格式必须满足低耗散、低色散、各向同性的要求。
文中对声源的传播计算采用线化欧拉方程, 其直角坐标下的具体形式如下:
|
(1) |
式中
|
式中，u′、u分别为流向速度的脉动量和平均量(其余变量定义类似), ut, vt, wt为湍流速度3个分量。
对于控制方程(1), 空间离散采用了9点5阶的DRP格式, 为了抑制数值振荡, 在空间离散格式的基础上添加了9点模板的人工黏性格式[9], 时间推进采用了高精度大时间步长的龙格库塔格式HALE-RK64[9], 远场边界条件采用了理想匹配层无反射边界条件[9]。
2 SST-SAS模型Menter综合考虑了k-ω湍流模型具备高精度模拟逆压梯度流动的特点, 以及k-ε湍流模型对初始参数不敏感的优势, 通过混合函数方法提出了剪应力输运模型(shear stress transport, SST)。SST模型的无量纲形式为:
|
(2) |
|
(3) |
模型中的参数详见参考文献[10]。
基于SST模型的SST-DES(detached-eddy simulation)方法是由Strelets等人提出的[10]。SST-DES模型中湍流尺度参数定义为lk-ω=k1/2/βkω, 其中k方程为:
|
(4) |
在SST-DES方法中, lk-ω由min(lk-ω, CDESΔ)表示, 其中Δ定义为网格单元的最大边长。常数CDES通过混合关系式获得:
|
(5) |
式中, CDESk-ω=0.78, CDESk-ε=0.61。
在靠近壁面边界层时, 由于lk-ω远小于网格单元尺度, 此时模型处于RANS的工作模式, 而当远离物面时, lk-ω增大到大于CDESΔ时, 模型开启LES的工作模式。模型方程中的生成项SP和耗散项SD可分别表示为: 
SST-DES方法通过比较湍流尺度与当地网格尺度相比较的方式进行湍流工作模式的切换, 相比于SA-DES模型相比, SST-DES模型对计算网格的依赖有所降低。
Menter为了改善网格依赖性这个缺点, 提出了SST-DDES(delayed detached-eddy simulation, DDES)模型, 将SST-DES模型中的滤波尺度重新进行了构造。在SST模型中存在着直接指示附面层的混合函数F1、F2, Menter等人利用该参数重新定义了湍流尺度的控制函数FDES, 其作用是保证亚格子模型在边界层内关闭。
|
(6) |
在SST-DDES模型中Fsst可以是F1或F2, 将k方程写成如下形式:
|
(7) |
当FSST≡0时, (7)式转化为最初的SST-DES模型。
两方程SST-DDES模型相比于SST-DES模型, 减小了网格的敏感性, 但其在LES区的滤波尺度仍然包含了当地网格尺度, 即在时间步长分辨率固定的前提下, 空间流动尺度分辨率完全取决于网格分布。
随后, Menter基于SST湍流模式提出了尺度适应模型, 与DES方法不同的是, 尺度适应模型[10] (scale adaptive simulation, SAS)对流场的分区不再依赖网格自身的尺度分布, 而是由所构造的冯·卡门长度、湍流尺度以及当地流场结构共同决定。其中ω方程变换为
|
(8) |
为保持附面层内仍保持SST湍流模式, Menter将上方程变换为原始ω方程和QSAS的组合:
|
(9) |
式中，Lvk=κ|U′/U″|, 
从湍流黏性系数μt表达式来看, 尺度适应模型在远离壁面处起到类似大涡模拟的作用, 此时滤波尺度仅由当地的流动结构决定。
3 气动噪声混合预测方法集成非定常CFD求解方法和线化欧拉方程, 搭建了气动噪声混合预测方法, 其具体框架如图 1所示。
|
| 图 1 混合预测方法框架 |
分为以下5个步骤:①在进行气动噪声计算之前, 针对CFD计算和CAA计算的特点分别准备2套网格, 分别为命名为CFD网格和CAA网格; ②基于CFD网格, 通过求解可压缩的NS方程获得非定常湍流脉动场信息; ③将CFD网格上的湍流脉动场信息传递到CAA网格上; ④根据湍流速度信息提供给线化欧拉方程的源项, 离散求解获得声场。
4 算例验证——M219空腔噪声预测M219空腔模型是美国军方用来研究武器舱空腔流动及气动噪声的验证模型之一。该模型已在英国的ARA风洞[11]中进行了一系列的实验测试, 其相应的实验数据常被用作研究武器舱空腔流动和气动噪声的基准数据。该模型长深宽比为5:1:1, 其中空腔长度L=0.508 m, 深D=0.101 6 m, 宽W=0.101 6 m。空腔底部分布了10个监测点, 具体几何尺寸如图 2所示。
|
| 图 2 几何示意图 |
模型的计算坐标如图 3所示, x轴方向为空腔长度方向, y轴方向为空腔宽度方向，z轴方向为空腔深度方向。风洞试验状态为:总压P0=99 612.06 Pa, 总温T0=305.06 K, 来流马赫数Ma=0.85, 单位尺度雷诺数Re=13.46×10+6, 采样频率f=6 000 Hz。文中所采用的计算模型尺寸以及计算状态与风洞实验均保持一致。对于空腔的计算策略为:首先进行定常流场计算, 然后将其作为初始流场进行非定常计算, 非定常计算采用作者所在的研究团队自主开发的混合RANS/LES方法程序(尺度自适应模型SST-SAS)进行计算, 程序的可靠性已在参考文献[12]中进行了验证说明。CFD计算时间步长为8.33×10-5 s。整个CFD计算域为空腔长度的5倍, 网格总数为2 100万(具体如图 5所示), 附面层第一层网格法向距离为5×10-6, 经计算y+≈1.0。
|
| 图 3 模型坐标系说明图 |
|
| 图 4 周期T内的瞬时涡量云图 |
|
| 图 5 M219空腔计算网格图 |
图 4给出了三维空腔一个周期T内瞬时涡量云图。
由于腔内外流速的差异, 气流经过空腔前缘时, 在空腔上方形成剪切层。剪切层在沿流向的发展过程中, 由于能量不足而产生不稳定运动, 进而发生脱落、分离。脱落涡在向空腔后壁发展的过程中, 在腔后缘再附着或与腔后壁相撞, 向空腔内流动输入能量, 从而在腔内形成多个小旋涡。从旋涡与空腔的相互作用来看, 脱落涡在空腔后缘的附着或与空腔后壁的碰撞, 产生强烈噪声, 形成一次声波; 由于受空腔内前后的压力差的作用, 分离流通过空腔底部向上游运动, 与腔前壁发生碰撞激发空腔前缘新的涡生成、发展与脱落, 而后与后壁相撞, 形成二次声波。在满足一定条件时, 腔内流动将形成自激振荡现象, 诱发腔内产生多峰值的强烈噪声。
基于上述针对三维空腔发声机理的计算分析, 下面分别采用RANS/LES+LEE方法和统计声源与声传播(stochastic noise generation and radiation, SNGR)方法开展了三维M219空腔模型的气动噪声预测。通过计算分析, 判断这2种类型的气动噪声混合预测方法对于空腔这种流激振荡类型的气动噪声问题的适用性。该模型声场计算采用无量纲时间步长Δt=0.000 5。流场计算的CFD网格和声场计算的CAA网格如图 5所示, CFD网格总数2 100万, CAA网格总数3 500万。图 6为CFD网格与CAA网格湍动能分布图。
|
| 图 6 空间截面的湍动能分布云图 |
对于空腔M219模型的湍流速度生成与噪声的传播, 文中分别采用了SNGR和RANS/LES+LEE 2种方式。图 7为采用SNGR方法求解得到的湍流速度分布云图。图 8为采用RANS/LES方法获得的湍流速度分布云图。从这2幅图中可以看出:无论哪种湍流速度生成模型, 其湍流速度主要集中于空腔后壁附近, 与湍动能分布范围基本一致, 并且其无量纲量级均在±0.15左右。
|
| 图 7 湍流速度分布(SNGR) |
|
| 图 8 湍流速度分布(RANS/LES) |
图 9为某一时刻的声传播云图, 图 9a)为SNGR计算结果, 图 9b)为RANS/LES+LEE计算结果。
|
| 图 9 声传播云图 |
从图中可以看出:无论哪种方法其计算的声压量级基本一致, 声源主要集中在空腔后壁面, 腔体后端产生的噪声不断向腔体左上方传播。
为了对比SNGR和RANS/LES+LEE方法计算空腔噪声的可适性, 在空腔底部等比分布了10个监测点, 从前壁到后壁监测点分别命名为Kulite20到Kulite29, 具体如图 3黄色点所标记。图 10到图 12分别为空腔底部K20、K25、K29这3个监测点的声压曲线。
|
| 图 10 Kulite20的声压曲线 |
|
| 图 11 ulite25的声压曲线 |
|
| 图 12 Kulite29的声压曲线 |
图 13至图 15分别为空腔底部10个监测点的声压级频谱特性曲线图。
|
| 图 13 Kulite20的声压级频谱曲线对比结果 |
|
| 图 14 Kulite25的声压级频谱曲线对比结果 |
|
| 图 15 Kulite29的声压级频谱曲线对比结果 |
图中实线分别为本文采用RANS/LES+LEE方法及采用SNGR方法计算的结果; 图中虚线为Rossiter的经验关系式f=(U∞/L)((m-a)/(M∞(1+(γ-1)/2M∞2)-1/2+(1/κ)))计算出的前4阶模态频率, 分别为159 Hz、371 Hz、582 Hz、794 Hz。通过计算结果与实验结果的对比, 从图中可以看出:(1)SNGR方法无法预测出不同模态的声压级峰值, 但总体的频谱曲线变化趋势与实验结果相当。(2)RANS/LES+LEE方法能准确预测出空腔自激振荡的固有模态, 各阶模态的声压级峰值比实验结果略大, 但总体的频谱曲线变化趋势与实验结果吻合的很好。
图 16为空腔底部中心剖面观测点处的声压均方根prms分布曲线图, 其中(prms=sqrt(∑(pi-p0)2)/n, i=1, …n)。
|
| 图 16 压力脉动均方根曲线图 |
通过对比分析可以看出:SNGR方法所计算的声压均方根最大, 与实验结果对比有所差距, RANS/LES+LEE方法相对于SNGR方法更加接近于实验值。从计算结果与实验结果的声压均方根曲线的整体趋势可看出:前壁到空腔中部, 这一区域的声压强度变化趋于平缓, 从中部到后壁这段区域值增长迅速。这主要是由于气流流经空腔前缘, 在空腔上方形成剪切层, 剪切层在空腔前半部分还处于相对稳定的状态, 跨过腔中部, 气流未深入腔内部触及腔底面, 故腔前prms部和中部声压变化较为平缓; 剪切层发生分离, 并与腔后壁碰撞, 产生强烈的压缩波, 导致后壁处声压强度明显。
5 结论通过上述的对比分析, 可以得出以下结论:
1) SNGR这类采用重构湍流速度模型的混合噪声预测方法不适合模拟空腔这类有着强烈自激振荡现象的气动噪声问题, 它更适合模拟类似于翼型后缘、缝翼远场处的宽频气动噪声问题;
2) RANS/LES+LEE方法由于湍流速度采用的是非定常求解手段, 为此对于空腔内部非定常流动、存在固有频率的噪声问题更适合。
| [1] | Paliath Umesh, Premasuthan Sachin. Large Eddy Simulation for Jet Installation Effects[R]. AIAA-2013-2137 |
| [2] | Busch E Rebecca. Aeroacoustics of a High-Fidelity CFD Calculation of a Counter-Roataing Open Rotor in Take-Off Conditions[R]. AIAA-2013-2202 |
| [3] | Pouangue A Fosso. Subsonic Jet Noise Simulations Using both Structured and Unstructured Grids[J]. AIAA Journal, 2015, 53(1): 1024-1035. |
| [4] | Colin Y, Wlassow F. Installation Effects on Contra-Rotating Open Rotor Noise at High-Speed[R]. AIAA-2014-2971 |
| [5] | Wang L. Detached-Eddy Simulation of Landing Gear Noise[R]. AIAA-2013-2069 |
| [6] |
洪志亮, 高鸽, 景晓东, 等. 一种预测平板尾迹噪声的时域无网格方法[J]. 航空学报, 2015, 36(11): 3501-3514.
Hong Zhiliang, Gao Ge, Jing Xiaodong, et al. A Grid-Less Time Domain Method for Plate Trailing Edge Noise Prediction[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(11): 3501-3514. (in Chinese) |
| [7] |
龙双丽, 聂宏. 飞机起落架气动噪声数值仿真与实验[J]. 航空学报, 2012, 33(6): 1002-1013.
Long Shuangli, Nie Hong. Simulation and Experiment on Aeroacoustic Noise Characteristics of Aircraft Landing Gear[J]. Acta Aeronautica et Astronautica Sinica, 2012, 33(6): 1002-1013. (in Chinese) |
| [8] |
招启军, 徐国华. 基于CFD/Kirchhoff方法的直升机旋翼高速脉冲噪声分析[J]. 计算物理, 2006, 23(2): 137-143.
Zhao Qijun, Xu Guohua. A Numerical Study of High-Speed Impulsive Noise of Helicopter Rotors with the CFD/Kirchhoff Method[J]. Chinese Journal of Computational Physics, 2006, 23(2): 137-143. (in Chinese) |
| [9] | Fang Q, Hu X D, Li D K. Absorbing Boundary Conditions for Nonlinear Euler and Navier-Stokes Equations Based on the Perfectly Matched Layer Technique[J]. Journal of Computational Physics, 2008, 227: 4398-4424. DOI:10.1016/j.jcp.2008.01.010 |
| [10] | Strelets M. Detached Eddy Simulation of Massively Separated Flows[R]. AIAA-2001-0879 |
| [11] | Henshaw M J de C. M219 Cavity Case, Verification and Validation Data for Computational Unsteady Aerodynamics[R], RT0-TR-26, 2000 |
| [12] |
余培汛, 白俊强. 尺度适应模型用于模拟射流流量对空腔压力脉动的抑制[J]. 声学学报, 2015, 40(1): 71-81.
Yu Peixun, Bai Junqiang. Suppression Effect of Jet Flow on Pulsating Pressure of Cavity Using Scale-Adaptive Simulation Model[J]. Acta Acoustic, 2015, 40(1): 71-81. (in Chinese) |
2. Laboratory of Aeronautical Acoustics and Dynamics, Aircraft Strength Research Institute of China, Xi'an 710065, China
