2.2 一般网格生成方法_河流数值模拟与信息化应用-QQ阅读男生轻小说网

上QQ阅读APP看本书，新人免费读10天

设备和账号都新为新人

2.2 一般网格生成方法

为利用不规则网格进行数值模拟计算，需建立一个能够描述计算域几何形状、物理属性和边界条件的计算模型。对于较简单的问题，输入的数据较少，所需要的输入数据能很方便地计算出来并手工处理然后用键盘输入。但对于中到大尺度模拟计算问题，手工生成网格数据就变得容易出错甚至不可行。而且在生成网格拓扑数据时，即使生成正确无误的模型，但若模型需要修改，如计算区域扩大或缩小，则需要从新生成模型，从而必然造成大量的重复劳动。对于二维问题，特别是非恒定流问题，必须以图形形式进行分析。因此，有大量的网格自动生成算法用于生成计算网格（Huyakorn&Dugoon，1975），这些算法一般都能编辑输入数据并自动生成整个或部分网格并包含如下功能：网格划分、局部网格优化及其引起的单元重新分布的优化、自动给单元和节点重新编码，并具有自动纠错功能。这些功能能够节省用于建立网格的大量工作并减少计算费用。网格生成程序一般能编辑所有的输入数据并自动生成所有或部分计算网格属性数据（拓扑关系数据、节点坐标与高程、单元平均高程等）。许多程序采用交互式操作方式并使用图形输入输出设备，如数字化仪、交互式绘图仪及图形终端，另外还可执行更多的任务以增加更大的灵活性和更有效的离散化，如：网格绘制、局部网格加密及连续的单元光滑以及对单元和节点进行编号，建立优化的拓扑关系。

数据处理模块的最重要部分为：自动生成所有或部分计算网格以及网格的编码及建立拓扑关系，这两种功能有可能大大减少建立网格的工作量并减少计算量。Haber（1981）估计，进行水流数值模拟时，80%的时间和费用花在采用传统手工方法获取数据、生成计算网格以至得到计算模型上。其估算即使存有误差，但由此也可知采用自动网格自动生成方法可大大节省工作量和计算费用。网格自动生成和编码算法各不相同，取决于处理程序和分析程序的类型。

一般地，网格生成分为两步：①将计算域剖分为形状相对简单的一个或多个计算子区域，该步可手工完成或利用计算机图形设备和软件交互式地完成（Haber，1981）；②将各个计算子区域再剖分为有顺序的单元和节点的集合，包括节点坐标、各单元的节点列表与其他一些相关信息，如结点地面高程和单元属性数据的获取。可用于第二步剖分的技术有：① 自动不规则三角网生成技术；②网格平滑技术；③坐标转换技术。这三种技术没有优劣之分，而且经常结合起来建立计算网格（万洪涛，2000）。

2.2.1 无结构网格

三角网自动生成技术可用于任何几何形状的计算区域，该类方法的主要缺点为生成的网格较差（shephard，1993），光滑技术常与自动三角网生成方法结合起来以修正所生成的三角网格。三角网自动生成算法可分为两类。

第一类三角网自动生成算法（Tracy，1976；Sadek，1980）统称为阵面推进法。该法第一步为在计算域边界适当位置上布设节点。这些节点组成边界多边形。通过对边界多边形进行操作，切去锐角或用多边形内部的点替代多边形上的点，逐步生成所需要的三角网。例如，Tracy（1976）提出的算法中，多边形的每个顶点都有一个朝向多边形内部的角，当该角小于90°时，则用线联结与其相邻的两个节点组成三角形，并且该点从多边形顶点中消除。对多边形中的所有锐角都如此处理，然后，对一个角度低于180°的顶点，参照该点及其两相邻节点的空间位置，在多边形内部适当位置放置一新节点，生成两个相邻的三角形。如果该步导致某些顶点的内部角度低于90°则采用第一步的方法生成新的三角形和多边形。以上操作一直持续执行，直到组成多边形的顶点只有三个，则该三个顶点为最后生成单元的顶点。

第二类三角网自动生成算法（Frederick and others，1970；Cavendish 1974；Bykat，1976）的第一步也为在边界多边形合适位置上与子域内部布设适当的节点。在子区域内部的所有节点既可用手工方法也可用自动方法设置。按照一定的规则将多边形边界和多边形内部的点相连以生成三角网，使生成的三角形应互不交叉设置且尽可能为等边三角形。

在非结构网格自动生成算法中，主要的方法有Delaunay方法、修正的四分树方法以及阵面推进法。

1.Delaunay三角形

在模拟复杂几何形状的流场时，网格的选取和生成是一个十分困难的问题。尽管有贴体网格生成的方法，但由于网格的安排是有序和按一定结构的，因此不可避免地会出现该密的地方不能做到很密，不该密的地方又变得很密，在采用自适应网格时更有这样的问题。另外在有结构网格中，有时要生成有一定次序的网格是非常困难的，以至于网格形状不能得到保证。无结构网格的生成，特别是三维情况，是十分耗机时的烦琐工作，寻找通用的高速有效的生成方法是急需解决的。

