一种适合于旋翼前飞非定常流场计算的新型运动嵌套网格方法

王 博，招启军，徐 广，徐国华

（南京航空航天大学 直升机旋翼动力学重点实验室，江苏 南京210016）

0 引 言

旋翼是直升机最主要的气动部件，其气动特性不仅对直升机的飞行性能起着举足轻重的作用，还对振动、噪声等特性产生重要影响。因此，旋翼的空气动力学特性是直升机研究领域的重要研究内容之一，是研发具有我国自主知识产权直升机旋翼设计的前提和基础。目前，关于旋翼气动特性的研究方法主要分为两类：试验方法和数值模拟方法。采用试验方法可以较准确地获得旋翼气动特性，但是存在成本高、周期长等缺点。而数值模拟方法可以有效避免这些缺点，因此，开展旋翼气动特性数值模拟研究具有重要现实意义。

目前常见的旋翼气动特性计算方法主要分为涡方法［1－2］和计算流体力学（CFD）方法［3－6］。涡方法是一种相对成熟的方法，在传统桨叶气动特性的模拟、预测方面具有良好的效果，且具有很高的计算效率。然而，伴随着直升机技术的发展，旋翼的气动外形设计要求不断提高，涡方法由于其基于势流假设，很难精确模拟现代直升机旋翼的气动特性和流场，尤其是桨尖附近的细节流场特征。而基于Euler／N－S方程的CFD方法可以满足新型桨叶流场模拟的精度需要，因而成为了新的研究热点。然而，相对于固定翼飞机的CFD模拟，旋翼CFD方法含有很多技术难点，其中关键的一点是网格方案。旋翼除旋转运动外，还同时存在独特的挥舞、变距等运动，这些运动将对其自身的空气动力学和动力学特性产生重要影响，并直接导致旋翼贴体网格生成的困难，因此目前通常采用嵌套网格方法。但由于随着方位角的不断变化，桨叶与背景网格之间的相对位置也在改变，所以在采用嵌套网格方法计算时必须不断地刷新网格间的嵌套关系，这对嵌套网格方法的效率和鲁棒性提出了严峻的挑战。

使用嵌套网格时需要在网格重叠部分进行插值，因此必须判断网格间的单元对应关系，而判断方法的效率和鲁棒性将直接影响嵌套网格技术的实用性。为解决这一关键技术，国内外学者提出了很多搜索判断方法［7－12］。传统的表面向量法［7］需要预设洞边界，并且在判断每个给定点时都需要寻找与该点最近的交接面的面单元，因而效率不高。射线法［8］对给定点做射线，根据射线与边界交点的数量判断是否在洞内，同样由于需要在判断每一个点时循环计算交接面上的面单元导致计算效率很难提高。洞映射方法［9］建立与原网格对应的辅助均匀笛卡尔网格，用辅助网格近似原网格的挖洞曲面，该方法效率和自动化程度都很高，但是当参考网格较密时会消耗较多的内存。而基于洞映射方法和射线法的目标射线法通过建立垂直于某一方向的辅助网格，判断点在沿该方向的射线与边界交点的数量判断点的属性，同时具有效率高和内存消耗小的优点，但由于混合了洞映射方法和目标射线法算法复杂度有所增加。国内的学者也有提出多种洞单元识别方法［10－12］，获得了一定的效果。

通过吸收洞映射法思想，结合旋翼流场计算网格特点，本文提出了一种新的“透视图”挖洞方法用以解决前飞旋翼流场计算中笛卡尔背景网格和结构网格之间的嵌套关系，与前面提到的多种方法相比较，该方法同时具备节省内存、算法简单、高效、可靠等特点，适合于旋翼悬停／前飞状态的计算。Inverse Map［13］是嵌套网格技术中又一重要技术方法，用于寻找覆盖于给定点在曲线网格上贡献单元，本质也是判断网格间的单元对应关系，因此将此方法推广到Inverse Map的建立过程中，相对于传统建立方法有效地提升了效率。最后，本文采用提出的嵌套网格方法对旋翼悬停／前飞流场进行了计算，验证了该方法的有效性。

