网格生成技术.doc

I
目 录
1 概述 1
2 结构网格 3
2.1 贴体坐标法 3
2.2 块结构化网格 11
3 非结构网格 16
3.1 概述 16
3.2 阵面推进法 16
3.3 Delaunay三角划分 19
3.4 四叉树（2D）/八叉树（3D）方法 21
3.5 阵面推进法和Delaunay三角划分结合算法 22
4 其他网格生成技术 23
4.1 自适应网格 23
4.2 混合网格 25
4.3 动网格 26
4.4 曲面网格 27
4.5 重叠网格 28
5 网格生成软件 29
5.3 Gambit 29
5.2 ICEM CFD 30
5.1 TrueGrid 32
5.2 Gridgen 34

33
1 概述
计算流体力学作为计算机科学、流体力学、偏微分方程数学理论、计算几何、数值分析等学科的交叉融合，它的发展除依赖于这些学科的发展外，更直接表现于对网格生成技术、数值计算方法发展的依赖。
在计算流体力学中，按照一定规律分布于流场中的离散点的集合叫网格（Grid），分布这些网格节点的过程叫网格生成（Grid Generation）。网格生成是连接几何模型和数值算法的纽带，几何模型只有被划分成一定标准的网格才能对其进行数值求解，所以网格生成对CFD至关重要，直接关系到CFD计算问题的成败。一般而言，网格划分越密，得到的结果就越精确，但耗时也越多。1974年Thompson等提出采用求解椭圆型方程方法生成贴体网格，在网格生成技术的发展中起到了先河作用。随后Steger等又提出采用求解双曲型方程方法生成贴体网格。但直到20世纪80年代中期，相比于计算格式和方法的飞跃发展，网格生成技术未能与之保持同步。从这个时期开始，各国计算流体和工业界都十分重视网格生成技术的研究。上个世纪90年代以来迅速发展的非结构网格和自适应笛卡尔网格等方法，使复杂外形的网格生成技术呈现出了更加繁荣发展的局面。现在网格生成技术已经发展成为CFD的一个重要分支，它也是计算流体动力学近20年来一个取得较大进展的领域。也正是网格生成技术的迅速发展，才实现了流场解的高质量，使工业界能够将CFD的研究成果——求解Euler/NS方程方法应用于型号设计中。
随着CFD在实际工程设计中的深入应用，所面临的几何外形和流场变得越来越复杂，网格生成作为整个计算分析过程中的首要部分，也变得越来越困难，它所需的人力时间已达到一个计算任务全部人力时间的60%左右。在网格生成这一“瓶颈”没有消除之前，快速地对新外形进行流体力学分析，和对新模型的实验结果进行比较分析还无法实现。尽管现在已有一些比较先进的网格生成软件，如ICEM CFD、Gridgen、Gambit等，但是对一个复杂的新外形要生成一套比较合适的网格，需要的时间还是比较长，而对于设计新外形的工程人员来说，一两天是他们可以接受的对新外形进行一次分析的最大周期。要将CFD从专业的研究团体中脱离出来，并且能让工程设计人员应用到实际的设计中去，就必须首先解决网格生成的自动化和即时性问题，R.Consner等人在他们的一篇文章中，详细地讨论了这些方面的问题，并提出：CFD研究人员的关键问题是“你能把整个设计周期缩短多少天？”。而缩短设计周期的主要途径就是缩短网格生成时间和流场计算时间。因此，生成复杂外形网格的自动化和及时性已成为应用空气动力学、计算流体力学最具挑战性的任务之一。
单元（Cell）是构成网格的基本元素。在结构网格中，常用的2D网格单元（图1.1）是四边形单元，3D网格单元（图1.2）是六面体单元。而在非结构网格中，常用的2D网格单元还有三角形单元，3D网格单元还有四面体单元和五面体单元，其中五面体单元还可分为棱锥型（楔形）和金字塔形单元等。

（a）三角形
（b）四边形
图1.1 常用的2D网格单元

（a）四面体
（b）六面体
（c）五面体（棱锥）
（d）五面体（金字塔）
图1.2 常用的3D网格单元

现有的网格生成技术一般可以分为：结构网格，非结构网格和自适应网格，此外还有一些特殊的网格生成方法，如动网格，重叠网格等。本文将重点介绍结构网格和非结构网格，因为这两种是CFD研究中应用最为广泛的两种网格生成技术。

图1.3 结构网格（左）和非结构网格（右）示意图