原则上讲，无结构网格可以具有任意结构，然而这样的结构实现起来是很困难的。而且生成的网格应该没有折叠，每个三角形的形状应当尽量接近正三角形，应当避免钝角接近180°或锐角接近0°的三角形出现，由于流场的区域往往是复连通区域，还必须避免网格落在域外。

图2.2 Delaunay三角形

1850年Dirichlet提出了一个思想：设平面上离散地分布着N个点P₁，P₂，P₃，…，P_N，选其中一个点P_i，则平面上离P_i最近的所有点构成一点集，形成一区域V_i，在V_i内所有的点P满足以下条件：

这里∀i≠j表示对于所有i≠j的情况，即j=1，2，…，i-1，i+1，…，N。这里V_i称为Voronoi区，相邻的Voronoi区的中心点P_i相连就形成三角形，它们称做Delaunay三角形。Dirichlet指出：这些Delaunay三角形是不会互相重叠的。在图2.2中有四个节点P₁、P₂、P₃、P₄、φ₁Q₂Q₁φ₂、φ₂Q₁φ₄、φ₄Q₁Q₂φ₃、φ₃Q₂φ₁为四个Voronoi区，而P₁P₂P₃及P₁P₃P₄即为两个Delaunay三角形，显而易见，Q₁、Q₂就是这两个三角形外接圆的圆心。

在生成上述三角形时需要注意两点：一是除给定的边界点以外，所有的点都是逐步加入的，因此需要给出加点的条件；二是由于三角形单元是不断地形成和消除的，因此为了使这一过程有比较高的效率，程序设计必须充分优化。加点的原则是：对于每一个点，P₀与其相连接的点为P_j（j=1，2，…，M），可以设一个数dP₀，即

三角形中心点的dP₀值定为所在三角形三个顶点的dP₀值的平均值，当该中心点到三顶点的距离d_m（m=1，2，3）满足d_m>αdP₀，且该中心点到三顶点的最短距离s满足s>βd_m时，则该点即为应加入的点。这时所在的三角形将破坏，这个三角形的相邻三角形，若它的外接圆包含这一中心点，则这两个三角形的公共边将被打破，从而重新形成新的三角形（见图2.3）。利用上述方法，可对三角形反复加点，直到不能再加点为止。这样形成的三角形不一定是很匀称的，但可以采用光滑方法对它们进行光滑。

图2.3 Delaunay不规则三角网的判定（Gable，1995）

2.四叉树-Delaunay方法

四叉树方法的基本思想是用树状的矩形网格覆盖计算域，然后将各矩形划分为三角形。采用这种方法生成的网格中，边界处的网格形状较差。前面已介绍过的Delaunay方法的优点是生成的网格质量较高，且效率较阵面推进法略高；但这种方法在区域外也有网格生成。阵面推进法的基本思想是首先生成边界的网格，然后逐渐向里推进形成空间网格；阵面推进是一个局部过程，因而执行过程可从边界点的任一点开始进行，新点的引入是在推进过程中自动生成的。

三角网的生成分为两步：区域内布点及其三角化。目前Delaunay方法已流行成为最常用的不规则网格自动生成方法之一。但严格说来，Delaunay并不是一个网格生成方法，它只不过是一个给空间结点提供连网一种的方法。因此采用Delaunay方法时必须先提供布点方法给研究区域产生结点分布，代表的算法是先产生边界点序列，边界点根据Delaunay准则产生初步网格，再插入其他点生成其他网格，然后对整个三角形网格进行优化从而使整个网格符合Delaunay准则。

四叉树布点方法是一种较好的布点方法，将四叉树布点方法和Delaunay连网方法相结合，汪承义（2000）提出了四叉树-Delaunay网格生成方法。其基本思想是运用四叉树方法对区域布点，从而使整个区域点的分布与控制条件的复杂程度相关联，再运用Delaunay方法对整个区域的点连网，从而得到质量高的网格，最后对整个网格进行光滑处理，并对三角形网格加入断线与边界控制。使整个网格变化平缓而更适宜计算。此方法无论计算域的边界如何复杂都可处理，并且生成的网格形状好、连通性好且变化平缓。方法的基本流程如下：

（1）运用四叉树布点方法进行布点。

（2）对产生的点进行数据分区。对数据进行分区是为了提高数据的查找速度，提高运算效率。一种快速的数据分区方法是将点集划分为N/K个相等的单元，每个单元中平均包含K个离散点（通常K=4），建立一个一维数组，存储每个分区单元中第一个点的索引，同时建立一个点链表，存储位于同一分区单元中所有点的信息以及其下一点的索引。对于要确定某个点落在三角网的某个三角形内，为了提高三角形检索的效率，需对三角形进行检索。三角形的检索方法是按形心坐标进行数据分区，这样要判断点P落在三角网的哪个三角形内，只需根据P点的平面坐标计算点P落在哪一个分区单元内，将该分区单元内的三角形取出逐一判断是否位于三角形内即可。

（3）根据边界控制产生初始三角网。

（4）把点集逐步插入到三角网。

（5）对产生的三角网进行优化。

（6）对优化后的三角网进行平滑处理。

（7）重复（5）、（6）过程三次。

