基于HDCS-E8T7格式的气动噪声数值模拟方法研究进展

姜屹，邓小刚，毛枚良，刘化勇

(1．中国空气动力研究与发展中心 空气动力学国家重点实验室，四川 绵阳 621000; 2．中国空气动力研究与发展中心 计算空气动力学研究所，四川 绵阳 621000; 3．国防科学技术大学，湖南 长沙 410073)

基于HDCS-E8T7格式的气动噪声数值模拟方法研究进展

姜屹1，2，邓小刚3，毛枚良1，2，刘化勇1

(1．中国空气动力研究与发展中心 空气动力学国家重点实验室，四川 绵阳 621000; 2．中国空气动力研究与发展中心 计算空气动力学研究所，四川 绵阳 621000; 3．国防科学技术大学，湖南 长沙 410073)

·编者按·

为提高《空气动力学学报》学报质量，促进空气动力学学科发展，2014年3月编委会确定今年在《学报》上开设“高精度数值模拟方法”、“高超声速气动力/气动热”两个热点专栏。经过半年的征稿、审稿工作，共录用稿件20余篇，从本期开始分三期连续刊载。在专栏的组织落实中，得到了广大科研工作者的大力支持，在此表示诚挚感谢!

发展了气动噪声高精度数值模拟方法，空间导数离散采用了HDCS-E8T7格式及精度与之匹配的边界格式，时间离散方法采用了高精度多步龙格库塔方法和高精度隐式双时间步方法，发展了高精度对接边界算法，将方法推广至多块对接网格以满足解决复杂几何构型问题的需要，采用隐式大涡的概念处理可能出现的湍流问题。在此基础上，研究了几何守恒律对计算结果的影响，展示了复杂网格中高精度计算满足几何守恒律的重要性，完成了等熵涡、双圆柱散射、串列柱翼构型和喷嘴射流等典型噪声问题的求解，所得计算结果展示了所发展的模拟方法具有良好的预测精度和解决复杂构型气动噪声问题的潜力。

混合型线性耗散紧致格式;计算气动噪声;几何守恒律;隐式大涡模拟;复杂外形

0 引 言

气动噪声是人类生活和工业应用中常见的问题，会对人类健康和军事飞行器及潜艇的隐秘行动等多方面产生影响［1］，因此，抑制气动噪声是一个非常重要课题［2］。然而，由于实际问题中流动现象的复杂性，目前有效的噪声预测手段是把复杂噪声系统分解为不同的子系统，分别研究其中的噪声机制，如射流噪声和流动/物体干扰噪声。

近十多年以来，计算机技术迅猛发展，计算机资源得到了极大的丰富，为数值研究噪声问题提供了强大的硬件支撑。在射流噪声数值模拟研究方面，Bogey等［3-4］数值模拟了马赫数和雷诺数对射流噪声的影响，同时围绕准确预测射流噪声在数值技术和数值处理方面开展了大量的研究工作，比如边界条件和大涡模拟(LES)的亚格子模型。Morris研究了入口边界条件对噪声预测的影响［5］。Bodony和Lele则研究了亚格子模型对射流噪声计算的影响［6］。Bogey和Bailly在一系列文章中系统地研究了入口条件［7-8］和亚格子模型［9］对射流噪声计算的影响。最近，对射流噪声的研究已经把喷嘴考虑在计算域内，预测不同喷嘴的射流噪声，研究不同喷口形状对射流噪声的影响。其中包括双头喷嘴射流［10-12］，具有锯齿边缘的双头喷嘴射流［13］以及含有V型边缘的亚声速射流［14-19］。Engblom［20］等针对V型锯齿单股冷喷和热喷射流，采用雷诺平均Navier-Stokes（RANS）方法数值研究了锯齿内弯角度与射流流动的关系，尽管预测了远场噪声频谱的变化趋势，但是，他们的数值结果存在明显的不足［21］。事实上，对于射流噪声预测，RANS方法还不能胜任。Tide等在文献［21］中采用SST模型预测射流噪声，只能够得到不同位置噪声大小的变化趋势，而噪声绝对值与实验值差别则很大。采用LES方法可以很好地克服这个问题，比如Andersson［22-24］等采用LES方法预测射流噪声时，得到的数值解与实验值差异为±3d B。在流动与物体干扰噪声数值模拟研究方面，串列圆柱-翼型结构是研究这类噪声的基本模型，不同的计算方法被用于研究该干扰构型产生的气动噪声，包括RANS方法［25-26］、LES方法［25-27］和分离涡模拟（DES）方法［28-30］，而采用RANS方法预测串列圆柱-翼型干扰结构噪声没有得到可信赖的结果。Min等［31］通过实验和RANS方法研究了串列圆柱-翼型结构中，圆柱与翼型的相对位置对气动噪声的影响，他们的计算值与实验值差别非常大。Boudet等［27］通过对比发现采用LES计算串列圆柱-翼型干扰构型流动问题，不但可以得到比RANS更好的噪声频谱，而且可以得到更好的平均流场。出现这种现象的原因是RANS方法难于模拟这种湍流尾迹与壁面相互干扰的非定常现象，而这种现象是产生气动噪声的主要原因。

上面的研究实践表明，采用大涡模拟方法处理湍流效应非常有助于得到可靠的气动噪声预测结果。事实上，气动噪声问题与通常的动力学问题明显不同。首先，气动声学问题通常是非定常的波动问题，而且其脉动的频率非常宽，除了需要求解声源区域的流动特性外，还要探讨远场声波的指向性和声能量分布，这要求数值方法在全场都有一致的精度［32］；其次，声学脉动量与流动量相比非常小，通常比流动量低5～6个量级，这有可能使声学小量被数值误差掩盖［33］。为了正确和有效地计算气动噪声问题，气动噪声计算格式需要非常低的色散和耗散误差。由于高阶紧致有限差分格式的高分辨率和灵活性，它对于计算气动噪声而言是一种最具吸引力的方法［34］。