2 结构网格
结构网格是正交的，排列有序的规则网格，网格节点可以被标识，并且每个相邻的点都可以被计算而不是被寻找，例如i，j这个点可以通过i+1，j和i-1，j计算得到。采用结构网格方法的优势在于它很容易地实现区域的边界拟合；网格生成的速度快、质量好、数据结构简单；易于生成物面附近的边界层网格、有许多成熟的计算方法和比较好的湍流计算模型，因此它仍然是目前复杂外形飞行器气动力数值模拟的主要方法，计算技术最成熟。但是比较长的物面离散时间、单块网格边界条件的确定以及网格块之间各种相关信息的传递，又增加了快速计算分析的难度，而且对于不同的复杂外形，必须构造不同的网格拓扑结构，因而无法实现网格生成的“自动”，生成网格费时费力。比较突出的缺点是适用的范围比较窄，只适用于形状规则的图形。其发展方向是朝着减少工作量，实现网格的自动生成和自适应加密，具有良好的人机对话及可视化，具有与CAD良好的接口，并强调更有效的数据结构等。
结构网格主要分为常规网格、贴体坐标法（Body-Fitted Coordinates）和块结构化网格。常规网格是网格生成技术中最基本、也是最简单的，本章重点介绍后面两种方法。
2.1 贴体坐标法
在对物理问题进行理论分析时，最理想的坐标系是各坐标轴与所计算区域的边界一一符合的坐标体，称该坐标系是所计算域的贴体坐标系。比如直角坐标是矩形区域的贴体坐标系，极坐标是环扇形区域的贴体坐标系。贴体坐标又称适体坐标、附体坐标。
从数值计算的观点看，对生成的贴体坐标有以下几个要求：
① 物理平面上的节点应与计算平面上的节点一一对应，同一簇中的曲线不能相交，不同簇中的两条曲线仅能相交一次；
② 贴体坐标系中每一个节点应当是一系列曲线坐标轴的交点，而不是一群三角形元素的顶点或一个无序的点群，以便设计有效、经济的算法及程序。要做到这一点，只要在计算平面中采用矩形网格即可，所以贴体坐标系生成的是结构网格；
③ 物理平面求解区域内部的网格疏密程度要易于控制；
④ 在贴体坐标的边界上，网格线最好与边界线正交或接近正交，以便于边界条件的离散化。
生成贴体坐标的过程可以看成是一种变换，即把物理平面上的不规则区域变换成计算平面上的规则区域，主要方法有微分方程法，代数生成法，保角变换法三种。
2.1.1 微分方程法
微分方程法是20世纪70年代以来发展起来的一种方法，基本思想是定义计算域坐标与物理域坐标之间的一组偏微分方程，通过求解这组方程将计算域的网格转化到物理域。其优点是通用性好，能处理任意复杂的几何形状，且生成的网格光滑均匀，还可以调整网格疏密，对不规则边界有良好的适应性，在边界附近可以保持网格的正交性而在区域内部整个网格都比较光顺；缺点是计算工作量大。该方法是目前应用最广的一种结构化网格的生成方法，主要有椭圆型方程法、双曲型方程法和抛物型方程法。
2.1.1.1 椭圆型方程
以求解椭圆型偏微分方程组为基础的贴体网格生成思想最早是由Winslow于1967年提出的。1974年，Thompson、Thames及Martin系统而全面地完成了这方面的研究工作，为贴体坐标技术在CFD中广泛应用奠定了基础。此后，在流体力学与传热学的数值计算研究中就逐渐形成了一个分支领域——网格生成技术。文献中所谓的TTM方法就是指通过求解微分方程生成网格的方法在（TTM系上述三人姓的首字母）。用椭圆型方程生成的贴体网格质量很高，而且计算时间增加不多，不仅能处理二维、三维问题，而且还能处理定常和非定常问题，该方法成功实现了双流道泵叶轮内三维贴体网格的自动生成。
用椭圆型方程生成网格时的已知条件是：
① 计算平面上ξ，η方向的节点总数及节点位置。在计算平面上网格总是划分均匀的，一般取Δξ=Δη=1（0.1或其他方便的数值）；
② 物理平面计算区域边界上的节点设置，这种节点设置方式反映出我们对网格疏密布置的要求，例如估计在变量变化剧烈的地方网格要密一些，变化平缓的地方则应稀疏一些。
需要解决的问题是：找出计算平面上求解域的一点（ξ，η）与物理平面上一点（x，y）之间的对应关系，如图2.1所示。

（a）物理平面
（b）计算平面
图2.1椭圆型方程生成网格的问题表述图示

