11.2 非结构网格

11.2 非结构网格

非结构网格通常由二维三角形或三维四面体组成。但如今，用多种单元类型构造非结构网格越来越常见。例如，以六面体或棱柱体离散边界层，流动区域其余部分用四面体填充，并在六面体或棱柱体与四面体之间用金字塔形单元过渡，故称混合单元网格。结构六面体网格的优点是：对于高拉伸黏性网格，能保持壁面法向精度，而且与四面体网格相比，单元、边和面数量更少。非结构四面体网格的优点则是：能快速离散图11.3、11.4这样的复杂几何，且只需很少的用户干预。混合网格力求结合两者的优点。

非结构网格的节点与单元近似随机排列，即无法通过编号直接识别相邻单元或网格点，见图3.3。这使其具有极大的几何灵活性，因为网格无需遵循预先规定的拓扑。此外，根据物理解自适应加密或粗化网格，非结构网格也比结构网格容易得多。

用于 CFD 的非结构网格生成方法，通常基于以下两种方法之一：

• Delaunay 方法；

• 波前推进方法。

这两种方法也可结合使用。根据基础方法的不同，分别称为波前推进 Delaunay 格式 $ [43–45] $，或波前 Delaunay 格式 $ [46–48] $。除这两种标准技术外，还有所谓的气泡堆积算法等方法 $ [49–54] $。这里仅介绍 Delaunay 方法与波前推进方法的基本特征；两者的综述例如见文献 $ [20, 55, 56] $。

Delaunay 方法主要指连接网格点的一种特定方式，在二维形成三角形，在三维形成四面体。Delaunay 剖分最重要的特征是：任意三角形的外接圆，或任意四面体的外接球，内部均不包含其他网格点。空外接圆准则意味着，二维情况下所有三角形的最小角得到最大化，即最大最小角剖分，因此可获得较好的网格规则性。

波前推进方法的基本思路是，从已知的边界离散，即表面网格出发，逐个单元地生成网格。单元的开放表面构成波前；通过在波前前方布点，构造新三角形或四面体。这样，波前逐步穿过计算域，直至所有空腔均被填满。布点由所谓的背景网格控制。波前推进方法的优点是点分布平滑，并自动保证边界完整性，但速度比 Delaunay 方法慢。

11.2.1 Delaunay 剖分

Delaunay 剖分基于 Dirichlet 于1850年提出的方法，将空间唯一划分为一组紧密覆盖的凸区域 $ [57] $。对于给定点集，每个区域包含某个点周围的空间，其中任意位置到该点的距离，均小于到其他点的距离。这些区域形成多边形，三维则形成多面体，称为 Dirichlet 镶嵌或 Voronoi 图 $ [58] $。若用直线连接共享 Voronoi 图某条边或某个面的区域所对应的点对，即得到 Delaunay 剖分 $ [59] $。所得三角形或三维四面体集合覆盖点集的凸包，见图11.13。Delaunay 剖分是 Voronoi 图的对偶。二维 Voronoi 多边形的顶点是相应三角形的外接圆圆心；三维则是四面体的外接球球心。因此，每个三角形的外接圆或每个四面体的外接球，内部均不包含点集中的其他点。

如前所述，Delaunay 方法规定了网格点的连接方式，点的位置必须由其他技术确定。因此，构造 Delaunay 网格的一种常用方法，是向初始剖分中逐个插入节点，再局部重新剖分，使其满足空外接圆或空外接球准则。增量插点策略可描述为以下步骤：

1. 离散物理区域的边界。

2. 生成包含全部边界节点的初始 Delaunay 网格：可以直接剖分边界节点，如图11.14a；也可以构造一个包围区域的四边形，三维为六面体，再划分为三角形或四面体，如图11.14b。

3. 若尚未完成，则采用满足 Delaunay 条件的技术，将全部边界节点插入初始剖分。

4. 将不满足某种尺寸或质量指标的全部三角形或四面体列成清单，并排序，使最差的单元排在最前。