为了把HDCS-E8T7格式应用于复杂外形，除了满足几何守恒律，还需要多块结构网格技术，因为一般情况下很难为复杂外形生成单块网格。对于多块结构网格而言，首要问题是实现网格块之间信息高精度地传播。Thomas［61］在拼接网格上，发展了多维高精度插值方法，实现相邻网格块信息的传递，该方法在网格生成、计算效率等方面比对接网格具有明显的优势，但高精度多维插值方法比较复杂，涉及的插值模板点很多，如对于二维3阶插值，该文使用了12个点，特别是对于复杂情况，所希望得到的虚拟点可能并不在计算域中，实际计算精度并不能保证。Kim等［62］基于对接网格提出广义特征对接边界处理技术，能够消除几何奇性对计算精度的影响，实现网格块之间流动信息的高精度传递。邓小刚等基于流场信息沿特征线传播的原则，建立了对流信息沿特征方向传播的高精度对接边界算法［63］。在这种高精度对接边界算法的基础之上，HDCS-E8T7格式被成功地推广应用于复杂外形中数值模拟［64］。

本文主要介绍七阶HDCS-E8T7格式预测气动噪声的进展，重点将集中在格式应用于气动噪声预测的数值模拟方法、几何守恒律对数值计算的影响以及其对典型气动噪声问题的模拟情况。文中下一节将详细介绍数值方法，第二节展示了几何守恒律对数值计算的影响，典型的数值算例将在第三节中给出。

1 数值方法

1.1 控制方程

考虑曲线坐标系下的Navier-Stokes方程：

当θ＝1是Navier-Stokes方程，θ＝0是Euler方程，S为源项，它只在双圆柱散射算例测试中用到。

粘性应力项可以写为：

热传导项为：

其中γ为比热比，μ为粘性系数（可以通过Sutherland公式求解）。状态方程和能量方程为：

1.2 空间离散方法

七阶HDCS-E8T7格式被用于离散控制方程（1）。考虑无粘项的离散：

和它的半离散形式：

HDCS-E8T7内点格式为七点模板，因此在左侧j＝1，2，3和右侧j＝N，N-1，N-2需要三点边界和临近边界格式。临近边界点j＝3和j＝N-2采用六阶格式：

图1 HDCS-E8T7格式的修正波数Fig.1 Modified wavenumber of the HDCS-E8T7

临近边界点j＝2和j＝N-1采用四阶格式：

边界点j＝1和j＝N的四阶格式为：

边界格式中，j＝3／2和j＝N－1／2的五阶插值为：

同时j＝5／2和j＝N－3／2采用的五阶耗散插值为：

对于粘性项的计算，首先求解原始变量u、v、w、T的导数以得到应力张量和热传导通量。然后再一次使用七阶HDCS-E8T7格式得到粘性通量导数。

1.3 时间积分方法

采用了两种时间积分格式：显式格式和隐式格式。显式时间积分采用Hu提出的Runge-Kutta方法［66］。如果R表示右端项，控制方程为：

则Hu的五步Runge-Kutta方法通过式（14）从时间i0（第n步）积分至i0＋Δi（第n＋1步）。

如果计算网格拉伸严重，显式时间推进方法稳定性不足，隐式时间推进方法非常必要。本文的噪声预测体系中的隐式时间推进格式为基于牛顿迭代［67］的隐式双时间步方法［68］。二阶双时间步方法可以表示为：

其中M（U）＝∂R（U）／∂U，τ为虚拟时间。第一个子迭代时，p＝1，Up＝Un，而当p→∞时，Up→Un＋1。隐式时间推进算子中，空间离散采用二阶格式，而对右端R的离散采用HDCS-E8T7格式。

1.4 网格导数计算方法

根据邓小刚等在文献［40］中的研究，面积守恒律在均匀网格中自动得到满足，但是，在曲线网格中则有可能不能满足。如果面积守恒律没有得到满足，则在复杂曲线网格计算中可以会出现数值不稳定，甚至是计算崩溃现象。本小节，我们将根据满足面积守恒律的原则，讨论网格导数的计算方法。守恒网格导数（2）具有以下的离散形式：

为了使高精度有限差分格式满足面积守恒律，邓小刚等在文献［40］中提出了守恒网格导数算法CMM。CMM方法是通过使用与求解通量导数相同的格式计算守恒型网格导数（16）实现的，即δI＝δII。文献［40］中证明了CMM方法可以很容易应用于通量导数算子δI不分裂的高阶有限差分格式，但是对于将δI分裂为迎风算子δ＋I和δ－I的格式则很难。虽然守恒型网格导数（16）中的内层算子δIII对几何守恒律没有影响，邓小刚等在文献［40］中推荐关系式δIII＝δII。不久之后，他们从几何的角度解释了关系式δIII＝δII的意义［48］，并且提出了可以显著提高不规则网格中计算准度的对称型守恒网格导数算法SCMM［48］。为了满足面积守恒律，本文网格导数计算采用了SCMM方法［48］。半节点网格坐标等相关变量由高阶精度插值得到，其中内点δII和δIII的八阶显式插值为：

1.5 高精度对接边界处理方法

众所周知，针对复杂外形生成单块网格非常困难，一般需要多块网格技术，而且，网格线由于外形的原因通常会拐折。图2给出了这种网格的一个简单例子，从图中可以看出，网格线在奇点（图中以黑点表示）处突然发生拐折，显而易见，在这些点上左侧的网格导数不等于右侧的网格导数，这样就产生了奇性问题［62］。为了避免网格奇性，我们根据Kim等在文献［62］中的思想，把计算域沿着奇性线分解为两块。之后，在对接点和临近对接点为HDCS-E8T7格式构造特殊的对接边界格式。如图2所示，奇性点被分解为两个虚拟点。对于左侧和右侧网格分别计算导数，我们期望对接边界格式为不跨过对接面的单边格式，而且对接边界格式还需要满足几何守恒律以消除守恒误差即，δI＝δII＝δIII。为了满足这些要求，对接边界格式的构造分两步。

图2 沿奇性线分割计算域Fig.2 Decomposing a computational domain along the singular line