我们如果把（x，y）及（ξ，η）都看成是各自独立的变量，则上述问题的表述就是规定了一个边值问题，即已经知道了边界上变量（x，y）与变量（ξ，η）之间的对应关系（相当于第一类边界条件），而要求取在计算区域内部它们之间的关系。
① 从物理平面上来看，把ξ，η看成是物理平面上被求解的因变量，则就构成了物理平面上的一个边值问题：即已知道物理平面上与边界点（xB，yB）相应的（ξB，ηB），要求出与内部一点（x，y）对应的（ξ，η）。在数学上描述边值问题最简单的椭圆型方程就是Laplace方程。根据Laplace方程解的唯一性原理，可以把ξ，η看作物理平面上Laplace方程的解，即：
（2.1）
同时在物理平面的求解区域边界上规定ξ（x，y），η（x，y）的取值方法，于是就形成了物理平面上的第一类边界条件的Laplace问题。
② 从计算平面上来看，如果从计算平面上的边值问题出发考虑，则情况就大为改观，因为在计算平面上可以永远取成一个规则区域。所谓计算平面上的边值问题，就是指在计算平面的矩形边界上规定x（ξ，η），y（ξ，η）的取值方法，然后通过求解微分方程来确定计算区域内部各点的（x，y）值，即找出与计算平面求解区域内各点相应的物理平面上的坐标。实际上用椭圆型方程来生成网格时都是通过求解计算平面上的边值问题来进行的。为此需要把物理平面上的Laplace方程转换到计算平面上以ξ，η为自变量的方程。
利用链导法以及函数与反函数之间的关系，可以证明：在计算平面上与式（2.1）相应的微分方程为：
（2.2）
其中
（2.3）
从数值求解的角度，偏微分方程（2.2）的求解没有任何困难，它们是计算平面上两个带非常数源项的各向异性的扩散问题。由于参数α，β，γ把（x，y）耦合在一起，因而两个方程需要联立求解（采用迭代的方式）。在获得了与计算平面上各节点（ξ，η）相对应的（x，y）以后，就可以计算各个节点上的几何参数（xξ，xη，yξ，yη，α，β，γ）。
2.1.1.2 双曲型方程
如果所研究的问题在物理空间中的求解域是不封闭的（如翼型绕流问题），此时可以采用双曲型偏微分方程来生成网格。用双曲型偏微分方程来生成二维网格的方法是Steger和Chaussee于1980年提出的，随后，Steger和Zick将该方法推广到三维情况。这种生成方法通常是物面出发，逐层向远场推进，适用于没有固定远场边界网格的生成，在二维情况下，其控制方程为：
（2.4）
第一个方程控制网格线的正交，第二个方程控制网格单元尺度的分布，Ω为单元面积分布函数。在η=0（物面）上给定网格节点分布作为初值，然后沿η方向逐层推进生成网格。其优点是不用人为地定义外边界且可以根据需要直接调整网格层数；缺点是由于双曲型方程会传播奇异性，故当边界不光滑时，会导致生成的网格质量较差。所以，该方法通常用于生成对外边界的位置要求不严的外流计算网格或嵌套网格。
2.1.1.3 抛物型方程
采用抛物型方程来生成网格的思想是由Nakamura于1982年提出来的，这种方法生成网格的过程为：从生成网格的Laplace或Poisson方程出发，对方程中决定其椭圆特性的那一项作特殊处理，从给定节点布置的初始边界（设为η=0）出发，在ε=0及ε=1的两边界上按设定的边界条件（即节点布置），一步一步地向η=1的方向前进。其优点是概念简单，通过一次扫描就生成了网格而不必采用迭代计算；同时又不会出现双曲型方程的传播奇异性问题。
2.1.2 代数生成法
代数生成法实际上是一种插值方法。它主要是利用一些线性和非线性的、一维或多维的插值公式来生成网格。其优点是应用简单、直观、耗时少、计算量小，能比较直观地控制网格的形状和密度；缺点是对复杂的几何外形难以生成高质量的网格。
2.1.2.1 边界规范化方法
所谓边界规范化方法（Boundary Normalization）就是指通过一些简单的变换把物理平面计算区域中不规则部分的边界转换成计算平面上的规则边界的方法，这些变换关系式因具体问题而异。下面通过一些例子来说明。
① 二维不规则通道的变化：如图2.2所示一个二维渐扩通道的上半部，给定了不规则的上边界的函数形式为y=x2，1≤ x ≤2。则可采用下列变换把上边界规范化：
（2.5）
这里xt为上边界节点的x值。对于一条边界为不规则的二维通道，只要规定了不规则边界上y与x之间的关系式，都可以用这种方法来进行变换。