5. 在清单中第一个单元的外接圆心或外接球心 $ [60-62] $（见图11.15），或在 Voronoi 图的边上 $ [63, 64] $，放置一个新点，然后局部重新剖分。检查每个新单元，若不符合尺寸或质量指标，则加入清单。

6. 若清单中仍有单元，返回步骤5。

7. 删除区域外的单元，并恢复边界。

8. 检查网格质量，见11.2.5节；必要时光顺网格及／或交换边。

也可在步骤3之后，先完成步骤7的边界恢复。

可见，Delaunay 网格生成的核心任务，是向有效剖分中插入新点，即步骤3和5。为此已提出多种算法，最常用的分别由 Green 和 Sibson $ [65] $、Bowyer $ [66] $、Watson $ [67] $ 提出。下面介绍 Watson 算法。

为完成步骤4和5，需要采用适当指标评估单元，可采用单元质量，例如最小角或长宽比，见文献 $ [68, 69] $，也可采用单元尺寸，例如体积或边长。单元质量可直接根据几何形状评估，但尺寸指标必须定义为计算域内空间位置的函数，可采用以下方法：

• 指定解析函数，例如尺寸与距物体的距离成正比；

• 利用初始剖分，将边界上的尺寸分布插值到内部，例如见文献 [20，第1章，第17—20页]；

• 在基于四叉树、三维为八叉树结构的背景网格上插值 $ [70-75] $；

• 在区域内部指定点源、线源等，例如见文献 $ ^{[20]} $ 第1章第20—22页。

此外，也可借助波前推进技术布点 $ [45, 46] $。

Watson 算法

Watson 插点与重新剖分方法包括以下步骤 $ [67] $：

1. 定位包含待插入点 P 的单元，见图11.15。

2. 找出外接圆或三维外接球包含点 P 的全部单元，见图11.16左侧。

3. 从剖分中删除所有上述单元。

4. 将凸空腔边界上的点与点 P 连接，形成新单元，见图11.16右侧。

若为每个单元存储以下信息，就可得到特别适合 Watson 算法的数据结构：

• 构成单元的节点编号；

• 指向与该单元共享面的相邻单元的指针；

• 外接圆或外接球的中心与半径。

采用这种数据结构，可高效完成步骤2的相关单元搜索。首先检查包含插入点的单元的邻居，再检查这些邻居的邻居，依此类推。若某个单元的外接圆或外接球不包含插入点，就停止沿该方向搜索。Delaunay 剖分的几何性质保证，该单元后面的相邻单元也不包含插入点；同时保证这种简单策略能定位所有相关单元 $ [67, 76] $。Watson 算法的计算量阶数为 $ N \log(N) $，其中 N 为网格点数，因此网格生成速度很快。

约束 Delaunay 剖分

Delaunay 剖分不会自动保留指定的边和面，例如物理区域边界上的边和面。所谓约束 Delaunay 剖分正是为此而提出 $ [77] $。图11.17展示二维边界边的恢复，根据具体情况，需要交换边或重新剖分。二维约束 Delaunay 剖分始终能产生有效网格，但三维情况无法保证。三维边界离散的恢复必须分两步：先恢复边界边，再恢复边界面 $ [56] $。某些情况下，必须在边界上插入额外的 Steiner 点 $ [78–81] $，否则空腔无法重新进行四面体剖分。恢复边界边及边界面后，即可删除流动区域之外的单元。

11.2.2 波前推进方法

波前推进方法由 Peraire 等人 [82,83] 及 Löhner 和 Parikh [84] 于20世纪80年代末首先提出。网格生成过程可概括为以下步骤：

1. 离散物理区域的边界，即生成表面网格。

2. 生成表示波前的边清单，三维则为面清单；初始波前为边界边或面。按边或面尺寸由小到大排序 [55]，这有助于生成尺寸平滑变化的单元。

