旋翼流动的块结构化网格自适应方法

肖中云，郭永恒，张 露，崔兴达

(中国空气动力研究与发展中心，绵阳 621000)

0 引 言

旋翼尾涡是流体在旋转桨叶作用下产生的一种特殊的旋涡流动现象，是由桨叶梢部脱出的集中涡和后缘脱出的尾迹面组成的流动演化系统，其中桨尖涡在尾迹流动中占据主导作用，是旋涡尾迹流场的骨架[1]。由于桨叶旋转、挥舞等非定常运动特点，桨尖涡的生成及对后续桨叶的干扰产生严重的非定常气动载荷，并诱发桨涡干扰（blade-vortex interaction,BVI）噪声。当前，用CFD 方法模拟旋翼尾涡仍是一大难点[2-4]，通常旋翼桨尖涡直径约为桨尖弦长的10%，而模拟桨尖涡的网格尺度则可能小到桨尖弦长的1%。由于空间网格尺度不足以匹配旋涡尺寸导致旋涡的数值耗散过大被抹平，影响到桨涡干扰相关的气动载荷、噪声预测的准确性。

自适应网格加密（adaptive mesh refinement, AMR）是解析局部流动特征的最有效方法之一，理论上可以根据流场特征对网格分辨率进行动态调整，是大范围内模拟旋涡空间演化的理想方法。在传统计算方法里面，多块对接结构网格由于有严格的JKL 指标关系，很难实现网格的局部加密。非结构网格在数据结构上比较灵活，但网格加密后单元质量下降及网格信息存储量大是其面临的难题。笛卡尔网格由于具有空间网格正交和易于进行网格自适应加密等特点得到了较好发展[5-6]，被誉为是CFD 技术的未来发展方向之一。笛卡尔网格自适应目前有单元自适应和块自适应（block structured AMR, SAMR）两种方法。单元自适应根据流场计算需要对网格单元进行加密；块自适应[7]指网格加密前后都具有块结构化的特点，待加密区域不是单个的网格单元，而是规则的立方体区域。目前发展的自适应加密第三方库中，采用单元加密方法的如libMesh[8]、P4est[9]等，采用块自适应网格的如CHOMBO[10]、PARAMESH[11]和SAMRAI[12]等。这些三方库将自适应网格加密和并行负载平衡等函数封装为独立模块，将其与所解决的具体物理问题和算法隔离开来，在流体力学、天文学、宇宙学、地质学等多个领域得到了广泛应用。

块结构笛卡尔网格结合了笛卡尔网格与结构网格的优点。首先利用笛卡尔网格自适应能力强的特点，实现空间任意区域的网格自适应加密；其次网格具有结构化特征，就使得网格单元的相邻关系尤其简单，并且结构网格的一些高阶重构计算方法能够得到继续使用，从而获得较高的流场计算效率。本文从上述特点出发，探索了块结构笛卡尔网格在旋翼计算中的生成方法，通过结构贴体网格加背景笛卡尔网格的双网格技术构建计算域，发展了背景网格的几何自适应技术和基于旋翼尾迹的流场自适应技术，采用该方法在UH-60A 旋翼上构建了悬停和前飞流场的自适应网格，表明当前方法的可行性。

1 块结构笛卡尔网格生成

1.1 笛卡尔网格划分方法

块结构笛卡尔网格划分方法如图1 所示。首先定义包围整个计算域的初始长方体网格块，该网格块内网格为三个方向均匀正交排列的笛卡尔网格，网格间距分别为（ Δx,Δy,Δz），并且约定每个方向的网格单元个数为偶数。网格生成是网格块不断细分的递归过程，当判定一个网格块需要加密时，该网格块在每个坐标方向上一分为二。这样经过一次网格细分，二维情况下一个网格块被分裂为四个网格块，三维情况下一个网格块被分裂为八个网格块。新产生的子网格块在三个方向上的网格维数同父网格保持不变，网格间距缩小为原来的1/2，当父网格块在某个方向的网格单元数为偶数的情况下，子网格块在这个方向上正好两个单元对应父网格块的一个单元。图1 显示了笛卡尔网格划分的前三层网格，其中第三层网格只对第二层的一个网格块进行了加密，这样第一、二、三层网格的网格块数分别为1、8、8 块，总的网格块数为17 块。

图1 块结构网格划分方法示意图Fig. 1 Schematic of block structured grid splitting