1 嵌套网格方法

1.1 “透视图”挖洞方法

该方法的基本思路是：遍历桨叶网格（细网格）中物面点并计算其在粗网格（背景网格）上的对应单元序号，将这些序号保存在相应的数组里，最后由这一数组即可在粗网格上重现对应的桨叶形状，所得到的形状即为洞的最小边界，并可用来标记洞单元。用同样的方法可以通过遍历桨叶网格外围的点获得洞的最大边界。

“透视图”方法主要分为以下四步骤：

第一步，建立桨叶网格以背景网格尺寸为分辨率的数组，并获取相关的网格参数；

第二步，循环网格中壁面点，根据点在背景网格的单元的序号更新已经建立的数组，最后获得物面的“透视图”，该图标明了挖洞的最小范围；

第三步，循环网格外表面点，根据点在背景网格的单元的序号更新已经建立的数组，最后获得物面的“透视图”，该图标明了挖洞的最大范围；

第四步，根据所建立的“透视图”进行洞点的判断。

以二维网格为例，如图1所示，在尺寸为N×M的笛卡尔背景网格上对一翼型网格挖洞，“透视图”方法的基本算法如下：

（1）建立尺寸为N×4的整数数组T，其中N为背景网格的I方向单元数，数组初始化为0。

（2）循环网格的物面点，根据P点坐标计算其所在背景单元的序号i，j。依次作如下运算：

（3）循环网格的外表面点，根据P点坐标计算其所在背景单元的序号i，j（因与物面点循环过程类似，图1并未给出）。依次作如下运算：

图1 二维翼型“透视图”法挖洞示意图Fig.1 Schematic of the top view map method for hole－cutting around 2Dairfoil

图2给出了循环物面点时，数组T中某一列的变化过程。当最终循环完所有物面点，数组T如图3所示。其中数字与图1中阴影部分相吻合。

图2 循环计算 A、B、C点时T（12，：）变化过程Fig.2 Change process of T（12，：）when cycling A，B，and C points

（4）标记洞单元，如图4所示，根据需要可改变算法来选择洞的大小：

a）根据T（i，1）、T（i，4）的值对第i列背景网格单元V（i，j）做标记，可得到最大范围的洞。判断条件如下：

b）根据T（i，2）、T（i，3）的值对第i列背景网格单元C（i，j）做标记，可得到最小范围的洞。判断条件如下：

图3 完成物面点循环后的数组TFig.3 Array Tafter completing cycling points on the wall surface

如图4（a）所示，按照最大范围策略对网格标记洞单元得到的最大范围的洞与网格外边界紧密贴合，这种挖洞策略适用于流场非线性程度较高的场合，例如前飞时桨叶尖部出现激波的位置。图4（b）中显示了采用最小范围内挖洞策略，洞边界恰好严密地包围了桨叶表面网格，这种方法适合于对洞大小限制较强的情况，例如桨叶片数较多时，桨叶根部间距较小的状态。

图4 不同挖洞策略Fig.4 Different strategies of hole－cutting

值得注意的是在生成“透视图”数组时可以选择一次性生成包含所有桨叶的数组，也可以分别生成每片桨叶的数组。采用后一种方法生成的不同“透视图”之间可以发生重叠，因此该方法可拓展用于带机身的旋翼／机身干扰流场计算中。

1.2 挖洞测试

为验证“透视图”法的执行效率和鲁棒性，设计了一组测试算例。

（1）在频率为2.6GHz的 Dual－Core PC上，对两片桨叶进行动态挖洞，背景网格尺寸为120×62×120均匀笛卡尔网格，桨叶网格尺寸为159×13×83，桨叶表面点共有10587个点，方位角由0°到90°，每3°挖洞一次。消耗时间如图5所示，从图中可以看出，当桨叶在不同方位角时挖洞消耗的时间基本相同。

