结合并行策略与距离减缩法的搜索算法及其在CAA中的应用

， ， ，

(1.中国飞机强度研究所 航空声学与动强度航空科技重点实验室，西安 710065；2.西北工业大学 航空学院，西安 710072)

1 引 言

对于计算气动声学CAA(Computational Aero-Acoustics)域和计算流体力学CFD(Computa-tional Fluid Dynamics)域之间的耦合问题，可将CFD计算结果作为CAA模拟的一个输入边界条件进行处理。由于两者离散格式分辨率存在区别，正确的边界条件对于数值模拟结果至关重要。采用混合CAA方法来预测气动噪声，在CAA网格上获得湍流速度的过程，涉及到CFD网格与CAA网格的时均和脉动等流场信息的传递。而流场信息从流场网格传递到声场网格的过程，需要经过网格点的搜索和数据的插值两个步骤。如何快速有效地实现数据从CFD网格到CAA网格的传递，是混合CAA方法的一个重要环节。

针对不同计算域间的信息传递问题中插值网格点的识别和插值单元寻址的关键技术，国外对此做了大量的研究工作，主要方法有表面法向量方法[1]、洞映射法[2，3]和目标X射线法[4]等，而这些方法主要应用在重叠网格中。我国学者在这些算法的基础上作了改进，提出了若干高效的搜索算法，其中，李亭鹤等[5]提出了一种新的分区重叠洞点搜索方法，即感染免疫法，实例验证表明了该方法具备有效的洞点搜索能力；田书玲等[6]发展了动态非结构重叠网格技术，提出了基于阵面推进的相邻单元搜索算法和适用于高阶空间离散格式的网格间插值关系定义方法，分别对定常M6机翼扰流和机翼下外挂物投放过程的非定常扰流进行了模拟，验证了方法的有效性；杨文青等[8]提出了距离减缩法，将三维搜索化为准一维搜索，大幅提高了搜索效率，算法通用性强，易于实现，适用于计算流体力学三维非定常计算中有相对运动的重叠网格系统。

通过对CFD网格与CAA网格特点的研究，以及对上述方法优缺点的分析，在文献[8]距离减缩算法的基础上，本文结合并行策略对该算法作了改进，实现了从流场网格到声场网格之间的信息快速传递，通过二维标模验证了方法的计算效率。在此基础上，搭建了气动噪声混合预测方法，对NACA0012翼型的后缘气动噪声进行了预测分析。

2 结合并行策略与改进的距离减缩法

由于距离减缩算法计算效率高、适用范围广(结构与非结构网格)及算法复杂度低等特点，在嵌套网格洞点识别中得到广泛应用。

以嵌套网格(图1)曲线网格点O在笛卡尔网格的点搜索为例，对距离减缩方法进行说明。从笛卡尔网格的任意一点开始，设为点A。首先，计算并比较A点周围的点与点O的距离，可得到距离最近的点B；接下来从点B开始，继续计算点B周围点到点O的距离，得到距离最近的点C；以此类推，可以由路线A-B-C-D-E-F-G-H-I寻找到目标点O。这种方法类似于一种梯度下山法，该算法在三维问题上的搜索效率更加明显。

对于计算气动声学的数值模拟问题，由于CAA网格规模大，采用上述距离减缩算法效率较低，不能有效满足实际的计算需求。为了改善这一问题，本文建立了耦合并行策略与改进的距离减缩方法结合的点搜索算法。

采用ICEM-CFD软件生成CFD与CAA网格，导出的数据格式为Plot3d格式。为了节省CAA网格点的搜索时间，在保证CFD网格与CAA网格的拓扑一致前提下(即网格块顶点位置一致，如图2所示)，采用并行策略将网格块编号按照网格点数分配到各个子进程中。

在实现这一策略过程中，当网格曲率变化较大时，由于CFD网格分布与CAA网格分布不一致，存在于CAA 网格块中的网格点有部分没有完全位于CFD 网格对应的网格块中。为了应对这种问题，在点搜索算法中需对这该情况进行特殊处理。

图1 距离减缩算法在嵌套网格中的应用

Fig.1 Application of distance reduction algorithm

图2 CFD网格与CAA网格的拓扑分布

Fig.2 Topology of CFD mesh and CAA mesh

以二维圆柱为例，如图4所示，将网格划分为33个网格块并依次编号，当计算过程中调用进程数为16时(主进程编号为0，子进程编号分别为1，…，15)，各子进程中分配到的网格块编号及对应的网格点总数列入表1。

表1 各子进程分配情况说明

Tab.1 Grid block allocation of child process

子进程数Block编号在子进程中的存储各子进程中总网格点数174181293737312282542517175152525627252572,18,242567823,30,31256793,4,52431101,17,3224311119,20,2624191210,11,2223671313,162450148,14,292439156,21,28,332568

图3 CAA块网格点落入位置出错示意图

Fig.3 Grid location error of CAA block

图4 二维圆柱Block编号分布示意图