根据上述网格划分特点，三维笛卡尔网格块采用八叉树数据结构进行管理。如图2 所示，八叉树的每个节点表示一个正方体的体积元素，顶部节点称为根节点，也是体积最大的节点。每个节点有八个子节点，这八个子节点所表示的体积元素加在一起就等于父节点的体积。没有后代的节点称为叶子节点，八叉树叶子节点代表了分辨率最高的情况。例如将空间某个区域的分辨率设成（ Δxt,Δyt,Δzt），那么覆盖该区域的网格块将不断细分，直到每个叶子网格间距将小于目标值为止。八叉树的叶子节点覆盖了整个计算域，图2 显示的叶子节点个数为15。一般情况下流场计算只需要在叶子节点上进行，当采用多重网格等加速收敛计算方法时，中间节点及根节点可以参与到计算中。

图2 八叉树数据结构示意图Fig. 2 Schematic of octree data structure

采用上述方式生成的块结构化笛卡尔网格有很多优点。首先是网格描述非常简洁、存储量极小，初始网格定义了网格块中心点坐标和三个方向的尺度及网格间距，新生成网格块层级和中心点坐标由上一级网格块得到，其他信息如网格点坐标、单元体积、单元面积等都可以直接用解析式给出。这样不需要存储每个网格点的信息，从而大大减少对内存的需求。其次是和结构化网格一样，块结构化的笛卡尔网格流场变量按指标顺序连续存放，有利于充分利用高速缓存提高存取效率，并且非常容易在三个方向上构造高阶插值模板单元，以很小的成本代价实现高阶格式。再次，块结构笛卡尔网格由粗网格不断细分得到，子网格和父网格之间构成多重网格的嵌套关系，可以根据此关系设计多重网格方法，提高流场计算的收敛速度。

1.2 多层空间填充Z 曲线

当前自适应网格方法可以简单地类比于给网格增加或减少补丁，加补丁就是给网格加密，揭补丁就是使网格稀疏。空间填充曲线技术是为了高效存储和管理空间网格的一种方法，它的目的是将多维空间数据转换到一维空间上，并通过转换后的一维空间索引值存储和查询多维数据。除了有降低维度的特性外，空间填充曲线还具有数据聚类特性，其特点是将空间上邻近的网格单元映射为一维曲线上尽可能接近的点，因此只需要访问查询点的邻近点，就能够获得网格单元的最近邻单元。常用的空间填充曲线有Z 曲线[13]、Hilbert 曲线和Gray 曲线等，其中Hilbert 曲线和Gray 曲线涉及到方向旋转，映射过程比较复杂，聚类特性更优。相比之下Z 填充曲线的映射过程比较简单，几何空间的网格相邻关系更容易确定而被网格方法所广泛采用。本文针对自适应网格特点将空间填充Z 曲线进行了多层表示，即填充曲线包含了叶子节点、它们的父节点及根节点等全部网格，将所有网格层展开后得到网格关系如图3 所示。保留非叶子节点网格的目的是能够自由地对网格进行加密或者稀疏变换，从图3 可以看到，通过延拓或者插值计算，LEVEL 1 和LEVEL 2 的流场信息交换是非常方便和直接的。图3 同时给出了表示多层网格的空间填充Z 曲线。其编码规则如下，编码首先从根节点开始，对于任何一个节点，首先判断其是否存在子节点，如果有则指向子节点，如果没有则指向同一级的邻居节点，邻居节点之间用Z 曲线连接，当子节点为Z 曲线最后一个节点时回到父节点。注意这里的节点不重复计数，如图3 中当出现节点5 指向节点1 时，由于节点1 已经存在于序列中而略过，因此5 的下一个节点为6。当采用上述方式索引以后，可以看到，任一节点的邻居节点可能是兄弟（相同父节点），也可能是堂兄弟（父节点兄弟的孩子），但总的来说，仍具有邻居节点位于空间填充曲线附近的聚类特性。

图3 多层网格的Z 填充曲线Fig. 3 Z filling curves of multi-layer grids

在用空间填充曲线进行排序以后，对网格进行并行分区、满足负载平衡就变得非常简单，由于每个网格块的网格单元数相同，因此满足负载平衡的并行分区只需要将一维的空间填充曲线进行等量剖分就可以了。当然也可以考虑到叶子节点和非叶子节点之间计算量的差异，给不同节点赋予不同的权重，采用加权平均的方法进行分区。图4 显示了将21 个网格块分配到5 个进程的分区情况。