图2.2 不规则二维通道
图2.3 梯形区域的变换

② 梯形区域的变换：如图2.3所示的一个梯形区域可以通过以下公式变换成计算平面上边长可以调节的矩形：
（2.6）
其中F2（x）和F1（x）分别为梯形上下边的y与x的关系式，a与b为调节系数（放大或缩小），而且F2（x）和F1（x）不必为直线，曲线也行（但与垂直x轴的直线只能有一个交点）。
③ 偏心圆环区域的变换：图2.4a所示的偏心圆环区域可以采用变换转化成为计算平面上的一个矩形（图2.4b）：
（2.7）
偏心圆环中的自然对流就可以用这类变换生成网格。

图2.4 偏心圆环的变换
图2.5 可用双边界法生成贴体坐标的区域

2.1.2.2 双边界法
对于在物理平面上由四条曲线边界所构成的不规则区域，可以采用一种具有通用意义的方法来生成网格，这就是“双边界法”（Two-Boundary Method）。如图2.5所示，设在物理平面上有一不规则区域abcd，其中ab，cd为两不直接联接的边界。首先选定这两条边界上的η值，设分别为ηb和ηt，于是该两边界上的x，y仅随ξ而异。这些因变关系应该预先取定，设为：
（2.8）
下标b与t分别表示底边与顶边。
为简便起见，计算平面上的ξ，η取在0-1之间，这里暂取η1=0，η2=1，则以上式子可写成：
（2.9）
为了确定在区域abcd内各点的ξ，η值，一种最简单的方法是取为上、下边界函数关于η的线性组合，即：
（2.10）
其中f1（η）=1-η，f2（η）=η，这样生成的网格，在物理平面的边界上网格线与边界是不垂直的，为了生成与边界正交的网格，f1（η），f2（η）需要取为三次多项式，且在式（2.10）中要增加两条边界上xb，yb，xt及yt对ξ的导数项。
图2.6a所示的梯形如果用双边界法转换，可取：
（2.11）
则按式（2.9）得：
（2.12）
这就相当于把y方向的长度规范化。这一变换所得出的物理平面上的网格线显然不与x=0及x=1两条直线正交。对于物理平面的计算边界上的节点设置为均分情形（为5×5的节点布置），用双边界法得到的物理平面上的网格如图2.6b所示。

（a）
（b）
图2.6 双边界法例题

2.1.2.3 无限插值方法
双边界法还可以看成是构造了一种插值的方式，即把上、下边界上规定好的xt（ξ），yt（ξ）及xt（ξ），yt（ξ）通过插值而得出内部节点的（x，y）与（ξ，η）间的关系。如果同时在四条不规则的边界上各自规定了（x，y）与（ξ，η）的关系，这种关系式是可以解析的，也可以给出离散的对应关系。设分别为xb（ξ），yb（ξxt（ξ），yt（ξ），xl（η），yl（η）及xr（η），yr（η），其中下标l，r表示左右，如图2.7a所示，则可以采用下列变换（插值）得到物理平面上计算区域内任一点（x，y）与（ξ，η）的关系：
（2.13）
式（2.13）所规定的插值可以把四条边界上规定的对应关系连续地插值到区域内部，插值的点数是无限的，因而称为无限插值（Transfinite Interpolation，TFI）。可以这样理解：如果把ξ=（0，1），η=（0，1）分别代入到式（2.13），可以得出四条边界的（x，y）与（ξ，η）的关系式，因而在0<ξ<1，0<η<1的范围内式（2.13）给出了求解区域内节点的位置据四条边界给定的关系进行插值的方式。应用无限插值方法生成的网格如图2.7b。

图2.7 无限插值方法生成的网格