图5 不同方位角挖洞所消耗时间Fig.5 Consuming time of hole－cutting at different azimuthal angles

（2）由于采用N－S方程进行计算时，需要物体表面网格尺寸满足一定的尺寸需求，这会导致网格数量大幅增加。为测试该挖洞方法的速度随桨叶网格尺寸增加的变化，分别生成了桨叶表面点数为1×10587～8×10587的桨叶网格。不同桨叶网格挖洞消耗时间如图6所示。可以看出，消耗时间随着桨叶表面网格点的变化基本上为线性增长。这是因为在建立“透视图”数组时所循环的点只针对桨叶网格的物面点和表面点。如对I×J的二维网格，需要计算的点的数量小于I×2，三维情况下小于I×J×2。

图6 挖洞时间与桨叶网格尺寸关系Fig.6 Relationship between hole－cutting time and blade mesh size

（3）为测试该挖洞方法的速度随背景网格尺寸增加的变化，分别生成了尺寸为120×124×120（2倍）、240×62×240（4倍）和240×124×240（8倍）的均匀背景网格。采用与前面相同的测试方法，在不同背景网格挖洞消耗时间如图6所示。从图中可以看出，背景网格尺寸增加对计算时间的影响明显弱于桨叶网格尺寸。这主要是因为对于均匀背景网格，查找背景网格单元序号的计算量并没有太大变化，而标记洞单元消耗时间占总消耗时间比重较少。对于非均匀直角网格，由于采用了二分法查找单元序号，计算量增加很小。

（4）在前飞状态下，由于桨叶同时存在挥舞运动、变距运动和摆振运动，运动规律较为复杂，这就需要网格方法具有较强的鲁棒性。作为方法验证，选取β＝×cos（ψ）为挥舞运动规律，θ＝×cos（ψ＋）为周期变距，ψ为方位角，应用所建立的方法进行挖洞，并对洞边界进行检测验证所建立的方法的可靠性。测试结果表明在所有给出状态下都可以正确的进行洞点识别。从图7中可以看出不同方位角时挖洞时间基本相同。

图7 计入挥舞及变距运动情况下挖洞消耗时间Fig.7 Consuming time of hole－cutting when considering flapping and pitching motions

1.3 隐含条件的处理及其意义

在使用“透视图”方法挖洞时存在一个隐含的条件，即在需要挖洞的区域内桨叶网格的密度必须大于背景网格。在实际应用中，为在保证计算精度需要的基础上减少计算时间，背景网格的网格点密度往往小于旋翼桨叶的密度。但作为可靠性考虑，对这种情况进行简单的讨论。如图8所示，当不能满足这一条件将有可能导致“遗漏”洞单元的情况有两种。第一种，由于翼型上表面和下表面在第9列上都没有网格点，最终会导致在第9列上遗漏了3个洞单元。而在第7列，由于下表面没有相应的网格点导致少识别2个洞单元。

当出现这种遗漏洞单元的情况，可以采用两种解决办法：

图8 挖洞可能出现的“遗漏”现象Fig.8 The possible missing phenomenon during the process of hole－cutting

（1）可以对物面网格点进行插值细分，保证物面网格点相邻两点间的最大距离小于背景网格上相邻两点间的最小间距即可。

（2）建立网格点密度较低的参考背景网格，在保证隐含的条件的基础上对参考背景网格挖洞，然后按照参考背景网格的洞单元对背景网格单元进行标记。

由2.2节的测试结果可以看出，采用第一种方法会导致计算时间的线性增长，但是在实际编程过程中这一方法实现较为简单，在原始挖洞时间较小的情况下可以采用这一方法。采用第二种方法增加了数据结构的复杂度，但是可以保证消耗时间增长不大，适合于对挖洞时间有较高要求的情况。

