网格生成技术

1. 网格生成技术数值模拟流体运动时，首先将流动区域离散成一定形状的网格，然后在网格节点上求解离散化的控制方程。

数值模拟的计算精度既与控制方程的离散格式的精度密切相关，也与网格结构和分布有关，为了尽量减少计算误差，保证解的稳定性，生成的网格至少满足下面的一些原则：①网格的贴体性。

计算网格应准确反映流动区域的边界形状，并且要能较容易的引入边界条件。

②网格的疏密变化。

在物理梯度大的地方，网格要密些，以提高对流动结构的捕捉能力，搞高数值模拟的精度。

另外，由于在数值模拟之前，人们对流动结构的主要特征不甚了然，如哪此地方会出现旋涡，水跃、激波会产生在什么地方等，因此，计算网格最好能根据计算结果的变化而进行调整，即网格具有自适应性。

由于流动边界和流体运动结构的复杂性，自动生成复杂流场的理想分布网格相当困难，网格生成所费人力也很高，即使在计算流体力学高度发展的国家，网格生成仍占一个计算任务全部人力时间的60%～80%，因此，网格生成技术成为了CFD（计算流体力学）中一个独立的分支，网格生成技术也是CFD中最活跃的研究领域之一。

目前，网格生成方法很多，根据网格拓扑结构可分为两大类：即结构网格和非结构网格。

1.1结构网格的生成结构网格中网格节点与邻点相连，连方式与节点的位置无关，如二维空间中的矩形网格、三维空间中的六面体网格。

对于简单的计算区域，很容易进行结构网格的剖分，对于复杂的流动区域，尽管可以采用阶梯形网格对边界进行近似处理，但是这种处理通用性差，且会影响计算精度，为解决这个问题，人们一般采用坐标变换技术生成计算网格。

坐标变换生成计算网格又称贴体坐标技术，其基本思想是通过数学变换将复杂的物理区域变换到规则的计算空间中，物理空间和计算空间一一对应。

目前生成贴体拟合坐标的方法可以分成代数变换和偏微分方程变换两大类。

代数变换网格生成是用代数公式，一般为显示，给出物理区域和计算区域之间的对应关系，常用的方法有保角变换（conformal mapping）、剪切变换(shearing transformation)和Hermit变换等。

代数变换网格生成方法应用范围有限，其原因是对于复杂的计算区域，代数变换较难实现，边界附近的节点控制也十分困难。

偏微分方程方法用微分方程将不规则区域变换成规则区域，其通用性较好，又有生成的网格均匀、网格疏密易于控制等优点，由此得到了普遍的应用。

Winslow于1967年提出用偏微分方程生成计算网格的思想，后来，Thompson，Thomas和Mastin对这一方法进行了全面而系统的研究，提出了著名的TTM方法。

在TTM方法中，计算网格控制方程中源项的各控制参数的选取没有一定之规，具体参数的选取与研究者的经验很有关系。

由于控制方程中源项影响网格的疏密分布，所以，源项中各控制参数的选取也成为TTM方法中的一个弱点。

Thomas 和Middlecoff于1980年提出了一种确定控制方程中源项的方法，该方法可以根据边界节点信息控制网格的疏密，但是求出的网格的正交性不太理想。

Hilgenstock提出了求控制方程源项的新思路，即根据网格与边界线的夹角来确定源项，大大改善生成网格正交性能。

现在Hilgenstock方法已在网格生成技术中得到了比较广泛的应用。

对于较为复杂的计算区域，在全域中生成一套计算网格比较困难，因为，此时局部网格会发生严重畸变，生成的网格不再适用于数值求解。

为此，人们提出了分区的思想，即将计算区域分解成若干几何边界比较简单的子区，在每个子区独立生成计算网格，同时也允许必要时在各子区使用不同的物理控制方程。

分区方法生成计算网格的优点是：⑴从局部到整体保证所生成的计算网格具有较高的品质，如较好的光顺性和正交性；⑵便于控制各子区的网格疏密；⑶易于实施并行计算。

从子区的连接方式看，分块网格可以分成两大类，即多块嵌套网格和多块搭接网格。

多块嵌套网格中，子区相互重叠，重叠区的网格线可以一致，也可以不一致。

如果重叠区的网格线不一致，则计算时，应用插值方法传递块与块之间的信息。

在多块搭接网格中，子区互不重叠，相邻块之间以公共边界相连接，块与块之间的物理量通过边界点传递，无需插值，边界上的动量自然守恒，在相同的条件下，其网格数量少于嵌套网格，但对于不可压缩流体运动问题，当用SIMPLE(或SIMPLER)方法求解流场时，要特别注意边界上压场的处理，如果处理不当，会在公共边界上形成压强间断，进而影响流场的求解。

实际网格生成时，这两种方法可以混合使用。

目前，多块网格生成技术已在空气动力学中得到了广泛的应用，尤其普遍应用于航空器周围绕流场的计算中，并出现了许多通用程序。

1.2非结构网格生成技术非结构网格中网格节点的连接方式与位置有关，每个节点的邻点数量与节点附近的网格结构有关，其单元形状一般为三角形（二维）或四面体（三维）。

非结构网格具有很强的几何灵活性，同时其数据结构的随意性有利于实现自适应策略，因此，近年来，非结构网格生成技术的研究和应用在CFD中很受重视，公开发表的非结构网格生成的文献十分丰富，关于非结构网格生成技术的学术讨论也时有召开。