2.1.3 保角变换法
保角变换，又保角映射、共形映射，是复变函数论的一个分支，是从几何学的角度来研究复变函数，将二维不规则区域利用保角变换理论变换成矩形区域，并通过矩形区域上的直角坐标网格构造二维不规则区域贴体网格。和其他方法相比，在变换过程中的需要引入的额外项数目最少，变换的偏微分方程相对简单。随着复变函数论和微分几何学的发展，保角映射的理论和方法得到进一步发展，其中基于Schwarz-Christoffel的保角变换具有更大的灵活性，在二维的边界处理中应用广泛。这种方法的优点是能精确的保证网格的正交性，网格光滑性较好，在二维翼型计算中有广泛应用；缺点是对于比较复杂的边界形状，有时难以找到相应的映射关系式，且只能应用于二维网格。
2.2 块结构化网格
上节介绍的贴体网格求解不规则几何区域中的流动与换热问题的方法，可以用来求解一大批不规则区域中的流场和温度场，其中TTM方法的提出大大促进了有限容积法、有限差分法处理不规则区域问题的发展。但由于实际工程技术问题的复杂性，仍然有不少不规则区域中的问题难以用贴体坐标方法解决。本节介绍另外一种有效处理不规则计算区域的方法——块结构网格（Block-Structured Grid）。
2.2.1 基本思想
块结构化网格又称组合网格（Composite Grid），是求解不规则区域中的流动与传热问题的一种重要网格划分方法。从数值方法的角度，又称区域分解法（Domain Decomposition Method）。采用这种方法时，首先根据问题的条件把整个求解区域划分成几个子区域，每一子区域都用常规的结构化网格来离散，通常各区域中的离散方程都各自分别求解，块与块之间的耦合通过交界区域中信息的传递来实现。于是，采用这种方法的关键在于不同块的交界处求解变量的信息如何高效、准确的传递。
采用块结构化网格的优点是：（1）可以大大减轻网格生成的难度，因为在每一块中都可以方便地生成结构网格；（2）可以在不同的区域选取不同的网格密度，从而有效低照顾到不同计算区域需要不同空间尺度的情况，块与块之间不要求网格完全贯穿，便于网格加密；（3）便于采用并行算法来求解各块中的代数方程组。
2.2.2 两种基本形式
块结构化网格可分为拼片式网格（Patched Grid）与搭接式网格（Overlapping Grid），前者在块与块得交界处无重叠区域，通过一个界面相接（图2.8a）；后者则有部分区域重叠（图2.8b），这种网格又称杂交网格（Chimera Grid）。

（a）拼接式
（b）搭接式
图2.8 块结构化网格的两种类型

在块与块的交界处网格信息传递的常用方法有D-D型（D-Dirichlet，即第一类边界条件传递）D-N型（D-Neumann）两种。在D-N传递中一个块在交界处给出第一类边界条件而另一块则在交界处给出第二类边界条件。下面对拼片式与搭接式网格来说明D-N型及D-D型的传递方法。
2.2.3 D-N型信息传递方法
为说明方便，以如图2.9所示两块的公共边界AB上信息传递方法为例来说明，为了求解块1（密网格块）中的离散方程，需要有一个东侧邻点的值。为此将块1的ξ方向的网格线延伸一格，与疏网格区的CD相交于S点。S点的值可以根据CD线上相关位置的插值得到。一般可以去线性插值直到三阶插值，以获得所需的变量值。对于变化剧烈的变量，高阶插值反而会导致不合理的结果，宜采用线性插值。

（a）
（b）
图2.9 界面上信息的D-N传递

类似地将粗网格块2的ξ方向网格线延伸一格，交块1中的网格线EF与Q，P点。对于粗网格这条延伸边界采用由密网格的密度（如热流密度）式通量插值以获得相应的粗网格边界上的值。假设粗网格延伸边界上PQ上的热流密度为qc，则有：
（2.14）
其中Δη为η方向的网格补偿，下标c及f分别表示“粗”（corse）与“密”（fine），Nf为密网格中位于P-Q范围内的控制容积界面面积进入P-Q的百分数。为保证界面上的守恒性，对密网格再CD线上得到的值还应根据下列界面上的守恒进行调整：
（2.15）
如果上述条件不成立，就可以对这些插值得到的值做总体修正。图2.10是采用上述方法计算得到的分叉扩散器中的流动。

图2.10 分叉扩散器的块结构化网格

2.2.4 D-D型信息传递方法
为了计算图2.11这种情形的流动与换热，可以采用图2.11所示的这种组合网格。这里两个圆柱面附近区域采用极坐标，其余部分则采用直角坐标。这两种网格时独立地设置，并不考虑相互间要正好联接起来，但彼此间要有重叠的区域。设极坐标区I的外边界为曲线aa，II的外界为圆弧bb，而直角坐标网格去II的外边界为cc。

图2.11 搭接式网格

在图2.12中画出了重叠区内两种坐标系节点插值情形。在重叠区内，一种网格系统边界节点的值，可以利用与之相邻的另一网格系中的四个节点的值按现行插值原则得出。例如对图2.12中的P点，有：
（2.16）
而对中的P点，则有：
（2.17）

（a）
（b）
图2.12 搭接式网格重叠区内的插值