从另一方面来看，由于只要满足隐含条件即可保证挖洞的需要，那么就可以通过适当减少物面点的判断来减少计算量，进一步减少挖洞过程消耗的时间。如图1中所示，实际上只需要循环A、C两点即可得到正确的T（12，：）。

1.4 非结构网格上的使用

前述方法只适合于背景网格为结构网格的情况，当背景网格为非结构网格时，由于网格单元间为无序排列，因此不能直接使用该方法进行挖洞。

为解决网格单元无序排列的问题，可以在对应的空间范围内建立结构化的参考网格，在计算前生成背景网格单元与参考网格单元之间的映射关系。通过参考网格映射可以将前述方法应用于非结构背景网格的挖洞过程。

1.5 透视图法在Inverse Map建立中的应用

嵌套网格方法中另一关键技术就是对贡献单元的搜寻，目前的搜寻方法基本都采用Inverse Map（IM）法［13］。在搜索某一点的贡献单元时，采用该方法可以将待搜索的贡献单元数量减少到几个甚至一个。目前在旋翼前飞流场 计算中主要采用刚性桨叶假设，只需要生成一次IM就可以满足使用的需要。但是当采用弹性假设时，桨叶网格会发生变形，必须重新生成。在这种情况下，由于需要多次生成IM，因此IM的生成效率变得非常重要。

将“透视图”法应用于IM的生成需要改变的地方有两点。

（1）数组T尺寸变为NI×MI×4，其中NI、MI为IM的尺寸。

（2）遍历桨叶上所有的网格点。保存在数组 中的参数由背景网格序号改为桨叶网格点在桨叶网格上的序号，所采用判断法则相同。

在建立好的数组中，T（i，j，1）为IM 上（i，j）单元在桨叶网格上对应的I方向序号的下限。T（i，j，2）为I方 向 序 号 的 上 限。T（i，j，3）为J方 向 序 号 的 下 限。T（i，j，4）为J方向序 号 的 上 限。 即T（i，j）对 应 的 贡 献 单元为B（ib，jb），其中T（i，j，1）≤ib≤T（i，j，2）、T（i，j，3）≤jb≤T（i，j，4）。

在计算量方面：首先，由于生成IM时需要循环桨叶网格上所有的单元，因此遍历并计算的点的数量大于挖洞所处理的点。其次，一般情况下IM为均匀网格，因此查找某一点在其上对应的单元序号可直接通过除法计算得到，因此总体计算量无明显变化。

2 流场求解方法

2.1 控制方程

对前飞流场数值模拟，采用以绝对物理量为参数的守恒的积分形式的NS方程为主控方程：

其中：

在以上各式中，q＝（u，v，w）T表示绝对速度，qr表示网格运动速度，E和H分别为总能和总焓，粘性相关项为Txx＝2μux－（2／3）μ▽q，Txy＝μ（uy＋ux），φx＝uTxx＋vTxy＋wTxz＋k∂T／∂x等。其中，V为控制体单元，S为单元面积，n表示单元表面法矢量，t为物理时间，μ、k、T分别为粘性系数、热传导系数和绝对温度。湍流模型采用B－L模型。

2.2 方程离散

本文时间推进采用显式五步Runge－Kutta迭代格式。空间离散采用格心形式的Jameson中心平均有限体积法，相邻单元交接面处的数值通量采用平均计算。另加入由二、四阶混合导数组成的人工粘性项消除中心差分具有奇偶失关联及高频误差的等缺点，避免了由无阻尼引发的非线性（如激波）数值振荡。为加快收敛速度还另外采用了隐式残值光顺等技术。

2.3 计算网格