首先，构造δI、δII和δIII的插值方法。为了避免网格奇性，需要构造单侧插值。如图3（a）所示，奇性点j被分解为空间位置相同，但网格导数值不同的两个虚拟点jl和jr。在奇性点j的左右两侧，分别计算jl和jr的网格导数。jl和jr的单侧插值为：

第二步，采用半节点型δI构造单边格式计算网格导数。如图3（b）所示，我们把由奇性点分解得到的两个虚拟点重新合并在一起，从而对接点j的六阶半节点型δI（等于δII和δIII）可以写为：

图3 奇性点附近计算点示意图Fig.3 Sketch of computational grids near the singular point

对于网格导数计算，方程（20）将分别应用于奇性点的左侧和右侧。由于HDCS-E8T7格式采用七点模板，该对接边界格式将应用于［j－2，j＋2］。对于奇性点左右两侧的信息交换，将采用特征对接边界方法［63］。时间推进之后，对接点上的变量将通过平均左右两侧的值进行修正：

其中星号表示修正后的值。这种高精度对接边界处理方法可以保持HDCS-E8T7在复杂网格中的计算精度和稳定性［69］。

2 几何守恒律对数值计算的影响

为了展示面积守恒律的影响，我们设计了一种新的网格导数的算法，它通过加权守恒网格导数和传统网格导数得到计算网格导数值，如其中一个网格导数为：

图4 双圆柱散射计算网格示意图Fig.4 Computed domain and mesh of two-cylinder scattering case

2.1 双圆柱散射

该算例是第四届计算气动声学专题研讨会的标准算例，选用此算例的目的是展示几何守恒律在多块网格中对计算的影响。Euler方程（1）右端的声源项只在能量方程中出现，它的形式为：

为了逐渐地把源项引入到计算域中，在计算开始阶段，一个立方律的修形函数加入到源项（24）中。壁面边界上，法向速度分量V＝ηxu＋ηyv＋ηzw为0，切向速度分量U＝ξxu＋ξyv＋ξzw、压力和密度通过五阶外插得到。在远场边界处采用了远场无反射边界条件［70］和海绵层技术［71］。时间推进采用1．3小节中的五步Runge-Kutta方法，无量纲时间推进步长为0．002（一个周期为125步），总共推进20 000步，其中最后2 000步用于统计。图4为计算网格示意图。计算采用了三套网格，它们的网格量分别为469× 409（粗网格），721×685（中等网格）和841×775（细网格）。

在中等网格中，表1列出了HDCS-E8T7和HDCS-E8T7A格式的最大面积守恒误差。如图5所示，HDCS-E8T7A（σ＝1）格式的最大面积守恒误差出现在奇点附近。

表1 HDCS-E8T7和HDCS-E8T7A格式的最大Ix和Iy误差Table 1 Maximal errors of Ixand Iyfor the HDCS-E8T7 and HDCS-E8T7A

图5 σ=1时HDCS-E8T7A格式的|Ix+Iy|误差Fig.5 SCL error（|Ix+Iy|） of HDCS-E8T7A withσ=1

图6画出了HDCS-E8T7和HDCS-E8T7A计算得到的脉动压力均方根值等值线。图中显示当面积守恒误差足够小时（σ＝10-5），声场没有被明显污染。但是当σ从10-5变大到10-3时，声场在面积守恒误差最大的奇点附近变得越来越差。如图7所示，当σ大于10-2时，整个声场都被污染，而且最差的地方仍然是在奇点附近。这清楚地说明了满足面积守恒律对HDCS-E8T7格式的应用非常重要。图8给出了HDCS-E8T7格式使用三套网格计算得到的圆柱壁面脉动压力。当网格加密时，数值解逐渐接近理论解。采用细网格得到的数值解与理论解非常吻合，这展示了满足几何守恒律的HDCS-E8T7格式在复杂网格中的高分辨率。

图6 HDCS-E8T7和HDCS-E8T7A计算得到的脉动压力均方根值云图Fig.6 Mean-squared fluctuating pressure contours of HDCS-E8T7 and HDCS-E8T7A

图7 HDCS-E8T7A格式计算得到的脉动压力均方根云图Fig.7 Mean-squared fluctuating pressure contours of HDCS-E8T7A

图8 圆柱壁面脉动压力均方根分布Fig.8 RMS fluctuating pressure on the surface of cylinders

2.2 圆柱绕流

采用圆柱绕流算例，测试几何守恒律对HDCSE8T7格式模拟湍流问题的影响。计算条件是来流马赫数为0．2，基于圆柱直径D的雷诺数为ReD＝3900。本文的大涡模拟方法是基于七阶HDCSE8T7格式的一种隐式大涡模拟方法（HILES）［72］，虽然该方法与单调积分大涡模拟方法（MILES）［73］的思想是一致的，但是它具有以下两点特殊的性质［72，74］：（1）空间离散采用的是具备固有耗散的七阶HDCSE8T7格式；（2）HILES可以消除可能污染流场的面积守恒误差。

无量纲计算域为（90×60×π），网格大小为180 ×180×45。图9展示了二维网格，网格在圆柱壁面附近进行了加密，法向最小的网格间距为0．001。在外边界，采用了基于一维无粘近似的远场边界条件［70］和海绵层技术［71］。固壁面采用无滑移边界条件，等温壁，五阶精度外插实现压力法向梯度为0。展向采用周期边界条件。

图9 圆柱绕流的计算网格Fig.9 Grid system of the cylinder test

表2列出了HDCS-E8T7和HDCS-E8T7A的最大面积守恒误差，其中可以看出当σ＝1时，面积守恒误差将达到最大。如图10所示，HDCS-E8T7A（σ＝1）的最大面积守恒误差出现在奇点附近。

表2 HDCS-E8T7和HDCS-E8T7A的最大Ix和Iy误差Table.2 The maximal errors of Ixand Iyfor the HDCS-E8T7 and HDCS-E8T7A