3. 选择清单中的第一条边或第一个面，在其中心上方沿法向放置新点 $ P_{0} $，见图11.18。向计算域内部推进的距离 d，由尺寸分布函数的当地值决定，例如见文献 [56]。

4. 以点 $ P_{0} $ 为中心，构造半径 r 的圆，三维为球；r 取决于当地网格尺寸。

5. 找出圆或球内部的全部点。

6. 若其中没有已有点，就用点 $ P_{0} $ 生成新单元；否则，按这些点到 $ P_{0} $ 的距离排序，即 $ P_{1} $、$ P_{2} $、$ P_{3} $。依次用各点构造单元，接受第一个不与其他单元相交、且满足指定质量指标的单元。

7. 删除当前波前边或面，将新形成的边或面加入清单，然后重新排序。

8. 返回步骤3，直至清单为空。

9. 检查网格质量，见11.2.5节；必要时光顺网格及／或交换边。

图11.19a—d 给出一个二维构型在波前推进过程不同阶段的网格。

区域内的单元尺寸分布，主要由边界上的点密度决定。此外，可在区域内部布置点、线、圆柱等不同类型的源，控制步骤3中的距离 d，例如见文献 [20，第1章，第20—22页]。可从基于四叉树、三维为八叉树数据结构的背景网格获得当地单元尺寸 $ [70-75] $。也可通过边界节点的 Delaunay 剖分插值尺寸 $ [46-48, 85] $，或采用笛卡尔背景网格 $ [86] $。四叉树或八叉树数据结构还可用于高效搜索附近点，即步骤5，并检查单元是否相交。文献 [3,87] 提出基于弹簧类比的简单波前相交检验。不同搜索算法的详细比较见文献 [88]。

与 Delaunay 剖分相比，波前推进方法能够自动保留边界离散，并更好地控制单元尺寸与网格平滑性。此外，Delaunay 方法对舍入误差更敏感。不过，点搜索及相交检验算法需要大量计算。CFD 中经常将波前推进与 Delaunay 方法结合，尤其在三维情况下 $ [43–48, 89] $。波前推进方法的较详细描述见文献 $ [90] $。

11.2.3 各向异性网格生成

11.2.1和11.2.2节介绍了生成各向同性三角形或四面体单元的两种方法，适用于模拟无黏流动。然而，精确求解 Reynolds 平均 Navier–Stokes 方程，需要在横跨边界层与尾迹的方向上具有很高的空间分辨率。因此，对于高 Reynolds 数流动，各向同性网格会产生过多单元。显然，黏性流动区域应采用各向异性、高长宽比的单元。可使用天然适于拉伸的四边形，三维为六面体或棱柱体，由此得到如今非常常用的混合单元网格，见下一小节。也可分解六面体或棱柱体，获得纯四面体网格；或者直接生成拉伸三角形、四面体网格。无论采用哪种方法，均应避免形成图11.20a这样的钝角三角形或四面体。文献 $ ^{[91]} $ 证明，这类单元会造成极大的空间离散截断误差，因此宜生成图11.20b、c这样的直角单元。不过，根据控制体类型的不同，直角单元也可能引起离散误差 $ [92, 93] $，参见5.2.3节。

拉伸 Delaunay 剖分

生成拉伸三角形或四面体单元，需要定义物理区域内的拉伸程度与方向。拉伸 Delaunay 剖分需要用背景网格提供这些信息。

可采用生成各向同性 Delaunay 网格时相同的插点技术，但空外接圆或空外接球准则，改为空外接椭圆或空外接椭球条件。椭圆或椭球沿当地拉伸矢量方向定向，其轴长比反映拉伸程度 $ [94-96] $。另一种方法是采用气泡堆积算法，并以椭球替代球 $ [51-53] $。

层推进方法

对于几何复杂的区域，各向异性 Delaunay 剖分的拉伸矢量难以指定。因此，Pirzadeh $ [3, 4, 97, 98] $ 及其他研究者 $ [99, 100] $ 提出的层推进方法更为常用。它与双曲型结构网格生成技术（11.1.4节）相似，也可视为经过修改的波前推进方法。