计算网格采用嵌套网格技术。直接采用上述“透视图”法来确定洞点以及洞边界，并用“透视图”法加速Inverse Map生成。以7A旋翼［14］为例，网格由两部分组成，桨叶网格采用C－O型网格，尺寸为159×30×83，背景网格采用笛卡尔网格，尺寸为166×140×166，背景网格上边界距离桨盘平面为3R，下边界距离桨盘平面为5R，周向边界距离桨毂中心轴为10R。为了比较准确地捕捉桨尖涡，以及减少数值耗散，在背景网格中桨叶运动附近区域进行了加密，该处的网格尺寸为0.1倍弦长。在背景网格上，采用Euler方程求解桨叶的流场和捕捉远场尾迹。

2.4 边界条件

桨叶表面的热力学和动力学边界条件分别取作法向导数为零，即＝0和＝0，在远场取基于Ri－emann不变量的远场边界条件。对于旋翼网格及背景网格之间的流场信息交换由三线性插值来完成。

3 旋翼前飞流场计算及结果分析

3.1 UH－60A旋翼

为验证本方法在新型桨尖上使用的可靠性，选取UH－60A旋翼算例，由于公开发表的文献中给出的试验状态数据不全，因此选择悬停状态试验［15］结果进行对比。该模型旋翼有4片桨叶，展弦比15.3。该旋翼的桨叶外形较7A旋翼更为复杂，桨叶尖部具有20°后掠；桨叶分为三段，由两种不同的翼型组成；负扭转分布为非线性负扭转，最大扭转为13°。计算状态为：Mtip＝0.628，θ0.75＝9°。

图9给出了0.775R和0.92R处的桨叶表面压强系数分布。从图中可以看出，旋翼靠近中间的部分和靠近外侧的计算结果与试验值吻合良好。图10给出了悬停效率的计算值与试验值的比较，可以看到计算结果的趋势和试验值一致。这表明所建立的挖洞方法结合基于N－S方程的求解器，能有效地模拟新型桨尖旋翼流场。

图9 UH－60A旋翼桨叶表面压强系数分布Fig.9 Distributions of pressure coefficient on the blade surface of the UH－60Arotor

图10 UH－60A旋翼悬停效率随拉力系数的变化Fig.10 Predictions of Figure of Merit versus CT／σfor the UH－60Arotor，compared with experimental data

3.2 HELISHAPE 7A旋翼

为了验证本文给出的"顶透视图"挖洞方法的高效和可靠性，计算了前飞状态旋翼的气动特性。选用有试验结果可以对照的HELISHAPE 7A旋翼［14］作为算例。该旋翼有4片桨叶，直径4.2m，展弦比15，实度σ≈0.0849，桨叶平面形状为矩形。HELISHAPE 7A旋翼的桨叶气动外形较复杂，桨叶分为三段，每段由不同的翼型组成，且有非常规的几何负扭转（－3.95°）。计算状态为：Mtip＝0.616，μ＝0.167。

为验证所建立方法的的有效性，分别对比了“透视图”法和从全局遍历方法查找对应点方法所消耗的时间。从表1中可以看出，采用“透视图”法后，无论是洞点单元标记还是Inverse Map生成速度都有了很大的提高。速度的提高主要因为计算时只需要循环计算桨叶网格点边界点在背景网格上的坐标一遍即可，避免了全局查找方法中的大量查找过程。其次，计算量中主要部分在于查找点P在的背景网格单元序号，这在本质上来说是一个有序数组的查找问题，可用成熟的查找算法加速。

表1 使用与不使用“透视图”方法所用时间的比较Table 1 Comparisons of consuming time between using and not using top－view map method

图11给出了三个不同方位角上桨叶桨叶表面压强系数与试验值结果［14］的对比。如图所示ψ＝90°时在桨叶尖部，结合该挖洞方法的N－S方程计算得到的桨叶上下表面压强系数与试验值较为接近；ψ＝180°时桨叶表面压强系数的分布，结果优于90°时计算结果。图11（c）给出方位角为270°时桨叶表面压强系数的分布。计算结果表明采用“透视图”法的流场求解器可以在保证计算结果正确性的前提下有效地提高旋翼前飞非定常流场的计算效率。