均匀流场作为初场开始计算，时间推进采用1．3小节中的二阶隐式双时间步方法，流动在无量纲时间i＝100时达到全湍流状态。统计结果将由i＝100到i＝200的计算结果得到。计算中无量纲的时间推进步长为Δi＝0.01。图11画出了采用HDCS-E8T7和HDCS-E8T7A计算得到的时均流向速度云图，从图中可以看出，面积守恒误差会污染计算流场。为了定量显示面积守恒误差的影响，我们在图12中画出了x／D＝1.06处时均流向速度的垂直分布。图中清楚地显示，随着面积守恒误差的增大，非物理振荡越来越明显。这种现象同样在图12中S1（1.06，1.6）和S2（1.06，1.7）两点的速度分布中得到了体现。图13给出了瞬时涡量云图。图中显示了基于HDCSE8T7格式的HILES模拟方法计算得到了非常细致的结果，特别是在圆柱附近的尾迹区，但涡量可能被面积守恒误差严重污染，因此，满足几何守恒律对高精度湍流模拟同样非常重要。

图10 HDCS-E8T7A（σ=1）的面积守恒误差|Ix+Iy|Fig.10 SCL error（|Ix+Iy|） of HDCS-E8T7A withσ=1

图11 时均流向速度云图Fig.11 Mean streamwise velocity contours

3 测试算例

3.1 等熵涡传播

流动设定为二维，且（u，v，p，T）＝（1，0，1，1）。在初始均匀场中加入一个涡核位于（xc，yc）的等熵涡，它满足以下条件：

图12 展向平均的时均流向速度分布（x／D=1.06）Fig.12 Spanwise averaged mean streamwise velocity distributions（x／D=1.06）

图13 瞬时涡量云图Fig.13 Instantaneous vorticity contours

其中α＝β＝0.8为涡衰减参数，ε＝0.3表示涡强度，τ＝r／rc，r＝［（x－xc）2＋（y－yc）2］1／2，rc＝1.0为涡核半径。这里T＝p／ρ是温度，S＝p／ργ是熵。

计算采用均匀网格，且计算域为-10≤x≤10，-10≤y≤10。为了考察格式对波传播的分辨能力，本文设计不同的网格，使得涡核点数（PPVC），即涡核内部x和y方向上的网格点数为6～12。根据这个原则，得到了7套不同的网格，网格量分别为61× 61，71×71，81×81，91×91，101×101，111×111和121×121，对应的PPVC分别为6、7、8、9、10、11和12。

计算开始于i＝0和初始位置（xc＝0，yc＝0），一直持续至涡核运动至（4，0）处，时间推进采用1．3小节中的五步Runge-Kutta方法。此时等熵涡还未到达外边界，远场无反射边界［70］对计算结果的影响可以忽略。为了定量地考察PPVC与计算结果的关系，我们比较了不同PPVC值的压力数值误差的L1模、L2模和L∞模。图14展示了压力误差随PPVC的变化情况，从中可以看出HDCS-E8T7的计算值展现了非常好的网格收敛性。根据L1模、L2模和L∞模，可以得到相应的数值精度，表3列出了计算得到的数值精度，从中可以看出，HDCS-E8T7格式的设计精度得到了数值验证。

图14 PPVC对压力误差的影响Fig.14 Effect of PPVC on the errors of pressure

表3 HDCS-E8T7格式的数值精度Table 3 Accuracy of the HDCS-E8T7

图15给出了PPVC＝6时，得到的计算结果。从图中可以看出，由于HDCS-E8T7格式的高分辨率，等熵涡已经被很完美地捕捉，结合表3列出的压力误差值，我们可以看出七阶HDCS-E8T7格式对线性波传播问题具有很强的捕捉能力，这对预测气动噪声是非常有利的条件。

3.2 串列柱翼构型噪声

采用串列柱翼构型展示本文的方法预测复杂外形中气动噪声的能力。计算结果将与Jacob在文献［25］中给出的实验值进行比较，采用的计算条件与实验的相同，即：NACA0012翼型（弦长为c＝0．1m）位于圆柱下游一倍弦长处，来流速度为U∞＝72m／s，基于圆柱直径d＝0.1c的雷诺数为Red＝4．8×104。在这个构型中，圆柱后缘湍流尾迹与翼型前缘碰撞，期间伴随着涡的挤压，破裂等现象，这也是气动噪声产生的主要原因之一。HILES将被用于模拟串列柱翼构型中的湍流流动现象，而对于远场噪声的预测，采用的是FW-H方法［75］。

图15 PPVC=6时，等熵涡计算结果Fig.15 Numerical solutions of vortex convection test with PPVC=6

图16 计算网格及实验测量点位置示意图Fig.16 Sectional view of the computational domain and sketch of the locations for the measurements

图16画出了网格示意图，总网格量约为1600万，同时为了与实验数据对比，图16还展示了Jacob等在文献［25］中给出的实验测量点的位置示意图。为了进行大涡模拟计算，翼型和圆柱展向都延拓了0．03m，且展向均匀分布的网格点数为45。网格在圆柱和翼型壁面附近进行了加密处理，壁面距离为1．0×10-5（通过翼型弦长c无量纲得到的值）。计算中采用的边界条件为：远场采用无反射边界条件［70］和海绵层技术［71］（图16中橙色部分）；第1．5小节中的对接边界策略用于应对复杂多块结构网格；壁面采用无滑移边界条件；展向采用周期边界条件。时间推进采用1．3小节中的二阶隐式双时间步方法，无量纲的时间推进步长为0．001，对应的真实时间步长约为1．5×10-6s，计算初始流场为均匀来流，共推进30 000时间步，其中最后20 000步用于统计平均，由于展向采用了周期边界条件，因此，统计时均变量值及脉动量均方根值时，对展向进行了平均。