图4 多层网格的并行分区Fig. 4 Parallel partition of multi-layer grids

1.3 网格加密判则与加密方法

采用一分为八的方式加密网格将导致网格量呈几何级数增长。为了让网格总量可控，本文采取的策略一是限制网格最大层数，一旦网格层级超过最大值则加密被终止，二是建立合适的加密准则，将网格加密局限在几何复杂、流动变化剧烈的局部区域，并且网格随流动特征变化的自适应加密。

限制网格最大最小层数的作法实际上是限制了网格单元的最大最小尺寸，最大尺寸规定的是远场的网格尺度，与计算域的大小相关；最小尺寸规定的是模拟流动特征的网格最小尺度，比如对于旋翼流动来说，网格加密是为了实现桨尖涡的模拟，因此可以根据桨尖弦长测算出桨尖涡的涡核直径，并由此估算出能够较好捕捉旋涡的网格最小尺度。

网格的局部加密包括了基于几何特征加密和基于流场特征加密两种类型。旋翼模拟一般将笛卡尔网格作为背景网格，笛卡尔网格与近场网格之间构成重叠关系，此时近场外边界网格尺度就作为重叠区笛卡尔网格加密的目标尺度。流场特征加密主要是基于旋翼桨尖涡的识别与局部加密，目前旋涡识别方法发展有很多种[14-15],常见的Q-判据定义方法如下：

其中 Ω表示旋转张量，S表示变形率张量。为了防止网格密度变化过渡剧烈，限定相邻笛卡尔网格块的网格尺度相差最大一倍。即当对某一网格块进行加密时，还需要判断该网格块的相邻块是否存在密度差大于一倍的情况，如果存在，则该相邻块也被列入到待加密列表中，以此类推，直到所有块都满足条件为止。需要注意的是，这里所谓的相邻块包括了面相邻块、棱线相邻块和角点相邻块。对于长方体来说，相邻块的数量最多情况下达26 个。图5 以二维为例对网格加密过程进行了说明，图5（a）显示了待加密网格块（用十字虚线表示），同时识别出周围三个网格块存在密度差大于1 倍的情况（桔色显示）；图5（b）对（a）中桔色显示的三个网格块进行了加密，同时又识别出周围两个网格块存在密度差大于1 倍的情况；图5（c）为最终的网格加密情况，可以看出每组相邻块的网格密度最大相差1 倍；图5（d）为加密后网格的空间Z 曲线填充情况。图5（d）只对叶子节点进行了空间填充的Z 曲线显示，显示多层网格的Z 填充曲线如图6 所示。从图中可以看到，当前网格一共包含了5 层，除了最后一层全部是叶子节点外，其他层凡是有向下箭头线的均为非叶子节点，即表明当前块有更细的网格剖分。图5（d）可以理解为图6 在平面内的投影，并且删除非叶子节点后得到的图形，这些非叶子节点不直接参与流场计算，但在方便网格组织与自适应加密中起到作用。

图5 网格加密过程示意图Fig. 5 Schematic of mesh refinement

图6 局部加密网格的Z 填充曲线Fig. 6 Z filling curves for local mesh refinement

2 算例与分析

2.1 旋翼近场贴体网格

为验证网格模型，本文选择生成UH-60A 旋翼的计算网格，模拟状态包括旋翼的悬停状态和前飞状态。UH-60A 旋翼共包括四片桨叶，旋翼半径R=8.178 m，桨叶展弦比15.5，桨叶具有非线性扭转、桨尖后掠等现代旋翼特征。桨叶贴体网格采用多块对接结构网格（见图7），其中桨叶剖面拓扑结构为O 型，桨叶表面网格为四边形单元，物面边界层区域采用大拉伸比网格进行刻画。

从图7 中可以看到，旋翼网格分为了4 个相对独立的子网格，由于桨叶在一周旋转运动过程中，即使是刚性假设桨叶，也需要进行旋转、挥舞、变距等运动，桨叶与桨叶之间的相对位置会发生变化，将单片桨叶网格独立作为一个子网格，这样就可以让桨叶网格随桨叶一起运动，每个时间步的子网格变化通过多次旋转变换得到。近壁区流场计算由贴体网格负责进行，背景网格深入到壁面的部分经过重叠挖洞以后被当作洞内点不参与计算。

图7 桨叶贴体网格Fig. 7 Body fitted grids of the blade

2.2 背景网格重叠边界的几何自适应

