1、2 三角剖分算法介绍2.1Delaunay 三角剖分理论基础2.1.1VORONOI 图1908 年,俄国数学家 M.G. Voronoi 提出了 Voronoi 图 1,它是计算几何学上非常重要的工具之一,被广泛应用于各个领域中。Voronoi 图是由平面区域中连接两邻点的线段的中垂线所形成的区域。它又叫 Dirichlet 图或 Thiessen 多边形 2。Voronoi 图是自然界中的宏观物体和微观物体以空间距离相互作用所形成的一种网格结构,应用的范围相当广泛,特别是在计算几何学领域的应用 3。Voronoi 图是一种关于平面或空间区域划分的基础数据结构。100 多年来,它被广泛应用在
2、与几何信息相关的各个领域。随着计算机科学技术的不断普及和发展,Voronoi 图的应用领域也在不断地扩大。对于 Voronoi 图的应用,以 90 年代以来的应用更为突出。Voronoi 图的本质属性是由空间实体几何唯一确定的,不是通过其他方法强加上去的。从 Voronoi 图的数字几何角度来看,它是针对平面 n 个离散点而言的,把平面分为若干区域,每一个区域包含一个点,该点所在的区域就是到该点距离最近的点的集合。设 p1,p 2 是平面上两点,L 是线段 p1p2 的中垂线, L 将平面分成两部分 LL和 LR,位于 LL 内的点 pl 具有特性:d(p l,p1)d(pl,p2),其中 d
3、(pl,pi)表示 pl 与 pi 之间的欧几里得距离,i=1, 2。这意味着,位于 LL 内的点比平面上其他点更接近点 p1,换句话说,L L 内的点是比平面上其他店更接近于 p1 的点的轨迹,记为V(p1),如图 2.1 所示。如果用 H(p1,p2)表示半平面 LL,而 LR=H(p2,p1),则有V(p1)=H(p1,p2),V(p 2)=H(p2,p1)。给定平面上 n 个点的点集 S,S=p 1,p2, ,pn。定义:V(pi)= (,ijjHp即 V(pi)表示比其他点更接近 pi 的点的轨迹是 n-1 个半平面的交,它是一个不多于 n-1 条边的凸多边形区域,称为关于 pi 的
4、 Voronoi 多边形或关于 pi 的Voronoi 域 。对于 S 中的每个点都可以作一个 Voronoi 多边形,这样的 n 个Voronoi 多边形组成的图称为 Voronoi 图。Voronoi 多边形的每条边是 S 中某两点的连线的垂直平分线,所有这样的两点连线构成一个图,称为 Voronoi 图的直线对偶图。如图 2.2 所示,虚线表示 voronoi 图,实线表示其对偶图。图 2.1 平面中两离散点 图 2.2 Voronoi 图Voronoi 图有以下一些性质:性质 1.n 个点的点集 S 的 Voronoi 图至多有 2n-5 个顶点和 3n-6 条边。性质 2.每个 Vo
5、ronoi 点点好是三条 Voronoi 边的交点。性质 3.Voronoi 图的对偶图实际上时点集的一种那个三角剖分,该三角剖分就是 Delaunay 三角剖分 4。性质 3 实际上告诉我们可以用构造 Voronoi 图的方法来求平面点集的三角剖分,但实际上这种方法较少使用,因为算法的效率不好。2.1.2Delaunay 三角剖分的基本概念1、基本概念(1)域分割给定平面上的 n 个不相重的散乱数据点,对每个散乱数据点构造一个域,使该域内的任一点离此散乱点比离其他散乱点更近,这种域分割就是上节介绍的 Voronoi 图,也称 Dirichlet 域分割,由上节可知,域边界其实就是连接两个相邻
6、散乱点的直线的垂直平分线。(2)Delaunay 三角化对平面上的散乱数据点进行域分割后,将具有公共域边界的散乱点对相连形成的三角剖分称为 Delaunay 三角剖分。(3)优化对三角网格进行优化,就是要使三角网格整体上尽量均匀,避免出现狭长三角形,也就是获得 Delaunay 三角化。2、三角剖分优化准则在三角剖分过程中,人们往往先用一种简单的方法构造散乱点的初始三角剖分,然后对其进行优化以获得 Delaunay 剖分。优化的方法取决于所采用的优化准则,平面三角剖分最常用的优化准则有 Thiessen 区域准则、最小内角最大准则、圆准则。Sibson 证明了这三个准则的等价性,并指出符合这三
7、个准则的三角剖分只有一个,即 Delaunay 三角剖分。(1) 最小内角最大准则对一个严格凸的四边形进行三角化时,有两种选择,最小内角最大准则就是要保证对角线两侧两个三角形中的最小内角为最大,如图 2.3 所示。图 2.3 对角线优化(2) 圆准则严格凸四边形中的三个顶点确定一个圆,如果第四个顶点落在圆内,则将第四个顶底与其相对的顶点相连,否则将另外两个顶点相连,这个准则称为圆准则。也就是说,符合圆准则的三角剖分中,任一三角形的外接圆内不应该包含其他点,如图 2.4 所示。图 2.4 空外接圆准则(3) 局部优化局部优化是指对任意一个凸四边形的对角线,依据某种优化准则做交还测试后所得到的三角
8、剖分。(4) 全局优化当三角形剖分 T 中每一条内边上的两个三角形所形成的凸四边形都满足局部优化标准时,称该三角剖分 T 满足全局优化。(5) 退化前面已经讲过 Delaunay 三角剖分满足圆准则,即任一三角形的外接圆内不能包含其他的点,如果规定三角剖分中任一三角形的外接圆内和外接圆上都不能有其余顶点,则称为标准的 Delaunay 三角剖分,在实际应用中,往往存在四点共圆的情况。我们把四点或四点以上的散乱数据点共圆的情况称为退化,相应的三角剖分称为准 Delaunay 三角剖分。3、二维任意区域内点集的 Delaunay 三角划分概念及基本定理将所有共边的 Voronoi 多边形的中心连接
9、起来所形成的三角网称为Delaunay 三角剖分 5。在区域 D 内的点集 V 的三角划分(记为 T(D,V))是指将 V 中的点用互不相交的直线段连接起来,使得域 D 内的每一个区域都是三角形,同时这些连线和三角形都在区域 D 内,包含在区域 D 内的边称为内边。二维任意区域点集的 Delaunay 三角划分 5(简记为 DATA)是指所有的内边都满足局部优化的 T(D,V)。4、基本定理下面的定理是闵卫东在文献 6中提出来的,其证明参见文献 6。定理 2.1:设 B,P 点在 AC 边的同侧,如果角 APC 小于角 ABC,则三角形 ABC 的外接圆内不包含点 P,参见图 2.5。图 2.
10、5 三角形外接圆 图 2.6 局部优化边 AB定理 2.2:设 AB 边为三角形 ABC 和三角形 ABE 所共有,如果四边形ACBE 的内角 CBE180或 CAE180,则 AB 边是局部优化的,见图 2.6。定理 2.3:设三角形 ABC 是 DATA(D,V)的一个三角形, P 是 V 中的一点且不是 ABC 的三个顶点,如果线段 PC 在域 D 内且它与 AB 边有交点,则ABC 的外接圆内部不包含点 P。定理 2.4:(DATA 的 Circle 准则)三角剖分 T(D ,V)是DATA(D,V)的充要条件是,对于 T(D,V)的任何一个三角形 ABC,V中的任何点 P,如果它与
11、A,B,C 相连的线段 PA,PB,PC 都在域 D 内,则P 不包含在三角形 ABC 的外接圆内部。定理 2.5:任何一个 T(D,V)都可以通过有限次的局部优化操作将其转化为 DATA(D,V) 。定理 2.6:DATA(D,V)在所有的 T(D,V )中平均形态比最大。以上几个定理实际上作为 Delaunay 三角划分呢算法的理论依据在许多以Delaunay 原理为基础的网格剖分算法中经常用到。5、三角剖分时需要满足的几个约束条件在进行网格剖分过程中,生成三角形时,必须满足一下三个约束条件:1. 三角形相互之间是不相交的,即两个三角形除端点外不应该有别的交点,此条件称为三角形相交约束。2
12、. 三角形相互之间是互不包含的,即任意一个三角形不能完全包含其他的三角形,此条件称为三角形包含约束。3. 三角形全部落在区域之内,此条件称为三角形合法约束。2.2Delaunay 三角剖分算法介绍Delaunay 三角剖分算法是最常用的一种平面点集的剖分算法,实际上,Delaunay 三角剖分算法是基于 Delaunay 原理的一类三角剖分算法的总称。总体来说,Delaunay 三角剖分算法可分为以下几种。2.2.1 翻边算法1977 年,Lawson 基于最小内角最大优化准则提出了一种二维点集 Delaunay三角剖分翻边算法 7。该算法的思想是:给定一个二维离散点集,首先对其进行一次初始三
13、角剖分,然后根据最小内角最大化准则判断初始三角剖分中形成凸四边形的共边三角形是否满足优化准则;如果不满足,就交换四边形的两条对角线。一次循环直到所偶的三角形都满足最小角最大优化准则为止。通过计算可以得到,翻边算法的时间复杂度为 O(n2)。根据前面的理论分析知道,运用翻边算法扑粉出来的所有三角形边都是 Delaunay 三角剖分的边,那么最终得到的三角网格一定是 Delaunay 三角网格。目前,二维的翻边算法发展比较完善,利用翻边法对二维离散点集进行Delaunay 三角剖分原理比较简单、得到的网格质量也较好。1989 年,Barry Joe基于空外接球准则提出了一种三维 Delaunay
14、三角剖分的翻面法 8。翻面法的原理跟翻边法的原理是一样的,只不过在三维空间利用翻面法会比较复杂,有待进一步研究。2.2.2 逐点插入算法随着 Delaunay 三角剖分的不断发展,大量剖分算法不断被提出来。1981 年,一种逐点插入算法被 Bowyer 和 Watson 提出来 910。该算法的基本思想是:首先构造点集的一个初始三角剖分,然后逐一地在当前三角剖分中插入一个新点;在插入过程中,要根据 Delaunay 三角剖分空外接圆住着呢进行三角网格优化,直到整个点集为空集为止。该算法首先构造一个极大三角形,使得要剖分的点集中所有点都落在该三角形里面;其次,逐一地插入各个新点,每插入一个点,必
15、须搜索位于外接圆内部的点的三角形;然后,将这些三角形从三角形队列中删除,可以形成多边形空腔;最后,将插入点和多边形空腔连接起来,构成若干个以插入点为共同顶点的新三角形。定位新插入点的位置是逐点插入算法的一个关键步骤,直接关系到剖分网格的好坏和剖分速度的快慢。由于包含插入点的三角形的外接圆位于插入点的附近区域,只要能够准确定位新插入点在当前三角剖分中的位置,便可快速检索到外接圆包围它的三角形。基于定位插入点位置的重要性,现在大量改进的逐点插入算法主要都体现在定位新插入点的位置上。由于逐点插入算法比较简单、易懂,可以将其推广到三维或者更高维空间的三角剖分。唯一不足的就是该算法时间复杂度较高,很多专
16、家学者在研究算法的时候也会将该因素重点考虑进去,从而设计出简单、高效的三角剖分算法。2.2.3 分割合并算法Hoey 和 Shamos 提出了一种用于生成 Voronoi 图的分割合并算法,由于Voronoi 图与 Delaunay 图互为对偶图,进而可以间接得到 Delaunay 三角网格 11。后来 Dwyer、Katajaine、L.J.Guibas 等人提出了用于直接生成 Delaunay 三角网格的分割合并算法。分割合并算法主要根据数学递归思想,递归将点集分割成规模相当的两部分子集,然后对分割出来的点集进行 Delaunay 三角剖分,再递归地将两个相邻Delaunay 三角剖分进行
17、合并组合,从而生成了整个点集的 Delaunay 三角网格。改进后的分割合并算法时间复杂度为 O(nlogn),目前最好的分割合并算法是Dewall 算法。该算法不仅可以应用于二维空间,还可对三维或更高维空间进行三角剖分,虽然分割合并算法实现的时间效率很高,但它的编程却非常复杂。以上三种算法都是比较典型的 Delaunay 三角剖分算法,应用范围比较广泛。除了这几种经典算法外,其实还有很多种其他的 Delaunay 三角剖分算法。在选择三角剖分算法的时候,必须从难易程度、效率高低和贴近实际与否等方面选择,其中逐点插入算法是目前最为简单、最为流行的一种 Delaunay 三角剖分算法,该算法思想
18、简单、易懂,可推广到维数更高空间的三角剖分。2.2.4Bowyer-Watson 算法Bowyer-Watson 算法是一种先形成粗原始网格,再逐步插入新点进行细化的算法,是一种逐点插入算法。假定已有初始剖分为 T,向已有网格内插入新点的 Bowyer-Watson 算法如下:Step1:插入一个新点 P 到现有的 Delaunay 三角剖分中。Step2:寻找并删除所有外接圆包含 P 点的三角元素,形成一个 Delaunay 空腔。Step3:连接 P 点和 Delaunay 腔壁上所有各点,形成新的 Delaunay 三角剖分。应用 Bowyer-Watson 算法进行网格剖分的整个过程简
19、述如下:Step1:家里一个覆盖整个计算区域的粗原始网格。Step2:采用 Bowyer-Watson 算法将边界点逐点插入原始网格。边界点由人工事先给定,并假设边界点的分布式合理的。Step3:删除计算区域以外的三角形,并确保边界面的正确三角剖分(进行拓扑相容性处理) 。Step4:用 Bowyer-Watson 算法逐步向计算区域内插入新点(对于给定点集,按一定顺序插入即可,若内点是自动生成的,则需要按一定的策略生成内点,直到所有内点都被插入或者网格达到一定的扑粉要求为止。Step5:对所生成的网格进行拓扑相容性检查、网格光顺等工作。3 矿井涌水量地质模型的实现3.1 初始粗网格的建立矿井
20、涌水量地质模型的建立,首先是要将煤层顶底板的曲面离散成三角形网格。本文中采用逐点插入法作为主要的算法。逐点插入法的第一步是要建立一个包含所有计算域顶点的初始粗网格。常用的建立初始网格的方法是包容盒法,一般可建立两种形式的包容盒,一种是建立矩形的包容盒 12,另一种是三角形包容盒 13。无论是给定点集或者是给定区域的三角剖分均可采用该方法。如果只给出边界区域,则内点需要自动产生。本文所建立的三角剖分是给定区域边界,自动生成内点。这种方法在使用过程中需要对生成的三角形进行边界拓扑相容性检测。因为拓扑相容性检测的步骤,本文给出一种无需建立三角形粗网格或矩形粗网格的方法,直接将给定区域的多边形划分成三
21、角形初始粗网格。该方法首先将给定边界按顺时针或者逆时针的方向将边界顶点进行排列,然后计算每一个顶点的角度大小,作为判断顶点角度的指标。随后每次选择最小的顶点,将其与相邻的顶点连成一个三角形,重新计算该顶点两侧的顶点角度大小,然后再选择最小角进行处理。由于多边形并非是完全的凸多边形,边界上很可能会出现凹顶点。在这种情况下,需要判断顶点的凹凸性。本文使用矢量面积法 14来判断多边形顶点的凹凸性。如果多边形边界以逆时针方向排列,那么为了判断某一顶点的凹凸性,可以用与该顶点相邻的两个顶点,共三个顶点来判断该顶点的凹凸性。计算这三个顶点形成的三角形的矢量面积,如果结果为正,那么该顶点为凸顶点,如果为负,
22、那么该顶点为凹顶点。计算公式可简单描述如下:设三个顶点为:A、B、C,那么:0.5()()()()SxAyCxByA判断完顶点的凹凸性之后,需要计算顶点的角度(弧度) 。可以使用相邻三个顶点形成的两个向量来计算顶点的角度。 、 的数量积可以计算出两向C量的夹角,再结合顶点凹凸性的判断,加上或者减去 180便是该顶点的角度。Delaunay 剖分的初始化步骤可用以下算法描述:算法 3.1Step1:判断顶点凹凸性,计算每一个顶点的角度。Step2:从顶点数组中找到角度最小的顶点,连接该顶点两侧的顶点,形成一个三角形。并计入三角形数组。Step3:重新计算该顶点两侧新顶点的角度,并从顶点数组中删去
23、该顶点。Step4:从新的顶点数组中找到新的最小角度顶点,重复 step2 至 step3。直到顶点数组中的顶点个数等于 3 个。Step5:连接最后三个顶点形成一个三角形,计入三角形数组。按照以上算法,最后会形成一个初始的粗三角形网格,或者称为多边形的一个粗三角划分。如图 3.1,3.2。图 3.1 任意多边形区域图 3.2 初始网格划分3.2 简单多边形的 Delaunay 剖分在建立了多边形的初始粗网格以后,便可以开始向计算区域内插入点进行Delaunay 三角剖分。本文采用自动布点技术向计算区域内插入点。采用的方法是:选择面积最大三角形的最长边的中点最为插入点,建立三角网格。这样需要设
24、置一个面积的控制量,以便控制剖分的程度。即三角形的个数或者三角网格的疏密程度。这是一个可以改变的量,可以方便地控制三角网格的大小。当插入一个点以后,需要搜索整个三角形数组,进行空外接圆检测,如果新插入点落在该三角形的外接圆内(包含边界) ,那么该三角形被添加到需要优化的三角形数组中,所有三角形检测完后,符合要求的三角形会形成一个多边形空腔,删除这些三角形的公共边,连接插入点与多边形空腔的各个顶点,形成新的多个三角形,加入三角形数组,删除被合并的三角形。这样一直到三角形面积符合控制要求结束。便完成了简单多边形的 Delaunay 剖分。算法可描述如下:算法 3.2Step1:设定三角形面积控制单
25、元,计算三角形数组中所有三角形的面积大小。Step2:选择面积最大的三角形的最长边作为新的插入点,并将该点加到顶点数组的后面。Step3:搜索整个三角形数组,找出所有插入点落在三角形外接圆中的三角形,计入优化三角形数组,并将三角形的边计入边数组,并且不能有重复边,若有重复边,便要从数组中删除该边。Step4:连接插入点与边数组中每条边的两个顶点,形成新的三角形,加入三角形数组。形成新的三角形数组。Step5:删除优化三角形数组中的三角形,更新整个三角形数组,完成一次点的插入。Step6:重复 step2、step3、step4,直到三角形的面积符合控制要求,结束点的插入,这样就完成了简单多边形
26、的 Delaunay 剖分。剖分实例如图 3.3。图 3.3 Delaunay 剖分在剖分的过程中需要注意几个问题。如果插入点落在多边形边界上,便会形成同一直线上的两条不同的边被同时计入边数组,当连接三角形的时候,便会出现三个顶点落在一条直线上的三角形,所以需要对这种三角形进行特别的处理。不然会出现错误。处理方法很简单,在连接三角形的时候,需要计算插入点与两个顶点连接形成的直线的斜率,比较这两条直线的斜率,如果斜率相等,那么表示这三个点处在同一条直线上,这样便不能连接三角形。算法表示如下:算法 3.3Step1:计算边数组中将要连接的那条边的两个顶点与插入点形成的直线的斜率。Step2:判断两
27、条直线的斜率。Step3:如果斜率相等,那么那条边的连接性置为 0,如果斜率不等那么直线的连接性置为 1。Step4:连接边的两个顶点与插入点,形成一个三角形,计入三角形数组。在计算斜率的时候,由于数据类型的限制,三角形的斜率不一定严格相等,所以当三点一线的时候斜率会有微小的误差,这一误差需要用一常量来作为判断斜率是否相等的限制,可以将这一误差设定为 1,用来帮助判断斜率。3.3 用 VB 实现的 Delaunay 剖分算法本文使用 Visual Basic 6.0 作为开发环境来实现 Delaunay 剖分算法。分为四个模块:1、 公共数据模块2、 初始化模块3、 剖分模块4、 数据输出模块
28、3.3.1 公共数据模块公共数据模块定义有程序中使用的数据类型,包括:1、 顶点类,存储有顶点的 x、y、z 坐标数据,以及顶点是否是边界格点的判断变量。2、 三角形类,存储有每个三角形的三个顶点编号,编号既是三角形顶点在顶点数组中的序号。3、 临时顶点类,该数据类型主要用于初始化时存储临时顶点数据,除了包括 x、y、z 坐标外,还有弧度值与该顶点在原顶点数组中的序号。4、 边类,该数据类主要用于剖分时存储被优化三角形的边信息,数组中的边将成为新生成的三角形的边。模块中还定义有全局顶点数组和全局三角形数组,全局顶点个数和全局三角形个数。全局顶点数组是一个只增不减的数组,顶点数会随着剖分的进行而
29、不断增加,且数据不会改变。全局三角形数组是一个动态数组,数组元素的个数和元素的取值都会随着程序的运行而改变。初始化时已事先为两个数组定义了元素个数。可随剖分的细分程度修改。3.3.2 初始化模块初始化模块用于对多边形区域进行初始的网格化。形成初始的三角形存入三角形数组。该模块中先定义临时顶点数组,该临时顶点数组继承全局顶点数组的全部数据。并存储每一个顶点在全局顶点数组中的序号。随后计算每一个顶点的角度,存入临时顶点数组,使用算法 3.1 判断顶点的凹凸性并计算顶点的角度。经测试矢量面积法能判断顶点的凹凸性,由此可以计算每一个顶点的角度大小,在程序中,角度值转换为弧度值计算。该模块中包括判断顶点
30、凹凸性的函数,初始化网格的过程和绘制网格的过程。这三个函数或过程在初始化过程中被调用。初始化网格的过程中,先找到顶点角度最小的角,将该顶点两侧的顶点连接形成一个三角形,插入全局三角形数组。并在临时顶点数组中删除该顶点。重新计算该顶点两侧的顶点角度大小,并更新临时顶点数组。在计算的过程中,需要考虑序号为 1、2、n-1、n 这四个顶点。因为数组数据类型的限制,当最小角顶点序号是这四个值时,顶点两侧的顶点序号并不连续,所以需要分开计算。用判断语句分别判断这四种情况,分别计算顶点两侧的顶点角度。其余情况可用循环语句计算。执行完初始化网格过程,全局三角形数组包含最初的三角形网格,全局顶点数组不变。调用
31、绘制网格过程绘制三角形。3.3.3 剖分模块该模块用于将多边形剖分成需要的 Delaunay 三角形网格。模块首先建立边数组,用于存储需要连接的边信息。首先计算现有三角形的面积,找到面积最大的三角形。然后取该三角形最长边的中点作为新的插入点,将插入点计入全局顶点数组。搜索整个全局三角形数组,找到插入点落在其外接圆中的三角形,并记录该三角形在全局三角形数组中的序号,将其三条边送入边数组等待优化。在将边送入边数组的过程中,需要判断边是否已经存在边数组中,如果已经存在则要在边数组中删除该边,如果不存在则存入边数组。搜索完全部三角形后,形成了一个记录需要优化三角形的序号数组和边数组,边数组中的边形成一
32、个连续的多边形空腔,插入点可能位于空腔内部,或者空腔的某一条边上。如果位于边上,则需要检测插入点落在哪条边上,而这条边则不能作为新三角形的边。使用插入点与边两个顶点的连线的斜率来判断插入点是否落在边上,这样便可除去三点一线的三角形。按顺序连接插入点与边数组中的每条边,形成新的数个三角形,将这些三角形送入全局三角形数组。如果符合要求的三角形只有一个,则要专门计算插入三角形的个数,形成两个新的三角形,删除一个旧的三角形。随后检测需要删除的三角形,从数组尾端开始向前逐个删除需要删除的三角形。这样便完成了一个点在多边形区域中的插入,更新了全局顶点数组和全局三角形数组。重新计算数组中每个三角形的面积,找
33、到面积最大的三角形,如果面积值大于某一设定值,则继续向其中插入点,直到满足要求为止。这样便完成了多边形区域的 Delaunay 剖分。3.3.4 数据输出模块该模块主要用于顶点数据和三角形数据的输出,以便在随后的计算中使用。将需要的数据输出到一个文本文件中,留待后用。程序界面如图 3.4 所示:图 3.4 程序运行界面3.4 对剖分算法的评价该算法可以将简单的多边形剖分成 Delaunay 三角网格。并且算法省去了建立初始大三角形或者矩形的步骤,直接将多边形划分成初始粗网格。然后进行剖分。这样就省去了进行边界拓扑相容性检查的步骤。和其他算法相比,该算法在时间复杂度上并无优势,剖分本质上是逐点插
34、入法,时间复杂度较高。该算法的一个缺点是,当三角形长边位于边界上时,不能保证该三角形是锐角三角形,也是就该三角形有可能是钝角三角形。目前尚无有效的办法解决此问题,这也是选择布点策略时产生的结果。3.5 网格节点编号优化3.5.1 节点编号优化算法有限差分分析时经常会碰到具有稀疏矩阵刚度的方程式。减少稀疏矩阵的带宽可以减少有限差分法的计算时间,和减少占用计算机的存储空间。在有限差分法的计算过程中,可通过节点重新编号来减少带宽。对于节点重编号的研究较多,最典型的有:Cuthill 和 Mckee15基于图论的节点分层编号方法,可在一定程度上减少刚度矩阵的带宽,这个方法后来被称为 CM 方法。Gib
35、bs,Poole 和 Stockmeyer16提出了更加有效的基于图论的减少稀疏矩阵带宽和外形的方法,此方法后来被称为 GPS 方法。Akhras 和 Dhatt17通过计算节点商方法来缩小带宽,该方法编程方便,此方法后来被称为 AD 方法。欧阳兴 18提出了前沿法与矩形法。CM 方法和 GPS 方法,采用了复杂的图论方法 19,没有估算算法本身的复杂度。实现比较困难,如果考虑节点编号算法本身时间消耗大,就会降低节点编号算法的效率。AD 法虽然实现比较方便,但是在网格重划和网格加密时节点编号优化速度慢。前沿法和矩形法对于板料单元和平面二维单元特别适用。黄志超,包忠诩,周天瑞 20等从分析最佳最
36、简单的矩形出发,总结出一种新的节点重编号方法,是一种改进的 AD 法。如图 3.5 是最典型的最佳二维网格节点编号,由此来分析其特性,总结出一般规律。图 3.5 二维网格 15 节点优化该 15 节点序号是最佳的编号法,它的最优同单元节点编号的最大差值为4。由该图的编号法可总结出以下规则:(1) 节点商是从小到大升序排列的;(2) 最小最大节点编号和也是从小到大升序排列的;(3) 最大节点基本上是按升序排列的,不过有编号重复的情况。由此可总结出一种新的节点重编号方法。算法 3.4:Step1: 计算出每个节点的相关单元节点编号总和;Step2: 计算出该节点所相关单元的最大节点编号与最小节点编
37、号;Step3: 计算出最小节点编号与最大节点编号之和;Step4:计算出该节点的节点商,将该节点的相关单元节点编号总和除以该节点的相关单元个数所得的值为节点商;Step5: 把节点商作为该节点重编号的依据,节点商小的编号也小,节点商大的编号也大,如果节点商相等则把最小最大节点编号和作为重编号的依据;step6: 重复(1)到(5)直到带宽不能再减小为止。经计算该算法能减少 50%左右的带宽。图 3.6 为一 495 个节点的网格形成的稀疏矩阵每一行的带宽分布图。图 3.6 节点带宽分布图3.5.2 节点编号优化算法的实现本文使用 VB 语言实现节点编号算法。输入数据为剖分完成后的节点数组和三
38、角形数组。该算法对应程序的计算水头模块。在公共数据模块中定义一新的数据类 OptVertex。该数据类包含节点在顶点数组中的序号 number,节点商 NodeTrade 和最大最小节点编号之和 maxmin 三个数据。在计算水头模块定义一该类型的数组作为节点优化数组,命名为OptVer。数组长度和全局顶点数组相同。算法首先将全局顶点数组的序号存入节点优化数组,然后开始循环。计算节点优化数组的每一个节点,搜索全局三角形数组如果三角形中包含该节点的编号,便将该节点相关单元数加 1,并根据三角形的三个顶点在全局顶点数组中的序号从节点优化数组中找到这三个顶点,并将其在节点优化数组中对应的序号累加存入
39、一变量中。同时要记录相关单元中节点序号的最大值与最小值。三角形数组搜索完成后,便得到了该节点的相关单元数、相关单元节点编号总和、最大节点编号和最小节点编号。接着计算该节点的节点商和最大最小节点编号之和,分别存入该节点的对应变量中。这样依次计算所有节点优化数组中所有节点的相关数据存入数组中。对计算完成的节点进行排序,使用选择排序法,将第一个位置的节点和后面所有的节点做比较,节点商大的序号也大,节点商小的序号也小。如果节点商相等,比较节点的最大最小节点编号之和,小的序号也小,大的序号也大。然后比较第二个位置的节点和后面的所有节点,依次排序下去。这样就完成了一次节点编号的优化。接着计算每一个节点的带
40、宽,初始带宽为节点的个数。然后计算每一个节点在稀疏矩阵中的带宽,计算方法为,从三角形数组中找到包含该节点的三角形,计算每一个节点和该节点之差的绝对值,也就是在稀疏矩阵的一行中包含数据的元素到该行对角线的距离的绝对值。计算出最大的带宽,然后和初始带宽比较,如果比初始带宽小,说明稀疏矩阵带宽减小了,应当进行下一次节点编号优化。将当前带宽作为下一次比较的带宽,继续循环重复以上步骤对节点编号进行优化排序。直到当前带宽不再改变为止。这样就完成了节点编号的优化。3.6 三角形网格有限差分法3.6.1 三角形均衡网格差分方程的建立三角形网格差分方程是计算网格中格点水头的重要步骤。也是有限差分方法的核心。差分
41、方程形成的系数矩阵是主要的数据结构。考虑以格点 i 为中心的多边形网格 Di 作为均衡区(如图 3.7) ,下面可根据达西定律和水均衡原理建立格点 i 的差分方程 21。图 3.7 多边形网格均衡区首先求出单位时间内通过 的周边流入 内的水量。iDi由图 3.7 可见,我们可以先求出通过三角形 ijk 内的边 和 流入 内的pbqiD水量。由达西定律有: 11nnjikiii ihhQTpbTbq(3.1)式中: 、 为流段 ij 和 ik 的平均导水系数; 、 、 、 为线段长度;ijik ijpbikq为从第 e 号三角形通过两线段 和 流入 内的水量。i pbqiD根据(式 3.1)的计
42、算方法,对格点 i 周围的所有三角形作类似计算,并求和得: 11 nnjikii ihhQTpbTbq侧(3.2)另外,由于本文所研究的计算区域为定水头边界,承压水,不考虑垂向补给量。所以在单位时间内,流入均衡区 内的总水量为:Q 侧 。由水均衡原理得iD均衡方程:(i=1,2,N)1s=niihQAt侧(3.3)式中: 为 的面积, 为给水度。iiDs这就是多边形网格的水均衡方程。由以上公式可以看出,在建立方程时,必须对格点 i 周围的每个三角形计算、 值,而在通常情况下,只知道格点的坐标,因此,为了方便,下面将pbijqk、 以及 均用格点坐标标出。ijiA4ijijbcpkij(3.4)
43、 4ikibcqji(3.5)并且记 ijkjikijby(3.6)ikjjikjicx(3.7)其中, =ij1()2ijijbc(3.8)接下来计算多边形 的面积 ,可以将 的计算分解成若干小三角形,分iDiAi别计算每个小三角形的面积。计算三角形 和三角形 的面积。pbibq2()16ijijkbcipb(3.9)同理可得 2()16ikijjbciqb(3.10)类似此法,在格点 i 周围的各个三角形中重复上述工作,即可求得 的值。iA由上述推导过程可见,要建立格点 i 的差分方程,只须对格点 i 周围的三角形逐一做类似的计算即可。3.6.2 三角形网格差分方程的解法由(式 3.2)和
44、(式 3.3)可知,在一个差分方程中,将涉及到格点 i 及其周围几个格点的水头值。在本文的算例中,由于一类边界格点是定水头边界,所以水头值是已知的。其余内格点是未知的,所以必须联立方程组求解。按前述方法,对每个内格点建立差分方程,若内格点有 N1 个,则共有 N1个方程,而这些方程包含的未知数也正好是 N1 个,所以方程数和未知水头数一样多。此外,如果把格点 i 的差分方程中的项 的系数排列在系数矩阵的主1nih对角线上,其系数矩阵为对角占优,故方程组存在唯一解,可用迭代法求解。SOR(Successive Over-Relaxation)法,即超松弛迭代法 22,是目前解大型线型方程组的一种
45、最常用的迭代加速方法 23。它是在高斯-塞德尔迭代法 24的基础上进行加速的,将前一步的结果 与高斯-赛德尔迭代法的迭代值 适当的kix 1kix加权平均,期望获得更好的近似值 。其迭代公式 2526如下:1i=1,2, ,n;k=0,1,2, 11() ()()inkkkkii jijixxba(3.11)SOR 法中 27的取值对迭代公式的收敛速度影响很大 28,它的好坏直接影响到加速的快慢。为了保证迭代过程的收敛,必须要求 02,超松弛迭代法取 12 29。算法 3.4Step1: 将内格点原始编号存入一个新的数组,获得系数矩阵大小。Step2: 搜索三角形数组,找到包含当前内格点的三角
46、形,计算当前内格点的多边形区域面积,周围格点的系数,存入数组中,直到所有内格点都处理完。Step3: 将数组中的系数数据按照内格点顺序存入系数矩阵中,如果格点为边界格点,将其计算结果放入右端已知项矩阵中。Step4: 使用 SOR 法解线性方程组,得到水头值,存入数组中。3.6.3 计算水头模块三角形网格差分方程的解法在程序的“计算水头”模块中实现。该模块前半部分为网格节点编号的优化。后半部分为三角形网格有限差分法的计算。该部分主要的程序代码为系数矩阵的形成。由于研究区域为定水头边界,水头值已知,且不随时间的变化而变化。所以这部分格点应作为已知量,放到线性方程组右端已知项中。因此,用于计算的未
47、知水头值个数即为内格点的个数。在已经优化的节点编号中需要除去边界格点。仅保留内格点。重新计算过的节点数组大小将只是内格点的个数。以此数组的大小来定义系数矩阵,右端已知项矩阵,初始水头值数组的大小。计算中包含的重要数据结构包括:Matrix(points,points) ,即系数矩阵,BMatrix(points) ,即右端已知项矩阵,x(points) ,即初始水头值向量,其中points 即为内格点的个数。以下描述程序的具体实现过程。(1) 首先形成新的节点数组,该数组不包含边界格点,仅包含内格点。这样,也即确定了系数矩阵的大小,重新定义系数矩阵的大小,方程右端已知项的大小。(2) 搜索全局
48、三角形数组,找到包含当前节点的三角形,逐个计算当前节点周围节点的系数,存入当前节点的系数数组中,计算当前节点的多边形区域面积,保存在变量中,待三角形全部搜索完毕,即获得当前节点多边形区域的面积和周围各个节点的系数。(3) 将数组中的各个节点的系数按照节点编号顺序存入系数矩阵对应位置,如果该节点为边界节点,便获得其水头值,将水头值与系数相乘加到右端已知项中。随后将已知项存入已知项矩阵对应位置。(4) 循环计算内格点的系数和右端已知项,存入系数矩阵和右端已知项矩阵对应位置,这样便完成了系数矩阵和右端已知项矩阵的形成。(5) 取得各个格点的初始水头值,存入初始向量数组中,将其和系数矩阵以及右端已知项
49、矩阵作为 SOR 法的输入数据。(6) 调用 SOR 法解方程组,迭代完成后即得到该时段末各个内格点的水头值。将其存入相应数组中。以上,便是按照水均衡法原理 30设计的矿井涌水量程序,程序获得水头值以后便可以交给相关部门预测涌水量,以及进行矿井水害防治等工作。6 10-108 工作面涌水量计算本文以霍州矿区白龙煤矿 480 水平,10-108 工作面作为研究区域。将区域做了简化,具体研究区域见图 6.1。并且假定工作面已经形成,然后开始涌水,工作面不随时间变化。图 6.1 研究区域边界及工作面位置外部多边形区域为本文研究区域的边界,边界为定水头边界,即水头不随时间变化,为一定值 550。含水层为 K2 灰岩,厚度为 10m,给水度 1/10000,渗透系数 K 为 2m/天。内部矩形区域为含水层下部巷道的位置,含水层中的水即会从矩形区域向下涌入巷道,该区域内部也为一定水头边界,为一定值480。渗透系数 K 为 5m/天,下渗路径长度为 5m,矩形区域面积为 10