为了展示HILES计算得到的流场，图17画出了中间截面上瞬时展向涡量和涡量大小等值线云图。从中我们可以看出圆柱后缘的尾迹碰撞翼型前缘后，破裂为更小的涡向翼型后缘传播。考察流场中典型位置的速度统计量，图18对比了计算得到的x／c＝-0.255和x／c＝0.25处流向时均速度与相应的实验值。在x／c＝-0.255处，计算得到的中心线附近速度时均值比实验测量值低，这导致了x／c＝0.25处，壁面附近的速度时均值高于实验测量值，这种现象在Boudet等［27］的计算结果中同样可以观察到。图19给出了x／c＝-0.255和x／c＝0.25处流向速度均方根值分布与相应的实验值，从中可以看出本文的方法很好地预测了湍流脉动。圆柱后缘点的可分辨能谱同样证明了计算结果的可靠性，如图20所示。图中St＝fd／U∞为斯特鲁哈尔数，f为频率，U∞为来流速度。从中可以看出，数值计算达到了分辨惯性子区中湍动能的尺度要求，在此区间湍动能与St－5／3具有线性关系［76］。图20所展示的湍动能衰减斜率表明HILES捕捉到的湍动能能谱是可信赖的。图21画出了使用密度着色的第二不变量的等值面，第二不变量Q为：

其中Ωij＝（ui，j－uj，i）／2和Sij＝（ui，j＋uj，i）／2分别为速度旋度的反对称分量和对称分量。由于HDCSE8T7格式的高分辨率，我们可以看到许多涡结构细节。图22给出了（x＝0．68m，y＝1．74m）处噪声功率谱密度。计算得到的声压峰值及其频率与实验值吻合得很好。但是，对于高频部分，本文的计算值与实验值相比偏高，这与Boudet等［27］的计算结果类似。Boudet等在文献［27］中对这种现象给出的一个解释是：FW-H积分面外的体积分没有用于计算气动噪声。不过，整体上还是很好地预测了串列柱翼构型辐射的宽频噪声。

图17 流场中涡量等值线云图Fig.17 Contours of instantaneous vorticity

图18 展向平均的流向时均速度分布Fig.18 Spanwise averaged mean streamwise velocity distributions

图19 展向平均的流向脉动速度均方根值分布Fig.19 Spanwise averaged RMS value of fluctuating velocity

图20 圆柱后缘点的可分辨能谱Fig.20 Resolved energy spectrum at a location behind the cylinder

图21 速度张量的第二不变量等值面Fig.21 Iso-surfaces of the second invariant of velocity gradient tensor for instantaneous vortex structure

图22 远场噪声功率谱密度图Fig.22 Power spectral density of the pressure perturbations measured in the farfield

3.3 喷嘴射流噪声

本文考虑的喷嘴与Andersson在文献［22］中的相同。喷嘴直径Dj为50mm，基于喷嘴直径和喷嘴出口速度的雷诺数ReD为5．0×104，喷嘴出口马赫数为Ma＝0．75。计算的射流出口静温Tj等于背景温度T∞，射流入口条件与Tide在文献［21］中的条件相同。

图23给出了三维模拟的计算域平面示意图。为了减小计算域出口和远场反射对流动模拟的影响，出口和远场都加了海绵吸波层［71］。计算域在射流出口向前方向上延伸了20倍喷嘴直径，射流入口处计算域直径为10倍喷嘴直径，出口为20倍喷嘴直径，海绵层的长度为40倍喷嘴直径。计算域采用大约具有1.8×108个网格点的多块结构网格离散。图24展示了计算网格系统。

图23 计算域示意图Fig.23 Sectional view of the computational domain

在射流入口边界处，给定了静压和总压。在射流出口边界处，即在海绵层的末端，给定了静压。对于外边界，采用了远场无反射边界条件［70］和海绵层技术［71］。固壁面采用无滑移边界条件，法向压力梯度为0和绝热壁。

图24 喷嘴射流网格系统Fig.24 Mesh system for the nozzle

计算开始时，先采用定常算法直到计算域里的整体质量流量率变化小于射流出口质量流量率的1%。在此基础上开始非定常计算，时间推进采用1．3小节中的二阶隐式双时间步方法，时间步长约为Δi＝4μs，共推进30 000步，最后20 000步用于流动统计平均。然后采用FW-H方法［75］计算得到远场观测点的脉动声压。图25给出了这些观测点的位置。气动噪声通过FW-H积分面上的压力计算得到，这个积分面取得足够大以确保它包括了所有的声源，即这个算例中的剪切层。

图25 远场噪声观察点位置Fig.25 Locations of the far-field receivers

图26画出了瞬时轴向速度云图。图27比较了本文计算得到的中心线上的轴向平均速度与Andersson等在文献［22］中报告的计算值和实验值。可以看出本文计算得到的速度核心区长度和核心区后的速度衰减更好。图28画出了采用射流喷嘴出口速度无量纲化的流向（urms）和法向（vrms）脉动速度均方根值。与实验值相比，本文的计算很好地捕捉了脉动速度urms和vrms的峰值位置。但是在速度核心区x／Dj＜6的位置，计算得到了强度非常低的湍流。Andersson等在文献［22］中的计算结果同样观察到了这种现象。图29比较了计算得到的和实验测量得到的观察点上的声压级，从中可以看出本文的计算值与实验值吻合地很好。

图26 瞬时流向速度云图Fig.26 Contours of instantaneous axial velocity

图27 中心线上轴向速度分布Fig.27 Distributions of axial velocity along the centerline

图28 urms和vrms沿中心线上的变化情况Fig.28 Variation of urmsand vrmsalong the centerline

图29 远场观察点处的声压级Fig.29 Sound pressure levels at observer locations

4 结束语

本文介绍了基于HDCS-E8T7格式的高精度数值模拟方法及其在气动噪声模拟方面的研究进展，得到以下结论：

（1）所发展的算法具有预期的高阶精度；

（2）在复杂网格中，几何守恒律误差对计算结果具有显著影响，消除它对模拟气动噪声问题尤其重要；

（3）采用隐式大涡概念处理湍流，是高精度模拟湍流噪声的有效途径之一。

从模拟的对象和对问题的研究深度审视，本文对气动噪声模拟研究才刚刚起步。随着工作向求解实际工程问题迈进，将面临着诸多挑战，其中最基本的挑战有：