采用层推进方法生成拉伸三角形或四面体网格，包括以下步骤，见图11.21：

1. 对全部边界表面进行三角剖分。

2. 计算边界节点处的近似法向矢量。

3. 沿表面法向布置网格点，生成厚度逐层增加的四边形单元层，三维为棱柱体单元层。

4. 将每个四边形划分为两个三角形，将每个棱柱体划分为三个四面体。三维连接模式需要特别注意。

5. 持续推进各单元列，直至波前与自身或其他波前相交，或单元长宽比接近1。

随后，以各向同性四面体填充流动区域的其余部分。通常采用波前推进方法，从最后一层边界单元构成的表面出发。图5.2和11.4给出按照这一流程生成的网格示例。

原则上，计算步骤2的边界节点法向矢量有两种选择。首先，若边界表面采用 NURBS（非均匀有理 B 样条）等解析函数描述 $ [101] $，则容易计算法向导数。文献 $ [102] $ 介绍了这类网格生成方法，但必须注意表面片接口，例如后缘处。

第二种、更常用的方法是利用表面三角剖分。二维情况下，只需对同一节点处相接的边界边法向求平均，因此比较容易。但三维情况下，这种简单平均并不合适，因为边界节点可由任意数量的三角形面共享，法向矢量会因剖分不同而产生偏差。在尖角和转折处，问题尤其严重。最优法向矢量的必要条件由所谓可见性准则给出 $ [103, 104] $：沿法向布置的一个点，必须对共享该边界节点的所有三角形面同等可见。已有多种方法据此构造法向矢量，例如见文献 $ [3, 103–106] $。为防止推进方向突变及网格线交叉，需要光顺法向矢量。通常采用加权 Laplace 型光顺 $ [104, 105] $，即

$$ \vec{n}_{i}=\omega\vec{n}_{i}^{(0)}+\frac{(1-\omega)}{\sum_{j}1/|\vec{r}_{ij}|}\sum_{j}\frac{\vec{n}_{j}}{|\vec{r}_{ij}|}.\tag{11.12} $$

式（11.12）中，$ \vec{n}_{i} $ 与 $ \vec{n}_{i}^{(0)} $ 分别表示新的和初始的节点法向，$ \vec{n}_{j} $ 表示相邻节点的法向矢量，$ |\vec{r}_{ij}| $ 为节点 i 与 j 之间的距离。权重因子 $ \omega $ 取决于表面曲率，凹区域取较小值，凸区域取较大值。每层分别进行光顺。通常，最初几层采用初始节点法向，之后逐步增加光顺程度，因为边界正交网格有助于降低空间格式的离散误差。

推进步长依据用户指定的首层厚度及法向拉伸率计算。例如，间距可按文献 $ ^{[3]} $ 确定：

$$ \Delta n_{k}=\Delta n_{0}\Big[1+a(1+b)^{k-1}\Big]^{k-1},\tag{11.13} $$

其中，$ \Delta n_{k} $ 表示第 k 层间距，$ \Delta n_{0} $ 是指定的首层厚度，$ 0.04 \leq a \leq 0.2 $ 和 $ 0 \leq b \leq 0.07 $ 为拉伸参数。还可在凸角或凹角处进一步修正推进步长 $ \Delta n_{k} $ [106]。也可沿边界层增长方向增大 $ \Delta n_{k} $，使 $ y^{+} $ 保持常数。

层推进算法的每一级，在二维生成一层四边形，在三维生成一层棱柱体；若表面采用四边形离散，则生成六面体。四边形可沿较短对角线容易地划分为两个三角形，但将棱柱体划分为三个四面体则需要谨慎：相邻棱柱体的公共面必须采用相同的划分方式。文献 $ [3] $ 提出一种基于整数编号的一致棱柱体分解方法；文献 $ [107] $ 给出更简单的格式，另一种分解方法见文献 $ [108] $。

法向推进与 Delaunay 的组合方法