根据网格生成思路不同，非结构网格生成方法可以分成：Delaunay三角形化方法、阵面推进方法、四分树（八分树）方法等几种，其中应用最为普遍的是Delaunay三角形化方法和阵面推进方法。

1.2.1 Delaunay三角形化方法Dirichlet(1850)研究发现，在平面上给定点集P，则s P可将此平面s分成互不重叠的网格，每个网格内含有给定点集中的一个点P，且网i格内的任意一点P到P的距离比到点集s P中的其它任何点j P都近，该i网格可以用集合j V 表示，即 },,,:{}{s j i j P j i i j P P P P P V ∈≠∀-≤-=网格所围的区域称为V oronoi 区域，把相邻的V oronoi 区域的形成点i P 相连时，则形成Delaunay 三角形。

以平面问题为例，如图所示，4321P P P P 为给定点集，点集将整个平面分为四个V oronoi 区域，即4321A A A A 、125A A A 、5236A A A A 和634A A A ，四个点连成了两个三角形321P P P 和431P P P 。

V oronoi 单元及对应的Delaunay 三角形有许多重要的特性，如i ）Delaunay 三角形的外接圆中只包含对应的V oronoi 顶点（In-circle property ）；ii ）给定点集，在所有可能形成的三角形中，Delaunay 三角形的最小内角大于所有可能形成的三角形中的最小内角（Minimum-Maximum property ）；iii ）给定点集，生成的Delaunay 三角形是唯一的。

Delaunay 三角形的这些特性，尤其是特性i ）已成为用Delaunay 三角形化方法生成非结构网格的主要判断准则。

A3A2A1A5A4A6P1P2P3P4图1. Voronoi单元示意图V oronoi 子区域可以用细胞生长方式类比，即假定空间（或平面）中有N 个细胞（点集），它们以相同的速度生长，当细胞壁相互接触时，细胞停止生长，最终形成的图案就是V oronoi图形。

V oronoi多边形又可以用动物的活动领地来解释，将离散的节点看成离散的食物源，则V oronoi多边形就是食物附近动物的活动领地。

所以，八十年代以前V oronoi多边形主要应用于气象学（如大气压或海拔高度的拟合）、社会统计学（如动物在领地中的行为）等领域，八十年代之后V oronoi多边形及相应的Delaunay三角形才被广泛地应用于计算流体力学之中。

Delaunay三角形化有许多算法，常用的有如下几种：⑴寻找V oronoi邻居，如Kennon方法；⑵Bowyer或Waston算法；⑶边交换方法（edge swap）; ⑷最长边传播路径方法（the Longest-Edge Propagation Path）。

寻找V oronoi邻居算法的基本思路是：根据Delaunay准则，即Delaunay特性i），对计算区域中的每一个规定节点寻找V oronoi邻居，然后将互为V oronoi邻居的点用直线相连，即形成了三角形剖分。

这种网格生成方法需对计算区域中的每一个规定点进行V oronoi邻居判断，搜索运算的工作量相当大，计算时间较长，很难推广到三维空间的网格剖分中。

Bowyer/Waston算法是Bowyer和Waston于1981年分别提出的，该方法以不断地加入节点的方式进行，具体思路如下：假设已对N个节点实现了三角形化，现增加一个节点，由于该节点的加入，原V oronoi图案遭到了破坏，必须进行修改，因此，首先判断哪些三角形不满足Delaunay性质，删除这些单元，并根据Delaunay 特性连接被删除单元的顶点和新增节点，即实现了新的三角形化。

Bowyer 算法所需时间可以用式N b N N a k K +log 表示，式中，N 是生成的节点数，k k b a ,是与空间维数k 有关的常量，该算法可以从数学上证明上最优的，并且很容易推广至三维空间中，因此在非结构网格的生成中得到了广泛的应用。

边交换方法主要用于网格的局部优化中，因此经常与其它方法一起使用。

最长边路径传播方法的思路是：根据三角形的最长边对区域中的三角形进行排序，如排成序列n t t t t ......210，其中三角形i t 是三角形1-i t 最长边的邻居，然后在所有三角形最长边的中心根据Delaunay 准则插入新点，重新生成网格，直到所有三角形满足给定要求（如所有三角形内角大于30°）为止。

从网格的生成方法易知，该算法的判断量（哪些三角形是相邻单元）计算量（需计算三角形的边长）是很大的，因此，将该方法推广到三维情形也有一定的困难。

上述四种算法均以Delaunay 准则为基础，但这些基本算法（算法4除外）有两个问题没有涉及，其一是边界的完整性如何得到保证，其二是网格的疏密如何控制。

为此，国外许多学者做了大量的工作，提出了许多具有实用价值的方法。

Weatherill 等提出，先在计算区域中插入边界节点，然后根据网格节点的分布函数插入内部节点，网格节点的分布函数由边界节点的分布信息确定或事先给定，也可以根据流场的计算信息确定（即采用网格自适应策略）。

边界节点插入后要进行边界完整性检查，以确保边界的贴体性要求。

Weatherill 方法充分利用了Bowyer算法的优点，是一种成熟高效的方法。

Rebay 也是先插入边界节点，然后在域内生成网格。

在生成域内部节点时，他提出了两种方案，一、对已生成的三角形根据节点分布函数进行排序，在节点分布函数最大的三角形外接圆心处插入新点，直到节点分布函数满足要求为止；二、对三角形进行分类，如下：I域内三角形II域外三角形i 符合要求的三角形ii 不符合要求的三角形a活动三角形b非活动三角形活动三角形指的是与其相邻的至少有一个满足要求的三角形或域外三角形，非活动三角形则反之。