（1）网格生成的困难。对接网格可能难于应对工程实际问题中的复杂几何特性，需要对复杂外形适应能力更强的网格生成方法，比如拼接网格和重叠网格等，以及相应的分区边界算法。

（2）湍流模拟问题。采用隐式大涡概念可以有效地处理湍流，但求解效率还不足，目前还难以在工程实际问题中应用，需要发展计算效率更高的湍流求解方法。

（3）流固耦合的气动噪声问题。目前求解都是流动非定常的气动噪声问题，对于飞行器几何构型变化（包括位置变化和变形）引起的动态非定常噪声问题尚未涉及。

在下一步工作中，我们将围绕这些挑战，继续完善算法，在解决工程气动噪声问题中不断进步。

［1］WANG M，FREUND J B，LELE S K．Computational prediction of flow-generated sound［J］．Annu.Rev.Fluid Mech．，2006，38：483-512．

［2］JORDAN P，COLONIUS T．Wave packets and turbulent jet noise［J］．Annu.Rev.Fluid Mech．，2013，45：175-195．

［3］BOGEY C，BAILLY C．Investigation of subsonic jet noise using LES：Mach and Reynolds number effects［R］．AIAA 2004-3023．

［4］BOGEY C，BAILLY C．Direct computation of the sound radiated by a high-Reynolds number，subsonic round jet［R］．The CEAS Workshop From CFD to CAA，Confederation of European Aerospace Societies Paper 2，2002．

［5］MORRIS P J，LONG L N，SCHEIDEGGER T E，et al．Simulations of supersonic jet noise［J］．International Journal of Aeroacoustics，2002，1（1）：17-41．

［6］BODONY D J，LELE S K．Large eddy simulation of turbulent jets and progress towards a subgrid scale turbulence model［C］．Proceedings of International Workshop on LES for Acoustics，DLR Göttingen，Göttingen，Germany，2002：1-12．

［7］BOGEY C，BAILLY C．LES of a high Reynolds，high subsonic jet：effects of the inflow conditions on flow and noise［R］．AIAA 2003-3170．

［8］BOGEY C，BAILLY C．Effects of inflow conditions and forcing on subsonic jet flows and noise［J］．AIAA J.，2005，43（5）：1000-1007．

［9］BOGEY C，BAILLY C．LES of ahigh Reynolds，high subsonic jet：effects of the subgrid modellings on flow and noise［R］．AIAA 2003-3557．

［10］VUILLOT F，LUPOGLAZOFF N，RAHIER G．Double-stream nozzles flow and noise computations and comparisons to experiments［R］．AIAA 2008-0009．

［11］FAYARD B，RAHIER G，VUILLOT F，et al．Flow field analysis for double stream nozzle：application to jet noise［R］．AIAA 2008-2983．

［12］FAYARD B，RAHIER G，VUILLOT F．Modal analysis of jet flow from a coaxial nozzle with central plug［R］．AIAA 2009-3355．

［13］VISWANATHAN K，SHUR M，SPALART P，et al．Flow and noise predictions for single and dual-stream beveled nozzles［J］．AIAA J.，2008，46（3）：601-626．

［14］SHUR M L，SPALART P R，STRELETS M K．Noise prediction for increasingly complex jets．Part II：Applications［J］．International Journal of Aeroacoustics，2005，4（3-4）：247-266．

［15］SHUR M L，SPALART P R，STRELETS M K，et al．Further steps in LES-based noise prediction for complex jets［R］．AIAA 2006-485．

［16］XIA H，PAUL G．Tucker and simon eastwood，towards jet flow LES of conceptual nozzles for acoustic predictions［C］．46th AIAA Aerospace Sciences Meeting and Exhibit Reno，Nevada，2008．

［17］XIA H，TUCKER P G，EASTWOOD S．Large-eddy simulations of chevron jet flows with noise prediction［J］．International Journal of Heat and Fluid Flow，2009，30：1067-1079．

［18］UZUN A，HUSSAINI M Y．Simulation of noise generation in near-nozzle region of a chevron nozzle jet［J］．AIAA J.，2009，47（8）：1793-1810．

［19］UZUN A，HUSSAINI M Y．High-fidelity numerical simulation of a chevron nozzle jet flow［R］．AIAA 2009-3194．

［20］ENGBLOM W A，KHAVARAN A，BRIDGES J．Numerical prediction of chevron nozzle noise reduction using WIND-MGBK methodology［R］．AIAA 2004-2979．

［21］TIDE P S，BABU V．A numerical predictions of noise due to subsonic jets from nozzles with and without chevrons［J］．Applied Acoustics，2009，70：321-332．

［22］ANDERSSON N，ERIKSSON L E，DAVIDSON L．A study of Mach 0．75 jets and their radiated sound using large eddy simulation［R］．AIAA 2004-3024．

［23］ANDERSSON N，ERIKSSON L E，DAVIDSON L．Investigation of an isothermal Mach 0．75 jet and its radiated sound using large eddy simulation and Kirchh of f surface integration［J］．Int.J.Heat Fluid Flow，2005，26：393-410．

［24］ANDERSSON N，ERIKSSON L E，DAVIDSON L．Large eddy simulation of a Mach 0．75 jet［R］．AIAA 2003-3312．

［25］JACOB M C，BOUDET J，CASALINO D，et al．A rod-airfoil experiment as benchmark for broadband noise modeling［J］．J.Theoret.Comput.Fluid.Dyn.，2005，19（3）：171-96．

［26］CASALINO D，JACOB M C，ROGER M．Prediction of rod airfoil interaction noise using the FWH analogy［R］．AIAA 2002-2543．

［27］BOUDET J，GROSJEAN N，JACOB M C．Wake-airfoil interaction as broadband noise source：a large-eddy simulation study［J］．Int.J.Aeroacoustic，2005，4（1）：93-116．

［28］GEROLYMOSG A，VALLET I．Influence of temporal integration and spatial discretization on hybrid RSM-VLES computations［R］．AIAA 2007-4094．