生成拉伸三角形或四面体网格的另一种方法，是将改进的波前推进布点方式，此处称为法向推进布点，与 Delaunay 剖分结合。Müller $ [109] $ 提出了二维方法，Marcum $ [110, 111] $ 及 Sharov 和 Nakahashi $ [106] $ 分别发展了三维算法。首先，对边界节点进行体剖分，将所得剖分用作背景网格，同时作为高效搜索结构。此时可以恢复边界 $ [106] $，但并非必需 $ [110, 111] $。随后，采用类似层推进的技术，在拉伸网格区域插入点，并利用背景剖分确定附近波前、检查相交。最后，以向外膨胀后的表面为新边界，采用标准 Delaunay 技术生成各向同性单元。

拉伸表面网格

在流动具有强方向性的区域，各向异性网格可大幅减少三角形面数及体单元数，同时保持求解精度，见图11.22。例如，RANS 模拟不需要在机翼或叶片展向保持很高的网格分辨率，机身轴向也如此。文献 [112,113] 报告，表面三角形数可降为原来的二分之一至六分之一。但是，如图11.22所示，必须正确确定表面单元的方向，尤其在高曲率区域；否则，真实表面轮廓表示不佳，会使流场解失真。拉伸表面网格可采用上述任一种各向异性网格生成方法；表面网格生成的一般讨论例如见文献 [16,75,114]。拉伸比与方向通常通过线源指定。另一种方式是先在参数坐标中生成各向同性网格，再用拉伸函数将网格变换到物理空间。

11.2.4 混合单元网格与混合拓扑网格

如今，用棱柱体、四面体、金字塔形单元等不同类型单元组成的网格离散流动区域，已经非常普遍。这类网格称为混合网格或混合单元网格。另一种思路是以结构区域和非结构区域共同组成网格 $ [115–122] $，称为混合拓扑网格（hybrid grid）。需要注意，有些作者虽使用 hybrid grid 一词，实际指的却是混合单元网格。采用两类网格的动机，都是在各流动区域使用最适合的单元或拓扑，以提高模拟精度、降低计算成本。不过，网格生成过程必须尽可能自动化，而且足够快速。

混合拓扑网格在实践中相对少见，因为需要耦合一个结构流动求解器与一个非结构流动求解器。编写两个 CFD 程序，尤其持续更新与维护它们，需要投入大量工作，因此必须设计专门的数据和程序结构。混合拓扑网格的典型应用，是以结构网格离散近壁区域，采用六面体，或采用半结构棱柱体网格，以便实现高阶空间格式、隐式方法 $ [123] $ 或多重网格 $ [124, 125] $；其余区域由非结构网格填充，从而在复杂几何区域获得明显优势。混合拓扑网格也可用于叶轮机械中的转子—静子相互作用：转子与静子周围结构网格之间的接口，由一层非结构网格组成，类似文献 $ [126] $ 的思路。转子每次移动后，重新生成该非结构网格，从而在接口处保持离散格式的精度与守恒性。

与混合拓扑网格相比，混合单元网格可能更灵活，尤其适用于复杂几何；它也更容易生成，特别是更容易加密。但是，流动求解器必须采用能统一、高效处理不同单元类型的数据结构和数值格式。第5章已经讨论了部分相关问题。下面简要介绍棱柱体／四面体，以及棱柱体／四面体／笛卡尔混合网格的生成方法。

棱柱体／四面体混合网格

Nakahashi $ [127-129] $ 最先提出一种网格生成方法，在近壁区域采用棱柱体，其余区域采用四面体。Kallinderis $ [130, 131] $ 以及 Connell 和 Braaten $ [107] $ 进一步研究了该方法。棱柱体单元的生成，采用前述层推进或法向推进的相同技术；四面体直接从最后一层棱柱体向外生成。若某个棱柱体的四边形面仍暴露在外，则采用金字塔形单元作为向四面体过渡的单元。

棱柱体／四面体／笛卡尔混合网格