图11 7A旋翼前飞状态桨叶表面压强系数分布Fig.11 Surface pressure coefficient of the 7Arotor in forward flight

4 结 论

结合旋翼流场实际计算情况，提出了一种快速确定嵌套网格中不同网格单元之间对应关系的方法，通过多种算例验证，得到如下结论：

（1）“透视图”挖洞方法采用简单的计算方法进行网格关系判断，可快速完成网格间嵌套关系的建立，满足旋翼悬停／前飞流场计算时所需要的动态挖洞要求。随着网格尺寸的增加，计算量增加较小，同时该算法具有节省内存、鲁棒性强的特点，具有重要工程实用价值。

（2）基于新的挖洞技术，建立了一个嵌套网格CFD方法，较为成功地应用于直升机悬停／前飞状态旋翼流场的计算，为进一步开展先进气动外形旋翼数值模拟研究打下了良好的基础。

（3）将“透视图”方法应用于Inverse Map的生成可以大幅度减少建立Inverse Map所需的时间。

［1］BAGAI A，LEISHMAN J G.Rotor free－wake modeling using apseudo implicit relaxation algorithm［J］.Journal ofAircraft，1995，32（6）：1276－1285.

［2］李春华，徐国华.悬停和前飞状态倾转旋翼机的旋翼自由尾迹计算方法［J］.空气动力学学报，2005，23（2）：152－156.

［3］KANG H J，KWON O J.Unstructured mesh Navier－Stokes calculations of the flow field of a helicopter rotor in hover［J］.JournaloftheAmericanHelicopterSociety，2002，47：90－99.

［4］解福田，宋文萍，韩忠华.低耗散格式在旋翼前飞气动噪声预测中的应用研究［J］.西北工业大学学报，2009（2）：151－156.

［5］牟斌，肖中云，周铸，等.直升机旋翼悬停流场的粘性数值模拟［J］.空气动力学学报，2009，27（5）：582－585.

［6］叶靓，招启军，徐国华.基于非结构嵌套网格和逆风格式的旋翼悬停流场数值模拟［J］.空气动力学学报，2009，27（1）：62－66.

［7］SLOTNICK J P，KANDULA M，BUNING P G.Navier－Stokes simulation of the Space Shuttle launch vehicle flight transonic flowfield using a large scale chimera grid system［A］.12th AIAA Applied Aerodynamics Conference［C］.Colorado Springs，Colorado，USA，1994.

［8］MEAKIN R L.Object X－rays for cutting holes in composite overset structured grids［A］.15th AIAA Computational Fluid Dynamics Conference［C］.Anaheim，California，USA，2001.

［9］ROGERS S E，DIETZ W E，SUHS N E.PEGASUS 5：an automated preprocessor for overset－grid computational fluid dynamics［J］.AIAAJournal，2003，41（6）：1037－1045.

［10］李亭鹤，阎超.一种新的分区重叠洞点搜索方法－感染免疫法［J］.空气动力学学报，2001，19（2）：156－160.

［11］李亭鹤，阎超，李跃军.重叠网格技术中割补法的研究与改进［J］.北京航空航天大学学报，2005，31（4）：402－406.

［12］杨文青，宋笔锋，宋文萍.高效确定重叠网格对应关系的距离减缩法及其应用［J］.航空学报，2009，30（2）：205－212.

［13］MEAKIN R.A new method for establishing intergrid communication among systems of overset grids［A］.10th AIAA Computational Fluid Dynamics Conference［C］.Honolulu，Hawaii，USA，1991.

［14］POMIN H，WAGNER S.Navier－Stokes analysis of helicopter rotor aerodynamics in hover and forward flight：Rotor wakes［J］.JournalofAircraft，2002，39（5）：813－821.

［15］LORBER P，STAUTER R，LANDGREBE A.A comprehensive hover test of the airloads and airflow of an extensively instrumented model helicopter rotor［A］.45th AHS Annual Forum［C］.Boston，Massachusetts，USA，1989.