3 非结构网格
3.1 概述
同结构化网格的定义相对应，非结构化网格是指网格区域内的内部点不具有相同的毗邻单元。非结构化网格技术主要弥补了结构化网格不能解决任意形状和任意连通区域的网格剖分的缺陷。因此，非结构化网格中节点和单元的分布可控性好，能够较好地处理边界，适用于复杂结构模型网格的生成。非结构化网格生成方法在其生成过程中采用一定的准则进行优化判断，因而能生成高质量的网格，容易控制网格大小和节点密度，它采用的随机数据结构有利于进行网格自适应，提高计算精度。从定义上可以看出，结构化网格和非结构化网格有相互重叠的部分，即非结构化网格中可能会包含结构化网格的部分。
非结构化网格技术从上世纪60年代开始得到发展，到90年代时，非结构化网格的文献达到了它的高峰时期。由于非结构化网格的生成技术比较复杂，随着人们对求解区域的复杂性的不断提高，对非结构化网格生成技术的要求越来越高。从现在的文献调查的情况来看，非结构化网格生成技术中只有平面三角形的自动生成技术比较成熟（边界的恢复问题仍然是一个难题，现在正在广泛讨论），平面四边形网格的生成技术正在走向成熟。而空间任意曲面的三角形、四边形网格的生成，三维任意几何形状实体的四面体网格和六面体网格的生成技术还远远没有达到成熟。需要解决的问题还非常多。主要的困难是从二维到三维以后，待剖分网格的空间区非常复杂，除四面体单元以外，很难生成同一种类型的网格，需要各种网格形式之间的过渡，如金字塔形，五面体形等。
非结构化网格技术的分类，可以根据应用的领域分为应用于差分法的网格生成技术（常常称为Grid Generation Technology）和应用于有限元方法中的网格生成技术（常常称为Mesh Generation Technology），应用于差分计算领域的网格要除了要满足区域的几何形状要求以外，还要满足某些特殊的性质（如垂直正交，与流线平行正交等），因而从技术实现上来说就更困难一些。基于有限元方法的网格生成技术相对非常自由，对生成的网格只要满足一些形状上的要求就可以了。一般来说，非结构网格生成方法可以分为以下几类。
3.2 阵面推进法
阵面推进法（Advancing Front Method）的思想最早由A. George于1971年提出，目前经典的阵面推进技术是由Lo和Lohner等人提出的。阵面推进法的基本思想是首先将待离散区域的边界按需要的网格尺度分布划分成小阵元（二维是线段，三维是三角形面片），构成封闭的初始阵面，然后从某一阵元开始，在其面向流场的一侧插入新点或在现有阵面上找到一个合适点与该阵元连成三角形单元，就形成了新的阵元。将新阵元加入到阵面中，同时删除被掩盖了的旧阵元，以此类推，直到阵面中不存在阵元时推进过程结束。其优点是初始阵面即为物面，能够严格保证边界的完整性；计算截断误差小，网格易生成；引入新点后易于控制网格步长分布且在流场的大部分区域也能得到高质量的网格。缺点是每推进一步，仅生成一个单元，因此效率较低。
在生成初始阵面和新的三角形单元时。需要知道局部网格空间尺度参数，这可以由背景网格提供，对背景网格的要求是它能完全覆盖计算区域。早期的阵面推进法采用非结构化背景网格，背景网格的几何形状与拓扑结构及其空间尺度参数通过人为给定，这种方法的缺点是人工介入成分多，不易被使用者掌握，生成的非结构化网格光滑性难于保证，进行空间尺度参数插值运算时需要进行大量的搜索运算，降低了网格的生成效率。
赵斌等利用编制的计算程序对环形（图3.1）和NACA009翼型（图3.2和图3.3）通过设置点源控制内部网格疏密进行了网格剖分。

图3.1 环形区域网格划分

图3.2 NACA009翼型的非结构网格
图3.3 NACA009翼型的非结构网格的放大图