在生成近场贴体网格以后，背景网格采用笛卡尔网格生成方法，即采用双网格技术构建流场计算网格。两组网格的计算域划分如图8 所示，考虑到当前旋翼半径R= 8.178 m，定义背景网格的外廓尺寸为x,y,z∈[-50 m, 50 m]的立方体区间，网格加密层级定义为9 级，这样最小网格单元尺度约为 Δx=0.048 m。图中立方体代表的是初始笛卡尔网格块，块内每个方向的网格单元个数为8 个。

图8 贴体网格与背景网格Fig. 8 Body fitted and background grids

背景笛卡尔网格具有能够自动生成的特点，从图8 所示的初始笛卡尔网格块开始，不断对网格进行细分，直到核心区域网格密度满足模拟要求为止。在这里背景网格与贴体网格之间构成重叠关系，要求重叠区域两组网格的网格尺度相当，因此该过程也被称之为网格几何自适应。由于背景网格要满足与贴体网格重叠插值的条件，这里以贴体网格的外边界作为特征对象，对背景网格块依次进行判断，如果特征对象与网格块相交或者完全落在网格块的内部，则判定该网格块为待加密网格块。对初始网格块来说，定义网格层级LEVEL 0，在经过一次加密后生成8 个子网格块，定义为LEVEL 1，以此类推，直到网格块的网格尺度小于特征对象的网格尺度为止。图9 显示的是经过几何自适应以后得到的笛卡尔背景网格分布，可以看到经过不断细分以后，网格密度大的区域向近场集中；此外网格呈块状分布，块内部网格均匀分布，相邻网格块的网格间距比值最大为2：1。图10 显示了背景网格和贴体网格交界区域的网格分布。可以看到背景网格在桨叶根部和尖部都进行了加密，在子网格外边界区域，加密准则要求背景网格与贴体网格的网格尺度相当，满足重叠插值对网格的要求。

图10 近物面区域的网格分布Fig. 10 Grid distribution in the near wall region

2.3 基于旋翼尾迹的流场自适应

除了前面的几何自适应以外，发展自适应网格的主要目的是根据流场特征对网格进行局部加密。为验证网格对旋翼流场的自适应能力，本文将通过尾迹模型方法得到尾迹分布作为流场特征，对背景网格进行局部加密。

2.3.1 旋翼悬停状态

悬停模拟状态为旋翼转速Ω= 255 r/m，拉力系数CT=0.0069。旋翼尾迹用Landgrebe 尾迹模型表示，具体公式如下：

其中：rtip、ztip分别表示桨尖涡的无量纲径向位置和旋翼轴向位置（用旋翼半径R无量纲化），Nb表示旋翼桨叶片数， ψw为桨尖涡尾龄角， θtw为桨叶线性扭转角，其他系数定义如下：

图11 为Landgrebe 尾迹模型计算得到的旋翼悬停状态下单片桨叶的桨尖和桨根尾迹分布，可以看到桨尖和桨根涡呈螺旋状向旋翼下方发展，桨尖涡的下降速度更快，并且在下降的同时还向内收缩，最终稳定在约0.7 倍旋翼半径位置处。

图11 悬停状态的桨尖和桨根涡线Fig. 11 Root and tip vortex lines of the blade in hover

加密准则是一旦尾迹线穿过网格块，则定义该网格块为待加密网格块，直到网格密度小于目标网格密度为止。这里定义目标网格密度为桨尖弦长的2%。图12 为悬停状态下的自适应网格分布，可以看到在分布有尾迹线的地方，网格均进行了局部加密。由于模拟旋涡的网格尺度要小于桨叶子网格外边界的平均网格尺度，因此，当地笛卡尔网格的分布更密更集中，远小于重叠区附近的网格尺度。

图12 悬停状态的自适应网格加密Fig. 12 Adaptive mesh refinement of the rotor in hover

2.3.2 旋翼前飞状态

前飞模拟状态为UH-60A 旋翼的C8534 前飞状态，该状态对应UH-60A 旋翼的大速度中等过载飞行状态：其中拉力系数CT=0.0069， 桨尖马赫数Mtip=0.642 ，前进比 μ=0.368， 桨盘倾角αs=-7.31°。

旋翼前飞状态尾迹用Beddoes 预定尾迹模型[16]表示，该模型给出了尾迹分布的代数表达式，为研究旋翼尾迹分布以及桨涡干扰BVI 现象提供了方便。Beddoes 预定尾迹模型中包含了尾龄角 ψw和桨叶方位角 ψb， μ为前进比。涡元的轴向位移z则由旋翼轴向速度和局部诱导下洗速度积分组成，则桨尖涡尾迹几何形状可表示为：