［29］CARANI M，DAI Y，CARANI D．Acoustic investigation of rod airfoil configuration with DE Sand FWH［R］．AIAA 2007-4016．

［30］CRESCHNER B，THIELE F，CASALINO D，et al．Influence of turbulence modeling on the broadband noise simulation for complex flows［R］．AIAA 2004-2926．

［31］JIANG M，LI X D，ZHOU J J．Experimental and numerical investigation on sound generation from airfoil-flow interaction［J］．Appl.Math.Mech.-Engl.Ed．，2011，32（6）：765-776．

［32］TAM C K W．Computational aeroacoustics：issues and methods［J］．AIAA J.，1995，33（10）：1788-1796．

［33］TAM C K W．Recent advances in computational aeroacoustics［J］．Fluid Dynamics Research，2006，38：591-615．

［34］COLONIUS T，LELE S K．Computational aeroacoustics：progress on nonlinear problems of sound generation［J］．Progress in Aerospace Sciences，2004，40：345-416．

［35］WANG Z J．High-order methods for the Euler and Navier-Stokes equations on unstructured grids［J］．Progress in AerospaceSciences，2007，43：1-41．

［36］EKATERINARIS J A．High-order accurate，low numerical diffusion methods for aerodynamics［J］．Progress in Aerospace Sciences，2005，41（3-4）：192-300．

［37］NONOMURA T，IIZUKA N，FUJII K．Freestream and vortex preservation properties of high-order WENO and WCNS on curvilinear grids［J］．Computers and Fluids，2010，39：197-214．

［38］THOMAS P D，LOMBARD C K．Geometric conservation law and its application to flow computations on moving grids［J］．AIAA J.，1979，17（10）：1030-1037．

［39］PULLIAM T H，STEGER J L．On implicit finite-difference simulations of three-dimensional flow［R］．AIAA 78-10，1978．

［40］DENG X G，MAO M L，TU G H，et al．Geometric conservation law and applications to high-order finite difference schemes with stationary grids［J］．J.Comput.Phys，2011，230：1100-1115．

［41］VISBAL R M，GAITONDE D V．On the use of higher-order finite-difference schemes on curvilinear and deforming meshes［J］．J.Comput.Phys，2002，181：155-185．

［42］DENG X，ZHANG H．Developing high-order weighted compact nonlinear schemes［J］．J.Comput.Phys，2000，165：22-44．

［43］DENG X G，MAEKAWA H．Compact high-order accurate nonlinear schemes［J］．J.Comput.Phys，1997，130：77．

［44］DENG X G，MAO M L，JIANG Y，et al．New high-order hybrid cell-edge and cell-node weighted compact nonlinear schemes［R］．AIAA 2011-3857．

［45］JIANG G，SHU C．Efficient implementation of weighted ENO［J］．J.Comput.Phys．，1996，181：202-228．

［46］DENG X，MAEKAWA H，SHEN Q．A class of high order dissipative compact schemes［R］．AIAA 96-1972．

［47］TAKU NONOMURA，DAIKI TERAKADO，YOSHIAKI A B E，et al．A new technique for finite difference WENO with geometric conservation law［R］．AIAA 2013-2569．

［48］DENG X G，MIN YAOBING，MAO MEILIANG，et al．Further studies on geometric conservation law and applications to high-order finite difference schemes with stationary grids［J］．J.Comput.Phys.，2013，239：90-111．

［49］PALIATH U，SHEN H，AVANCHA R，et al．Large eddy simulation for jets from chevron and dual flow nozzles［R］．AIAA 2011-2881．

［50］DAUDE F，BERLAND J，EMMERT T，et al．A high-order finite-difference algorithm for direct computation of aerodynamic sound［J］．Computers&Fluids，2012，61：46-63．

［51］LELE S K．Compact finite difference schemes with spectral-like resolution［J］．J.Comput.Phys．，1992，103：16-42．

［52］RIZZETTA D P，VISBAL M R，MORGAN P E．A high-order compact finite-difference scheme for large-eddy simulation of active flow control［J］．Progress in Aerospace Sciences，2008，44：397-426．

［53］MAO M L，DENG X G．Boundary schemes and asymptotic stability for high-order dissipative compact schemes［J］．ACTA Aerodynamica Sinica，2000，18（2）：165-171．（in Chinese）

毛枚良，邓小刚．高阶精度线性耗散紧致格式的渐近稳定性［J］．空气动力学学报，2000，18（2）：165-171．

［54］MAO M L，DENG X G，LI S，Spectrum characteristic of dissipative compact schemes and application to couette flow［J］．Chinese Journal of Computational Physics，2009，26（3）：371-377．（in Chinese）

毛枚良，邓小刚，李松．耗散紧致格式的频谱特性研究与应用［J］．计算物理，2009，26（3）：371-377．

［55］MAO M L，JIANG Y，DENG X G．Study of LDDRK schemes for DCS5 scheme［J］．Chinese Journal of Computational Physics，2010，27（2）：159-167．（in Chinese）

毛枚良，姜屹，邓小刚．基于DCS5格式的LDDRK算法［J］．计算物理，2010，27（2）：159-167．

［56］JIANG Y，MAO M L，DENG X G．Application of DCS5 scheme in CAA［J］．ACTA Aerodynamica Sinica，2012，30（4）：431-436．（in Chinese）

姜屹，毛枚良，邓小刚．DCS5在计算气动声学中的应用研究［J］．空气动力学学报，2012，30（4）：431-436．

［57］DENG X G，JIANG Y，MAO M L，et al．Developing hybrid cell-edge and cell-node dissipative compact scheme for complex geometry flows［J］．Sci.China.Tech.Sci.，2013，56：2361-2369．

［58］DENG X G，JIANG Y，MAO M L，et al．A family of hybrid cell-edge and cell-node dissipative compact schemes satisfying geometric conservation law［J］．submitted toComputers& Fluids．

［59］JIANG Y，MAO M L，DENG X G，et al．Numerical prediction of jet noise from nozzle using seventh-order dissipative compact scheme satisfying geometric conservation law［J］．Applied Mechanics and Materials，2014，574：259-270．