许厚谦和王兵提出了一种新的阵面推进算法——多点择优推进阵面法，它是在分析目前已有的几类生成方法的基础上改进而来。此算法的基本策略是推进阵面法，在映射平面的帮助下，实现对曲面的直接三角形网格划分，同时在划分结束后，能快速对格点进行Laplace松弛。与常规推进阵面法最大的差别在于：本方法在得到活动阵面的理想推进点时，不是由解析公式获得，也不是从曲面的几何信息中插值得来，而是从曲面上预先布好的点集中搜集到的一个最优点。此点集的规模很大，足够反映曲面的所有信息，最终的非结构网格格点数目仅占这个点集的千分之一左右。本方法的主要步骤:
① 生成背景网格。在背景网格中放置一定数目的源项，通过求解Poisson方程，可以实现网格尺度在整个流场域的自动分布。通过改变源项的数目、位置、尺度和强度，可以方便、有效地控制背景网格尺度的分布。和非结构背景网格相比，手工工作量小，且由于背景网格的尺度由求解方程获得，使得尺度分布更加光顺。
② 曲面边界离散。背景网格生成后，就可以对曲面边界进行剖分，得到初始阵面初始阵面如果质量不高，不仅影响网格质量，甚至会导致阵面推进失败。
③ 曲面离散成点集，存于二维数组。利用一系列纵横交错的线条来描述曲面，这些纵横线条被赋以整数参数来标识，根据曲面的定义，该点也对应一个空间点（x，y，z），同时网格交点也有了参数坐标（i，j），即将所有网格点的空间坐标存于3个二维数组中，即：x（i，j），yi，j），z（i，j），这种思想类似于将空间曲面表示成二元参数曲面。
④ 将边界离散点定位于参数平面（找出各点的参数i，j）。边界离散点已经事先确定，接下来需要找出各点在参数面上最接近的参数坐标（i，j）。这可以用低效率的逐一比较法，也可以用高效率的搜索法完成。如果（i，j）点对应的空间坐标与边界离散点的坐标不重合，则将边界离散点的坐标赋给x（i，j），y（i，j），z（i，j）。
⑤ 在空间曲面上进行阵面推进。首先选定一活动阵面，求出其空间长度S，两端点A，B的参数（iA，jA），（iB，jB）；其次从背景网格中求出该活动阵面中点的网格尺度L（作为理想等腰三角形之腰长），最后从活动阵面的中点M之参数（iM，jM）出发，搜索一个最优点P（iP，jP），其到两端点的空间距离最满足背景网格尺度要求，搜索算法的好坏直接关系到网格生成的成功与否及网格质量。因为空间曲面三角形与参数平面三角形的拓扑结构完全一样，因此可以在二维参数平面判断新三角形与阵面是否相交，及新三角形是否包含已有阵面点。
这种方法借助于参数平面的拓扑结构，对曲面直接进行三角形网格化，在划分结束后，还可以方便地对曲面网格进行Laplace格点松弛光顺。克服了传统曲面映射法造成的网格变形问题，且方法简单，在描述曲面时，仅仅要求能在曲面上布置一个网格点集，不需记录曲面的法向矢量、曲率等，适用面广（图3.4）。

图3.4 采用多点择优阵面推进法得到的运载火箭表面网格

3.3 Delaunay三角划分
Delaunay三角划分方法是在19世纪50年代Dirichlet提出Voronoi图的基础上发展而来的，是目前应用最广泛的网格生成方法之一。Delaunay三角形划分的步骤是：将平面上一组给定点中的若干个点连接成Delaunay三角形，即每个三角形的顶点都不包含在任何其他不包含该点三角形的外接圆内，然后在给定的这组点中取出任何一个未被连接的点，判断该点位于哪些Delaunay三角形的外接圆内，连接这些三角形的顶点组成新的Delaunay三角形，直到所有的点全部被连接。
要满足Delaunay三角剖分的定义，必须符合两个重要的准则：① 空圆特性：Delaunay三角网是唯一的（任意四点不能共圆），在Delaunay三角形网中任一三角形的外接圆范围内不会有其它点存在（图3.5）。② 最大化最小角特性：在散点集可能形成的三角剖分中，Delaunay三角剖分所形成的三角形的最小角最大。从这个意义上讲，Delaunay三角网是“最接近于规则化的“的三角网，具体的说是指在两个相邻的三角形构成凸四边形的对角线，在相互交换后，六个内角的最小角不再增大（图3.6）。

图3.5 空圆特性
图3.6 最大化最小角特性