其中，xtip、ytip、ztip为桨尖涡尾迹的笛卡尔坐标，用旋翼半径R无量纲化，坐标原点位于旋翼中心。rv为控制桨尖涡卷起及随涡龄角向内收缩的径向位置参数， μ为前进比，α为来流迎角， λi为当地诱导入流比，计算公式如下：

其中： λ0为 平均诱导入流比，E为尾迹倾斜角度参数。

图13 给出的是UH-60A 旋翼前飞状态下的桨尖尾迹分布三维视图，由Beddoes 预定尾迹模型计算得到。可以看到在前飞状态下，旋翼尾迹呈螺旋状朝旋翼后下方发展。研究表明，旋翼桨尖涡的涡核直径为桨尖弦长的1/10 量级，模拟桨尖涡所需网格的长度尺度约为桨尖弦长的1/100 量级，并且桨尖涡的分布范围很广，给传统数值模拟方法的网格生成带来了极大挑战——采用近场全局均匀加密方法导致网格量不可承受，必须发展网格当地加密的自适应生成技术。

图13 前飞状态的桨尖涡线Fig. 13 Tip vortex lines of the rotor in forward flight

网格最大加密层数设置为9 层，由于网格加密的密度要求高，尾迹线分布范围广，最终网格块数达到了4.4×104块，网格单元数达到2.25×107。加密运算以单个进程串行方式完成，在主频3.2G 的Intel I7-8700 CPU 芯片运算网格加密时间约为76 s。图14是旋翼前飞状态下的自适应网格分布，可以看到当前网格加密方法适应了尾迹线的分布情况，在旋翼的后方和下方区域对网格进行了加密，在桨尖尾迹线穿过的地方网格密度达到最大值。

图14 前飞状态的自适应网格加密Fig. 14 Adaptive mesh refinement of the rotor in forward flight

图15 用分层方式显示了旋翼附近的网格，位于中间的第6～7 层网格是在前5 层网格上的叠加，位于右边的第8～9 层网格又是在前7 层网格基础上的叠加；最大网格密度位于桨尖涡尾迹附近，为尾迹流动的精细模拟创造了条件。

图15 自适应网格的多层结构Fig. 15 Multi-layer structure of adaptive mesh refinement

图16 给出了旋翼前飞状态自适应网格的第8 层网格分布，网格最大层级为9 层，第8、第9 层代表了网格的核心加密区，可以看到这些密网格分布在旋翼附近及旋翼后下方区域。图中不同颜色代表网格块所在的不同分区，分区方法利用空间填充Z 曲线的聚类特性，使分在同一个进程的网格块保持相邻。同时需要注意的是由于空间填充Z 曲线在结束一个Z 遍历以后有跳跃的特点，前一个Z 的最后一个网格块与后一个Z 的第一个网格块并不相邻（见图5（d）），因此可能出现分在同一个进程的网格块不相邻的情况，有研究[17]表明在同一个进程内，不相邻的网格块最多为两组，这样保证了网格块不过于分散，降低了并行发送与接收方面的需求。

图16 块结构笛卡尔网格的并行分区Fig. 16 Parallel partition of the block structured Cartesian grids

3 结束语

本文探讨了旋翼计算中背景笛卡尔网格的生成方法，发展了块结构化笛卡尔网格的生成技术。将该技术应用在UH-60A 旋翼上的悬停和前飞流场，获得了网格密度分布合理的自适应网格，主要结论如下：

1）采用双网格建模思路，背景笛卡尔网格由于不需要适应物面边界，网格划分过程十分快捷高效。背景网格可以通过几何特征或流场特征对网格进行自适应加密，使当地网格密度满足重叠插值或流场计算的要求。

2）块结构自适应笛卡尔网格采用多层网格结构，高效实现了网格的加密和稀疏操作。用八叉数数据结构行管理网格块，相对于传统笛卡尔网格管理网格单元而言，整个八叉数所管理的节点数大大降低，提高了整个数据结构的管理效率。

3）采用多层空间填充Z 曲线遍历网格，利用该顺序对网格块进行编号和并行分区，分区方法利用空间填充Z 曲线的聚类特性，使分在同一个进程的网格块尽可能保持相邻，降低了并行发送与接收方面的需求。