Fig.4 Block number distribution of 2D cylinder

针对各子进程中网格块的处理，采用一种改进初始值选取方式的距离减缩搜索技术。结合具体情况，对于图5的CAA网格点A，寻找其在CFD网格中的所在网格单元(称为宿主单元)，具体实现步骤如下。

(1) 对所有CAA网格点采用flag=1进行标记。

(2) 对于flag=1的第一个CAA网格点，记为点A，从CFD网格区域中任取一点，记为B，计算点A与点B的距离lA B。

(3) 计算点B周围的4个点(对三维而言，点B周围为6个网格点)到点A的距离，选距离点A最近的点，记为C。点C到点A距离记为lA C。

(4) 比较距离lA B与lA C。如果lA B>lA C,则将点C标记为B,转向步骤(3)继续进行搜索；如果lA B≤lA C,则点B为距离点A最近的点，完成第一个点的搜索，转向步骤(5)对网格块中下一个网格点进行处理。

(5) 对于下一个CAA网格点A′，如果A′点与A点是相邻的，那么A′在CFD网格中宿主单元和A点的宿主单元相距很近。因此对于接下来的CAA网格点选取前一步的宿主单元作为初始网格点，循环步骤(1～3)，这样只需要寻找几步甚至是一步就可以找到A′点的宿主单元。

图5 CAA网格点A

Fig.5 Point A of CAA mesh

3 基于形函数的插值算法

在完成CAA网格点的搜索后，本文采用基于形函数的拉格朗日插值方法进行CFD信息的插值。由于计算网格采用的是结构网格，在曲线坐标系下网格会出现扭曲和变形等现象，不宜采用常规的双线性等插值手段。而基于形函数的插值方法，则是采用形函数将形状大小不同的单元通过标准化变换成统一的标准化单元，如图7所示。将网格单元进行标准化后，单元内部点的流场物理量可以通过网格单元顶点的物理量和单元顶点的插值形函数进行加权求和得到

(1)

图6 CAA网格点A寻点示意图

Fig.6 Diagram of pseudo mesh

图7 网格单元的标准化示意图

Fig.7 Standardization of grid cells

(2)

(3)

(4)

式中 (xC A A,yC A A,zC A A)为CAA网格的插值点坐标，(xi,yi,zi)是宿主单元顶点i的坐标。利用牛顿迭代法对其进行求解，获得局部坐标为

Vn+1 =Vn-F(Vn)/(∂F/∂V)n

(5)

4 搜索算法的效率验证

为了验证本文算法的网格点搜索效率，利用麦道公司30P30N二维三段翼标模进行验证分析。模型的CFD网格和CAA网格如图8所示，CFD网格量为1064356，附面层第一层网格长度为 1e -6，而CAA网格量为2425924，最小网格尺度为 0.001。

在CPU为Intel Xeon CPU E5-2690，主频2.6 G的计算机上，对比了传统的距离减缩算法和基于MPI并行策略的改进距离减缩算法的插值效率，具体结果列入表2。由表2可知，随着子进程数的提高，本文提出的算法比传统距离减缩算法，效率更高，所需时间更短。该方法为CFD网格与CAA网格间信息的快速传递奠定了良好的基础。

图9给出了缝翼处CFD网格及CAA网格的湍动能信息对比结果，两者在分布位置和量级上保持高度一致。从图9的插值结果可以说明，本文所建立的流场与声场信息传递方法具备处理复杂网格的能力。

5 CAA气动噪声计算流程

CAA方法的流程如图10所示，分为以下几个步骤。第一步，在进行气动噪声计算之前，针对CFD计算和CAA计算的特点分别准备两套网格，命名为CFD网格和CAA网格；第二步，基于CFD网格，通过求解可压缩的N-S方程获得流场信息；第三步,按1，2节中叙述的方法将CFD网格上的流场信息传递到CAA网格上；第四步,根据流场信息，在CAA网格上采用湍流速度生成方法[9]构造湍流速度场；第五步,根据湍流速度信息，数值求解带源项的声传播方程获得声场。

表2 不同算法的计算效率对比

Tab.2 Comparison of computational efficiency

子进程数本文算法/s传统距离减缩算法/s二分法/s04.025.12360.2412.08--31.17--70.73--150.52--

图8 计算网格示意图

Fig.8 Schematic diagram of computational grid

图9 前缘缝翼处的湍动能分布

Fig.9 Distribution of turbulent kinetic energy of slat

6 应用算例

应用本文方法计算NACA0012翼型后缘宽频气动噪声。NACA0012翼型参考长度为l=0.3048 m，马赫数M=0.1156,雷诺数Re=8.3e5。流场网格如图11(a)所示，物面法向网格最小距离为 3e -6，网格单元总数为124094。声场网格如 图11(b)所示，其中物面点有430个，第一层网格距离为0.0015，网格总数为182160。