Delaunay剖分是一种三角剖分的标准，实现它有多种算法，其中Lawson算法是最一种最经典的算法。逐点插入的Lawson算法是1977年提出的，该算法思路简单，易于编程实现。基本原理为：首先建立一个大的三角形或多边形，把所有数据点包围起来，向其中插入一点，该点与包含它的三角形三个顶点相连，形成三个新的三角形，然后逐个对它们进行空外接圆检测，同时用Lawson设计的局部优化过程LOP进行优化，即通过交换对角线的方法来保证所形成的三角网为Delaunay三角网。
Lawson算法的基本步骤是：
① 构造一个超级三角形，包含所有散点，放入三角形链表。
② 将集中的散点依次插入，在三角形链表中找出其外接圆包含插入点的三角形（称为该点的影响三角形），删除影响三角形的公共边，将插入点同影响三角形的全部顶点连接起来，从而完成一个点在Delaunay三角形链表中的插入。
③ 根据优化准则对局部新形成的三角形进行优化，将形成的三角形放入Delaunay三角形链表。
④ 循环执行上述第2步，直到所有散点插入完毕。
上述基于散点的构网算法理论严密、唯一性好，网格满足空圆特性，较为理想。由其逐点插入的构网过程可知，遇到非Delaunay边时，通过删除调整，可以构造形成新的Delaunay边。在完成构网后，增加新点时，无需对所有的点进行重新构网，只需对新点的影响三角形范围进行局部联网，且局部联网的方法简单易行。同样，点的删除、移动也可快速动态地进行。但在实际应用当中，这种构网算法当点集较大时构网速度也较慢，如果点集范围是非凸区域或者存在内环，则会产生非法三角形。

图3.7 离散点集合
图3.8 Delaunay三角剖分产生的网格

Delaunay三角划分的优点是具有良好的数学支持；生成效率高；不易引起网格空间穿透；数据结构相对简单；而且速度快，网格的尺寸比较容易控制。缺点是为了要保证边界的一致性和物面的完整性需要在物面处进行布点控制，以避免物面穿透。
3.4 四叉树（2D）/八叉树（3D）方法
Yerry和Shephard于1983年首先将四/八叉树法的空间分解法引入到网格划分领域，形成了著名的四叉树/八叉树方法。其后许多学者对该方法进行了完善和发展，提出了修正的四叉树/八叉树方法。
修正的四叉树/八叉树方法生成非结构网格的基本做法是：先用一个较粗的矩形（二维—）/立方体（三维）网格覆盖包含物体的整个计算域，然后按照网格尺度的要求不断细分矩形（立方体），使符合预先设置的疏密要求的矩形/立方体覆盖整个流场，最后在将矩形/立方体切割成三角形/四面体单元。图3.9是用一种基于四叉树方法生成的三角形网格示意图。

图3.9 一种基于四叉树方法形成的三角形网格

该方法的优点是：① 基本算法很简单而且树形的数据结构对于很多拓扑的和几何的操作（比如寻找邻近节点等）都很有效；② 可以很容易地把自适应网格细分合成进来（通过细化一个叶节点或删除一个子树）；③ 可以与固体模型很好地结合，因为在许多几何操作上它们用到了同样的基本思想；④ 网格生成速度快且易于自适应，还可以方便地同实体造型技术相结合。缺点是由于其基本思想是“逼近边界”且复杂边界的逼近效果不甚理想，所以生成网格质量较差。
3.5 阵面推进法和Delaunay三角划分结合算法
阵面推进法生成的网格具有质量好，边界完整性好的特点；而Delaunay三角划分法生成网格具有高效率和良好的数学支持的特点。在1993年Rebay首次提出这两种网格生成方法结合之前，它们一直是作为“竞争对手”来进行描述的，此方法一提出，就立刻吸引了众多学者对该方法进行研究，并提出了许多改进的方法。
阵面推进法的实施过程为：从边界网格出发，内部的点通过阵面推进法来生成，然后利用Delaunay算法对这些点进行逐点插入，重复以上过程直到网格的尺寸达到要求尺寸。阵面推进法的优点是网格的质量好、边界逼近效果好、网格生成效率高和有良好的数学支持；缺点是对于边界网格的依赖性较大，边界网格的质量直接影响网格划分的结果。

4 其他网格生成技术
4.1 自适应网格
为了使网格划分与数值求解过程结合，在90年代初提出了自适应网格技术（adaptive grid）。自适应网格方法是指计算中，在某些变化较为剧烈的区域，如大变形、激波面、接触间断面和滑移面等，网格在迭代过程不断调节，将网格细化，做到网格点分布与物理解的耦合，从而提高解的精度和分辨率的一种技术。自适应网格希望在物理解变动较大的区域网格自动密集，而在物理解变化平缓区域网格相对稀疏，这样在保持计算高效率的同时得到高精度的解。自适应网格技术能对预先划分好的网格按照用户所定义的误差准则，自动进行误差判断并进行网格疏密程度的调整，即以最适合于所求解问题的方式来布置节点及确定其间的联系，进而取得令人满意的数值模拟计算结果。自适应网格原则上只需定义一种描述问题几何特性的初始网格及可接受的误差水平，计算机可自动产生能够实现这一有效水平的网格，可以提高分析效率和计算结果的可靠性，除了在航空和气象领域应