2.3 河道贴体网格生成技术
河道贴体网格指的是能较好拟合复杂河道边界的网格，目前最具代表性的贴体网格有正交曲线网格，渐变三角形无结构网格，三、四边形混合网格等。
2.3.1 正交曲线网格
正交曲线网格的主要优点是即使对于比较复杂的河道，也能实现计算网格的贴体，且网格为正交网格，易于进行通量的计算。缺点有：描述物理现象的控制方程需进行变换，变换后的方程较复杂；对于河道急剧弯曲的情形，如鹅头型分汊河道，贴体与正交很难同时保证；进行网格局部加密时，需将加密区所在的行和列同时加密等。
1.网格变换方法
正交曲线网格可通过拉普拉斯变换、泊松变换或区域转换方程来实现。
（1）拉普拉斯变换。拉普拉斯（Laplace）变换的目的是将物理平面（x，y）上的不规则区域变换为计算平面（ξ，η）上的规则区域，并通过求解拉普拉斯方程，在物理平面和计算平面上生成一一对应的网格。
其控制方程为
控制方程可以变换为
在确定计算平面与物理平面对应区域的区域边界线和方程的边界条件后，通过求解式（2.22）和式（2.23），即可得到计算平面网格接点与物理平面网格节点的对应关系。
（2）泊松变换。同一几何形状或流动现象可以采用不同的坐标系来表达，坐标系之间可以有一一对应的关系相互转换。通常最基本的坐标系是笛卡尔直角坐标系（x，y），但对于某些复杂几何形状或流动现象采用其他坐标系（ξ，η）描述更为方便，因此需要进行坐标变换。对于曲线坐标系（ξ，η），J.F.Thompson发展了一种方法，他通过求解两个椭圆型偏微分方程（泊松方程）的解来寻找原坐标系统x，y与新坐标系统ξ，η之间的一一对应关系：
其中，P、Q（P、Q为0时为拉普拉斯变换）为一对“合并”或“聚散”函数（收缩因子），它主要用来调节网格线的位置变化，对变换结果有显著影响，其作用有两个：一是适当选择P、Q可使曲线网格的疏密根据需要分布；二是可以使变换后的曲线网格保持正交。
通过以（ξ，η）为自变量、以（x，y）为因变量的泊松方程可变换为
其中J=xξyη-xηyξ，如果使曲线网格保持正交，β取为0。
（3）区域转换方程。其主要原理是根据有势流的等势线和流线的正交机理生成正交曲线网格，它是在泊松变换基础上进行简化改进的一种方法。采用该方法生成的网格既能保持正交，又能控制网格疏密，而且计算简便。其核心是将式（2.24）、式（2.25）改写成如下形式：
其中
为了便于计算又不失一般性，取Δξ=Δη=1，则令
式中：gξ、gη分别表示x、y平面上曲线网格的长和宽，因此α与曲线网格的长宽比成比例。
将式（2.31）、式（2.32）改写成如下形式：
上二式右边即为式（2.24）、式（2.25）中的收缩因子P、Q，将P、Q代入式（2.26）、式（2.27）中即可得到平面（x，y）上正交曲线网格点坐标转换方程：
上二式中：Cξ、Cη为物理正交曲线坐标系中的拉梅系数。
在确定了计算平面与物理平面对应的区域边界线和方程的边界条件后，采用有限差分方法求解式（2.35）、式（2.36），就可自动生成物理平面（x，y）上正交曲线网格节点坐标。
（4）三种变换方法的比较。拉普拉斯变换能保证物理域网格线的贴体性、光滑性和边界正交性，但其致命的弱点是网格疏密无法控制，无法满足水力学计算所必需的合理的疏密分布。泊松变换通过增加和调节源项可实现网格的疏密控制，并保证网格的正交，但是其源项的确定和控制往往较为复杂，尤其对于复杂的边界条件情况。区域转换方程既能使生成的网格保持正交，又能控制网格疏密，而且计算和处理方便，因此对于河道模拟而言是一种难得的网格变换方法。
拉普拉斯变换方程是泊松变换方程的一个特例，而区域转换方程是基于泊松变换思想得到的一种简便、有效的变换方法。
2.边界正交处理方法
在对区域变换方程进行迭代求解过程中，为保证网格在边界的正交性，需对边界进行必要的处理。许多学者在此方面进行了研究，如固定一边，让另外一边变动，这样处理往往使部分边界正交性很好，而部分不好；又如通过计算域内靠近边界的内点向边界做垂线来确定边界点，这样能很好保证边界正交，但难以控制边界的疏密，有时导致网格线在边界上交叉或出现很小的网格步长。
根据多年实践经验的总结，我们提出了一种既简便又有效的边界处理方法。其主要思想为：首先对计算域边界线进行适当的光滑处理，得到基本能反映计算域外形的边界线；然后将光滑后的边界线看做“固定的轨道”，将要生成的网格边界点看做“带柄的滑轮”，让“滑轮”在“轨道”上有序滑动，且使“滑轮的柄”始终垂直于“轨道”。
实现方法如下：对光滑的边界线进行离散，令边界线不动，而网格的边界点可在边界线上以离散点为基础滑动。为保证边界网格的正交性，不需按传统的方法从靠近边界的内点向边界线做垂线来确定边界点，而是通过寻找内点与边界上离散点间的最短距离来确定边界点。这样处理的优点在于，无论如何迭代总能保证网格边界点落在不变的边界线上，且用最简单的方法来实现边界网格的正交。如图2.8所示，假设边界线j=N，为确定边界点（i，N）位置，只需从（i-1，N）到（i+1，N）间的边界离散点中（图2.8中圆形点）寻找与内点（i-1，N-1）距离最近的点即可，如图中P3点即为边界点（i，N）位置。
图2.8 寻找与边界垂直的点
对于边界的光滑，也采用了最为简单的方法，对于边界上的某点，将该点的位置向其前后两点连线的平分点移动即可，当然该方法适用于边界点较密的情形，否则会带来较大误差。该方法的最大优势在于能有效消除“锯齿”形不光滑边界。
3.应用实例
通过上面的分析说明，基于区域转换方程的贴体网格生成方法是适合于天然河道网格划分的一种简便、有效的方法。下面以长江下游镇扬河段的世业洲汊道段为例，将网格生成步骤说明如下：
（1）在CAD绘图软件中，绘制天然河道轮廓线，在不考虑水流漫堤情况下，一般以河道大堤堤线作为网格的初始边界线。
（2）沿河道从上游向下游布置控制断面（如图2.9中CS1～CS25），控制断面要能基本反映河道外形，断面可以为两点，也可以为多点（如CS7～CS19）。
图2.9 正交曲线网格生成图
（3）按指定的沿河宽方向节点间距，对各控制断面进行节点划分，连接两相临控制断面上相应的节点，得到沿河道方向的线段，按指定的沿河道方向断面间距进行断面划分，便可得到计算域的初始网格。
（4）采用有限差分方法迭代求解式（2.35）、式（2.36），便得到计算域的正交曲线网格。
2.3.2 任意平面域三角形网格
在进行三角形单元划分时，一般应遵循以下原则：①均匀原则：为了保证数值计算有较好的精度，要求三角形剖分是“正规”的，即要求每个单元中尽量不出现大的钝角，而且每个单元的外接圆和内切圆直径之比越接近越好；②不交叉原则：注意要使一个单元的顶点只能是相邻单元的顶点，而不是相邻单元边界上的非顶点，即任意两个单元的边与边不能交叉；③疏密原则：在计算域内，估计未知数变化剧烈的部分，可将网格分得密些；未知数变化缓和的部分，可将网格分得稀些。
非结构化网格生成方法主要有两种：Delaunay三角形化方法和阵面推进法。这两种方法均能实现网格生成的自动化，能够实现局部加密，能够加以推广生成三维非结构化网格。
Delaunay三角形化方法在20世纪八九十年代得到了迅速的发展，Delaunay三角化的依据来自G.F.Voronoi和Dirichlet分别提出的一种利用已知点集将平面划分成凸多边形的理论。这种理论的思想是：给定区域Ω及点集{Pi}，Pi位于Ω内部或者边界上，则存在每一个点的区域Di，区域Di内任意一点与点Pi的距离都比与{Pi}中其他点的距离近。这种划分方法将平面划分成了一系列不重叠的凸多边形，叫做Voronoi区域，并且使得
，且这种分解是唯一的。我们采用符号d（），表示欧几里得距离，则Di满足的关系可表示为
连接所有相邻凸多边形中的Pi点形成的三角形网格就是俄国数学家Delaunay提出的Delaunay三角网。而Delaunay三角网和Voronoi图互为对偶图（见图2.10）。
阵面推进法将边界离散信息视为扰动波，逐层向计算区域内部推进，直至扰动传递到计算区域内的每一点。这种方法网格生成过程直观，与Advancing Layer方法结合能在固壁附近生成数层高伸展比的三角形单元或四边形单元以更好地分辨边界层流动，能够实现方向性加密。但阵面推进法的缺点是需要人工给定背景非结构化网格上的网格空间尺度参数，网格的光滑性与疏密程度不容易控制，网格生成速度慢；可能出现畸形二角形网格；若对流动边界离散方式不恰当可能使网格生成过程崩溃，或者在流场计算后出现沿流动边界非物理压力振荡问题。
图2.10 Delaunay与Voronoi对偶图
在总结前人的基础上，基于阵面推进法和网格节点间距函数，我们提出了适合天然河道网格划分的方法。主要步骤为：首先确定计算域边界线，包括内边界及工程边界，然后对边界线进行离散，在沿河道进口向下游扩展直到计算河段出口，扩展过程中三角形边长由网格尺度函数确定，最后对生成的网格进行适当优化。
1.计算域边界的划分
首先绘制河道的轮廓线，一般由堤线和洲边界线构成。河道两岸堤线和河道进出口断面构成一逆时针的封闭线，而洲线则构成顺时针封闭线，还可选择工程轮廓线（不必封闭）或局部点作为网格需通过的点或线。
边界线上节点的疏密可通过下述方法控制：记节点的密度为k，当k小于、等于或大于1时，沿着一条直线的间隔将会逐渐变短、保持不变或逐渐变长，如图2.11所示。
图2.11 边界线的疏密控制
该直线上某点j的坐标可由下式计算：
为了使网格在生成过程中能对边界进行自动判断，有必要建立节点信息表（程序中用数组来表示），以记录所有节点的边界属性。如进口边界点用“2”表示，出口边界点用“3”表示，不过水洲滩、河岸及人工建筑物上边界节点用“1”表示，内部节点则用“0”表示（即非边界节点）。
2.网格疏密控制参数
当边界离散后，为控制区域内部网格的疏密过渡，引入区域内节点间距函数ρ（x，y），它可以通过利用边界离散点节点间距函数的加权平均计算得到，可通过定义边界点的节点间距函数来控制整个计算域的疏密，从而实现局部加密控制和疏密网格的均匀过渡。
式中：N为边界节点总数；ρ（xi，yi）为边界节点的间距函数；di为区域内任一点（x，y）距边界点（xi，yi）的距离；m为权函数，可取为2～3。
图2.12 判断某点是否在求解域内图示
3.判断一个点的位置是在域内还是域外
无论是在初始化三角形时把位于求解域外的三角形去掉，还是在自动加点时判断该点是否在某一三角形的外接圆内，都需要判断一个点的位置是在域内还是域外。参照相关文献，对复杂多连域发展出了一种判断方法，下面以图2.12为例，对其步骤加以说明。
（1）根据输入的边界网格点的数据，把每相邻两点连接成直线，用直线段来逼近边界的求解域，图2.12中内外边界的编号是连续的。设有一点P，坐标为xP，yP，现判断其是否在此求解域内。
（2）过x=xP做一垂直线，如果此直线不与该求解域任何直线段相交，则P点必在求解域外；如果此直线与求解域的一些直线段相交，则相交的直线段数目必为偶数。对图2.12中情形，这些直线段为19-20、29-30、34-35、87-88、79-80及66-67，相交线段数为6。
（3）计算各交点的y坐标。
（4）根据交点y坐标，将相交的直线段重新排序，从y值低的排列到y值高的。图2.12中应为66-67、79-80、87-88、34-35、29-30、19-20，这些线段以符号li（i=1，…，6）表示。
（5）判断P（xP，yP）是否位于任意两个直线段之间，而且其中位置低的那个直线段li值为奇数。如果两个判断都得出肯定的回答，则P点必在求解域内，否则P点不在求解域内。
4.节点与网格的同时生成
在引入节点间距函数后，可采用内部节点和三角形同时生成的前沿生成法生成Delaunay三角形。初始前沿为已离散的区域边界，记为非Delaunay边，优先考虑最长边，并尽可能保证新生成的边长度不断减小，以实现从大到小的过渡，具体步骤如下：
（1）从前沿队列中取出要扩展的有向边AB，记长度为a，如图2.13所示。
图2.13 三角形网格的推进
（2）以AB为底边，做AB的垂直平分线OD，OD的长度为h。
（3）根据上述的边界加权公式计算D的间距函数为ρ，并根据下列条件沿AB的垂直平分线调整OD的距离h：
（4）以D点为圆心，以h/3为半径做圆，对前沿队列中的所有节点进行判断，看是否有点落入圆中。若有，选择距D点最近的点C与AB边构成新三角形ABC，同时得到一条扩展边AC；若没有，则以D点为顶点构成新的三角形ABD，同时生成两条扩展边BD和AD。
（5）由新的扩展边开始重复上述步骤，直到前沿队列为空，整个计算域便布满三角形。
5.网格的优化
如果适当调整网格节点的位置，网格质量将进一步提高，这就是网格的优化。常用的一种网格优化方法称为Laplacatian优化方法。这种优化技术是通过将节点向这个节点周围的三角形所构成的多边形形心移动来实现的。如果Pi（xi，yi）为一个内部节点，Ni为与Pi相连的节点数，则优化技术可表示如下：
6.应用实例
图2.14为采用上述方法得到的天然河道三角形网格图，由中可见，所建立的三角形网格能很好地适应河道不规则的边界条件，可自由对重点关心部位进行局部加密，并保持疏密网格的良好过渡。
图2.14 三角形网格划分图
2.3.3 三四边形混合网格
对于主槽较窄、滩地较宽的河道，采用混合网格是最好的选择，将主河道划分为四边形网格，滩地划分为三角形网格，可以节省计算单元。
图2.15 单一河槽网格生成
1.四边形网格生成
（1）单一河槽网格生成。图2.15（a）为一概化的S形弯道河段。边界线由一系列控制点组成，其中上下游四个顶点组成河道进、出口的两个断面，通过给定沿水流方向的断面间距和垂直水流方向上的节点间距，就可生成贴体四边形网格，如图2.15（b）所示。
若要对局部河段网格进行加密，只需把局部河段按照上述方法调整断面和节点间距即可。
（2）多河槽网格生成。对于河道中不同的河槽，可按照以上方法生成不同的四边形网格，然后通过网格拼接技术，得到整个计算域的网格。
考虑到整体网格的一致性，各个不同汊道的有结构四边形网格，均采用无结构编码进行节点和单元的编号。
2.混合网格生成
对于主槽较窄、滩地较宽的河道，采用混合网格是最好的选择，一方面能保证主河槽在枯水条件下不会因为网格较粗而断流，另一方面可将滩地网格加粗，节省计算单元。可采用分区对接的方法生成混合网格，即首先在主河槽生成贴体四边形网格，然后在河道固定边界和四边形网格边界基础上，生成非结构三角形网格，然后进行拼接。在拼接过程中，在三角形网格和四边形网格的交界面上，不能出现非共享的边界点。
3.应用实例
图2.16给出了四边形网格与三角形网格的拼接效果，从中可见，四边形网格能保持较大的长宽比，即沿水流方向网格尺度远大于垂直水流方向，使得河道主槽中有较多的网格节点，又不占用过多的网格单元；而三角形网格能自由对任意区域进行网格划分，并方便地实现疏密网格的渐变。
图2.16 分块对接混合网格图