由于流场计算与声场计算对计算方法的要求不同，因此导致两者的网格分布不同。为了能高效可靠地将流场网格上的信息赋值给声场网格，采用本文1，2节所述方法进行处理，具体效果如图12所示。图12分别为CFD网格和CAA网格流向平均速度传递后的结果对比。可以看出,CFD网格上的流场信息能真实有效地传递到CAA网格上。

图13为采用第4节构建的CAA流程计算的翼型上半部分湍流速度分布。可以看出,两方向湍流速度均呈现正负交替的分布。采用带源项的声传播方程进行声源的传播过程求解，声场计算全场采用统一的无量纲时间步长Δt=0.0002,求解步数为50万步。图14为某一时刻的声压分布云图，图(a)为本文方法计算结果，图(b)为文献[9]计算结果。可以看出,两者计算的声压以翼型中心线呈反对称分布，声波以正负交递的形式向远场进行传播。图15为位于翼型后缘上方1.5倍翼型参考长度处的频谱曲线。其中，图15中Cale为本文计算结果，Ref为参考文献[9]计算结果。通过对比可以看出,本文计算结果除高频段略大于文献[9]，其趋势和声压级量级基本相当。

图10 CAA方法框架

Fig.10 Flowchart of CAA method

图11 NACA0012翼型计算网格图

Fig.11 Computational grid of NACA0012 airfoil

图12 流向速度分布

Fig.12 Distribution of flow velocity

图13 湍流速度分布

Fig.13 Distribution of turbulent velocity

图14 无量纲化的声压分布

Fig.14 Distribution of sound pressure

图15 (0.3048，0.4572)点处1/3倍频带的频谱曲线

Fig.15 1/3 octave frequency spectrum curve

7 结 论

通过上述数值计算结果，可得出以下结论。

(1) MPI并行策略与改进的距离减缩法结合的点搜索算法，比单纯的距离减缩法效率更高。

(2) 基于本文传递算法的CAA流程能有效模拟NACA0012后缘宽频气动噪声。

(3) NACA0012的后缘宽频噪声呈现相互交替现象向空间进行传播。

参考文献(References):

[1] Stotnick J P,kandula M.Buning P G.Navier-Stokes simulation of the space shuttle launch vehicle flight transonic flow field using a large scale chimera grid system[A].12t hApplied Aerodynamics Conference[C].1994.

[2] Steger J.Notes on Composite Overset Grid Schemes Chimera[D].University of California,1992.

[3] Suhs N E,Rogers S E,Dietz W E.PEGASUS 5:an automate pre -processor for overset-grid CFD[J].AIAAJournal,2013,41(6):1037-1045.

[4] Meakin R L.Object X-rays for cutting holes in composite overset structured grids[A].15t hAIAA Conputational Fluid Dynamics Conference[C].2001.

[5] 李亭鹤,阎 超.一种新的分区重叠洞点搜索方法—感染免疫法[J].空气动力学学报,2001,19(2):156-160.(LI Ting-he,YAN Chao.A new method of hole -point search in grid embedding technique[J].ActaAerodynamicaSinica,2001,19(2):156-160.(in Chinese))

[6] 田书玲,伍贻兆，夏 健.用动态非结构重叠网格法模拟三维多体相对运动绕流[J].航空学报,2007,28(1):46-51.(TIAN Shu-ling,WU Yi-zhao，XIA Jian.Simulation of flows past multi-body in relative motion with dynamic unstructured overset grid method [J].ActaAeronauticaETAstronauticaSinica,2007,28(1):46-51.(in Chinese))

[7] Park M A,Aftosmis M J,Compbell R L,et al.Summary of the 2008 NASA Fundamental Aeronautics Program Sonic Boom Prediction Workshop [J].JournalofAircraft,2013,51(3):234-245.

[8] 杨文青,宋笔锋,宋文萍.高效确定重叠网格对应关系的距离减缩法及其应用[J].航空学报,2009,30(2):205-212.(YANG Wen-qing,SONG Bi-feng，SONG Wen-ping.Distance decreasing method for confirming corresponding cells of overset grids and its application[J].ActaAeronauticaetAstronauticaSinica,2009,30(2):205-212.(in Chinese))

[9] Ewert R,Appel C,Dierke J,et al.RANS/CAA based prediction of NACA 0012 broadband trailing edge noise and experimental validation [A].15t hAIAA/CEAS Aeroacoustics Conference [C].2009.

[10] 余培汛，白俊强,郭博智，等.一种可用于机翼优化设计中的气动噪声预测模型[J].计算力学学报,2014,31(4):537-544.(YU Pei-xun,BAI Jun-qiang,GUO Bo -zhi,et al.A new aerodynamic noise model for aerodynamic optimization design [J].ChineseJournalofComputationalMechanics,2014,31(4):537-544.(in Chinese))

[11] Menter F R,Egorov Y.The scale -adaptive simulation method for unsteady turbulent flow predictions.part I：theory and model description [J].Flow,TurbulenceandCombust,2010,85(1):113-138.