［60］MAO M L，JIANG Y，DENG X G，et al．Noise prediction in subsonic flow using seventh-order dissipative compact scheme on curvilinear mesh［J］．Submitted toAdvancesin Applied Mathematics and Mechanics．

［61］THOMAS L G，et al．Multi-size-mesh，multi-time-step algorithm for noise computation around an airfoil in curvilinear meshes［R］．AIAA 2007-3504．

［62］KIM J W，LEE D J．Characteristic interface conditions for multiblock high-order computation on singular structured grid［J］．AIAA Journal，2003，41（2）：2341-2348．

［63］DENG X G，MAO M L，TU G H，et al．Extending weighted compact nonlinear schemes to complex grids with characteristicbased interface conditions［J］．AIAA J.，2010，48（12）：2840-2851．

［64］JIANG Y，MAO M L，DENG X G，et al．Extending seventhorder dissipative compact scheme satisfying geometric conservation law to large eddy simulation on curvilinear grids［J］．Submitted toAdvances in Applied Mathematics and Mechanics．

［65］TAM C K W，WEBB J C．Dispersion-relation-preserving finite difference schemes for computational acoustics［J］．Journal of Computational Physics，1993：107：262-281．

［66］HU F Q，HUSSAINI M Y，MANTHEY J L．Low-dissipation and low-dispersion Runge-Kutta schemes for computational acoustics［J］．J.Comput.Phys.，1996，124：177-191．

［67］JOHN M HSU，ANTONY JAMESON．An implicit-explicit hybrid scheme for calculating complex unsteady flows［R］．AIAA 2002-0714．

［68］GORDNIER R E，VISBAL M R．Numerical simulation of deltawing roll［R］．AIAA 93-0554．

［69］JIANG Y，MAO M L，DENG X G，et al．Extending seventhorder hybrid cell-edge and cell-node dissipative compact scheme to complex grids［C］．The 4th Asian Symposium on Computational Heat Transfer and Fluid Flow，Hong Kong，2013．

［70］POINSOT T，LELE S K．Boundary conditions for direct simulations of compressible viscous flows［J］．J.Comput.Phys.，1992，101：104-129．

［71］DANIEL J Bodony．Analysis of sponge zones for computational fluid mechanics［J］．J.Comput.Phys.，2006，212：681-702．

［72］JIANG Y，MAO M L，DENG X G，et al．Large eddy simulation on curvilinear meshes using seventh-order dissipative com-pact scheme［J］．Computers and Fluids（in press）．doi:10.1016/j.compfluid．2014．08．003．

［73］BORIS J P，GRINSTEIN F F，ORAN E S，et al．New insights into large eddy simulation［J］．Fluid Dyn.Res.，1992，10：199．

［74］JIANG Y，MAO M L，DENG X G，et al．Effect of surface conservation law on large eddy simulation based on seventh-order dissipative compact scheme［J］．Applied Mechanics and Mate-rials，2013，419：30-37．

［75］LYRINTZIS A S．Surface integral methods in computational aeroacoustics-from the（CFD）near-field to the（acoustic）farfield［J］．Int.J.Aeroacoust，2003，2（2）：95-128．

［76］XU C Y，CHEN L W，LU X Y．Large-eddy simulation of the compressible flow past a wavy cylinder［J］．J.Fluid Mech.，2010，665：238-273．

Progress of aeroacoustic simulation method based on HDCS-E8T7 scheme

JIANG Yi1，2，DENG Xiaogang3，MAO Meiliang1，2，LIU Huayong1
（1.State Key Laboratory of Aerodynamics，China Aerodynamics Research and Development Center，Mianyang 621000，China；2.Institute of Computational Aerodynamics，China Aerodynamics Research and Development Center，Mianyang 621000，China；3.National University of Defense Technology，Changsha 410073，China）

Noise prediction using high-order numerical method is one of the most interested research topics in CFD．HDCS-E8T7 scheme is a seventh-order hybrid linear compact dissipative scheme．This scheme overcomes the disadvantage that many high-order finite difference schemes may have when geometric conservation law（GCL）is not satisfied，the numerical instability phenomenon may be avoided when the HDCSE8T7 is applied to solve aeroacoustic problem in complex geometry．In the present high-order numerical noise prediction method，the spatial discretization adopts HDCS-E8T7 scheme and its boundary schemes have suitable accuracy comparing with the interior scheme．High-order multi-step Runge-Kutta and dual stepping scheme are employed for time integration，high-order interface boundary scheme is developed to extend the present method to multi-block point matched grid which meets the need of solving problem in complex geometry and the turbulence is handled by the concept of implicit large eddy simulation．Based on this method，the effect of GCL on the numerical solutions is studied．The importance of satisfying the GCL on complex grid for high-order simulation is illuminated．Some typical noise cases，such as vortex convection，scattering of acoustic waves by multiple cylinder，sound radiated by a rod-airfoil configuration and jet noise from nozzle，are investigated．The solutions of these tests show that the presented method has high prediction accuracy and potential application for handling aeroacoustic problem on complex geometry．

hybrid linear compact dissipative scheme；computational aeroacoustic；geometric conservation law；implicit large eddy simulation；complex geometry．

V211.3

Adoi:10.7638／kqdlxxb-2014.0103

0258-1825（2014）05-0559-16

2014-09-05;

2014-09-16

国家自然科学基金(11072259;11202226)

姜屹(1985-)，男，江西，助理研究员，研究方向：高阶精度数值计算方法研究与应用．E-mail：yijiang@mail．ustc．edu．cn

姜屹，邓小刚，毛枚良，等．基于HDCS-E8T7格式的气动噪声数值模拟方法研究进展［J］．空气动力学学报，2014，32(5): 559-574．

10．7638/kqdlxxb-2014．0103．JIANG Y，DENG X G，MAOM L，etal．Progress of aeroacoustic simulationmethod based on HDCS-E8T7 scheme［J］．ACTA Aerodynamica Sinica，2014，32(5):559-574．