（8）对得到的三角网插入特征断线。断线的处理主要分为以下几个步骤：①找出所有与断线相交的三角形并求出交点；②根据Delaunay点插入算法插入交点，并局部更新三角网的拓扑邻接关系；③更新断线，新的断线由一系列小断线组成；④重复①、②、③直到插入所有断线。

图2.4为断线处理的全过程，在图2.4（a）中，硬断线12与三角网中的三个三角形相交，交点为3、4。在图2.4（b）中，插入点3更新三角网。在图2.4（c）中，插入点4更新三角网，新的断线由13、34、42组成，从而得到了加入断线的三角网。

图2.4 断线的处理

3.光滑技术

光滑技术（Buell，1973；Herrmann，1976）为采用差分类型方法的迭代应用为有限元网格节点定位，通常与其他方法结合使用使生成的网格更为合理。

生成任意几何形状网格的最常用方法为Buell（1973）所描述的拉普拉斯方法，在该方法中，长方形网格的内部结点应满足方程

上二式中：（x_i，y_i）为内部结点i的坐标，（x_ij，y_ij）_j_=1，...，4为与其相邻的4个结点的坐标（见图2.5），该方法可推广于由任意四边形组成的非长方形网格中：

图2.5 网格结点及其邻结点

上二式中：N_i为与结点i相连的单元码；（x_nj，y_nj）和（x_nl，y_nl）为相邻单元n的结点坐标（见图2.5）。

因这些方程组成一个非线性代数方程组，他们通常用间接迭代方法求解，如果为大型计算网格则造成非常大的计算量，拉普拉斯变换生成的网格常不能拟合曲率大的边界且边界上的点距也不均匀，固定在边界曲率很大及边界节点距离不均匀的计算域上，用这种方法生成的网格很乱且看起来是斜的。

为了修正拉普拉斯方法对边界信息不敏感造成的不良后果，Herrmann（1976）提出了一种修正方法，该方法为基于8节点四边形元的局部二次多项式等参量变换，利用基于形状函数，8节点四边形的中心点坐标可通过以下公式计算：

上二式也可推广至非长方形网格中，则计算公式变为

以上式中变量的定义与前相似。这种公式称为等参量不规则网格生成方法。Herrmann利用以上四式进一步得到式中w为不小于0且不大于1 的权重，当w等于0时为拉普拉斯变换，当w等于1时为等参量网格生成方法。当0＜w＜1时，称为拉普拉斯等参量方法。虽然“纯粹的”等参量方法（w=1）能生成很好地拟合边界的计算网格，但其达到收敛的迭代数随w趋近于1而大大增大，因此现在网格质量和方法计算效率之间有所选择。

2.2.2 结构网格

结构网格要实现河道边界的贴体，一般需采用坐标变化技术。

坐标变换技术为单个多边形（二维空间中的单元三角形或正方形）与实际区域之间的坐标匹配，坐标变换技术为等参量制图的一个副产品——由Zienkiewicz和Phillips（1971）提出的拟合曲边有限元。采用这种方法整个的求解域粗略地分成为由等参量元组成的子域（见图2.6），利用多项式内插函数（形状函数），能得到特有的坐标变换图。则所希望的简单多边形分区可投影到实际空间以确定节点的Cartesian坐标。这种方法的主要缺陷在于：如果初始分区不够，其很难描述具有复杂边界的区域。

图2.6 等参量四边形网格

Gordon（1971）提出了一种采用混合函数（Blending Function）的制图方法来生成网格，该方法用区域内的表面来匹配区域边界的值，边界上的一套简单线性函数则“混合”起来在整个区域内形成更为复杂的内插函数，考虑到边界上的单元正方形，则映射函数φ被定义为

其中ξ和y为统一的参量坐标，内插φ得到精确拟合边界值的光滑表面。若函数中可描述四边形区域坐标，该函数则可用于转换实际区域内的所有网格。Barnhill等（1973）用近似正交的三角形网离散区域。Haber等（1981）声称线性超图（Transfinite Mapping）具有极好的特性，可用于网格生成，或利用这些制图函数以形成图形有限元预处理程序的基础。对于曲度很大的不规则区域，自动生成三角形网时，在区域边界以外也会有网格生成。Gorden和Hall（1973）引入了限制区域边界曲度的机理来解决这一问题，并以次数更高的制图函数用于拟合受限制曲线附近的网格线，另一种处理方法也将区域手工剖分成具有规则几何形状的小区域，得到的子域一般形状较规则；Haber等（1981）推荐采用该方法用于二维区域，而曲度限制技术可用于三维网格生成。

不少自动生成计算网格的坐标转换方法中，生成的计算网格往往不能正交。如美国J.F.Thompson等提出的曲线坐标变换法，是曲线坐标变换的一般表示方法，其ξ、η线的方向可任意布设，不能保证正交。正交曲线族出现于多种物理现象中，例如：有势流中的势函数和流函数、温度场中的等温线和热力线等。如果以这些物理现象为依托建立坐标网络，则可保证所建立的网格正交。

Willemse等根据势流理论中流函数与势函数必然正交的原则，导出了如下方程：