笛卡尔网格由与笛卡尔坐标轴对齐的正方形组成，三维则由立方体组成。即使几何复杂，也能容易且低成本地生成笛卡尔网格。但对于弯曲或不沿坐标轴定向的固体边界，笛卡尔方法的不足是流场解精度较差。在需要高拉伸、边界正交单元的 RANS 模拟中，这一问题尤其严重。因此，Melton 等人 [132]、Karman [133]、Smith 和 Leschziner [134]、Wang 等人 [135]、Delanaye 等人 [136]、Hitzel 等人 [137] 等作者，提出在壁面附近插入贴体四边形或棱柱体。采用金字塔形单元和四面体作为过渡单元，可在内部网格层与外部笛卡尔网格之间建立连续接口。另一方面，也可允许悬挂节点与悬挂边，以简化网格生成流程，见图11.23。只要处理正确，接口处的通量交换就能保持守恒 $ [136, 138] $。一种混合笛卡尔网格方法的详细描述见文献 $ [139] $。Karman $ [140] $ 提出一个有意思的思路：将笛卡尔体网格投影到边界，自动生成表面离散，再沿固体边界插入黏性网格层。

11.2.5 网格质量评估与改善

网格质量强烈影响模拟精度，对非结构网格尤为如此。第5章多处提及这一问题，例如5.3.3节的解重构。因此，所得网格中的单元应尽可能规则，单元尺寸及黏性区域的拉伸程度应在计算域内平滑变化。当然，网格还必须足够精细，以解析相关的几何和流动特征。

图11.24展示网格中应避免的四面体形状：钝角单元，亦见图11.20a；薄片单元，即四个点近乎共面；针状单元；楔形单元。唯一的例外是黏性区域内的楔形拉伸四面体，见图11.20。只有使用棱柱体才能避免这类单元。但由于最大流动梯度沿壁面法向，即楔形短边方向，数值误差不会造成过大影响。

（a）钝角单元

薄片单元

(c)

针状单元

（d）楔形单元

除单元角度 [111] 外，四面体网格还可采用以下参数作为质量指标 [20,141]；所列数值对应正四面体：

1. 外接球半径／内切球半径 = 3.0；

2. 最大边长／内切球半径 = 4.8990；

3. 外接球半径／最大边长 = 0.6125；

4. 最大边长／最小边长 = 1.0；

5. $\bar l^3/V=8.4797$，其中 $\bar l$ 为单元平均边长，$V$ 为体积；

6. $V^4/(\sum_f S_f)^3=4.585\times10^{-4}$，其中 $S_f$ 为各三角形面的面积。

校注：原书的正四面体指标数值存在疑点。第3项精确值为 $\sqrt6/4\approx0.612372$；第5项精确值为 $6\sqrt2\approx8.485281$。第6项原式具有长度六次方的量纲，因此不能作为尺度无关的形状指标；指数与数值均需另行核对。译文保留原书数值与指数。

为改善网格质量，可采用 Lawson 提出的交换边算法 $ [142, 143] $。交换边尤其适合消除网格中的薄片单元 $ [48] $。

网格光顺

另一种常用于改善网格规则性的技术是光顺，通过近似 Laplace 算子移动网格节点：

$$ \vec{r}_{i}^{n+1}=\vec{r}_{i}^{n}+\frac{\omega}{N_{\mathrm{A}}}\sum_{j=1}^{N_{\mathrm{A}}}\big(\vec{r}_{j}-\vec{r}_{i}\big),\tag{11.14} $$

其中，$ N_{A} $ 为与节点 i 相邻的节点数量，$ \omega $ 为松弛因子：内部节点取0.5—1.0，紧邻边界的节点取0.25，边界节点取0.0，参见文献 [48]。式（11.14）须迭代求解。文献 [47] 优先采用点 Gauss–Seidel 格式，而非点 Jacobi 格式，因为前者更能抑制负体积。图11.25展示采用 Laplace 算子（11.14）光顺非结构网格的示例。