有限元分析中的二维Delaunay三角网格剖分

有限元分析中的二维Delaunay三角网格剖分摘要本文从有限元分析出发，引出三角网格剖分的概念。

随后着重介绍了二维平面点集的Delaunay三角剖分。

给出了一些重要的Delaunay三角形的定理和性质，也体现出了Delaunay三角剖分的优点。

接着重点分析了构造二维Delaunay三角形的空洞算法，并用程序完成了它。

最后又分析了算法中的不足，并给出论文改进的方法。

关键词：Delaunay三角形，V oronoi图，网格剖分III1 第一章绪论1.1网格剖分的背景有限元分析是数学的一个分支。

其思想是将复杂的问题简单化，然后进行处理。

处理办法是将整个研究对象分成一些有限的单元，然后对每个小单元做相应的处理，最后整合起来去逼近原来的整个对象。

所以我们可以看到，有限元分析中将单元剖分的越小，得到的近似值就会越逼近真实值。

但是往往我们需要处理的对象很复杂，需要的计算量也很大，人工很难完成。

在早起年代，这个问题也阻止了有限元分析的发展。

近年来，随着计算机的发展，带动了一些需要大量计算的科学领域的发展。

有限元分析就是其中一种，因为当计算机取代人力之后，其快速的计算能力作用愈发凸显，人们只需要控制相应的算法即可。

作为最常用的处理手段，被大大的发展了之后，有限元分析也被应用于诸多方面。

早期的有限元分析主要应用与航空航天和地质、地球物理方面，现在越来越多的在工程分析计算和计算流体力学中看见。

图 1.1图 1.2常见的有限元分析可以分为六大步骤：问题及求解域的定义、求解域的网格剖分、确定状态变量及控制方法、单元推导、总装求解和结果解释。

上述步骤又可被分为三大阶段：前置处理、计算求解和后置处理。

而在前置处理中网格剖分作为最重要又最复杂的一个步骤，其处理结果制约着有限元的最后逼近结果。

网格剖分有很多形式：二维的主要剖分形状有三角形、四边形，三维的有四面体、六面体。

在有限元分析中网格剖分有如下要求：1、节点合法性。

指每个单元的节点最多只能是其他单元的节点或边界点，而不能是内点。

2、单元相容性。

指每个单元必须在求解域的内部。

3、良好的单元形状。

指每个单元尽量最好是正的，比如二维是正多边形，三维是正多面体。

4、自适应性。

是指在剖分域中曲率大或其他参数变化较大的地方剖分越密，单元越小，越平滑或其他参数变化不大的地方单元可以稍微稀疏。

这样，既可以提高计算收敛速度，又可以提高逼近精度。

网格剖分的对象越是复杂，剖分的要求越高。

对于复杂三维实体，现在还没有成熟的算法。

1.2 网格剖分的发展网格剖分是几何模型和数值计算之间的桥梁。

1974年网格剖分首次被Thompson等人用椭圆方程方法构造出来。

之后Steger等人又提出了用双曲型方程来构造出网格。

在20世纪90年代，非结构网格和自适应笛卡尔网格等技术相继被提出，大大推动了CFD的发展。

直到如今，网格技术已产生好多方法，非结构网格、结构网格以及自适应网格等，其中最为常用的是非结构网格技术和结构网格技术。

结构网格是指正交的排列规则的网格。

它的特点是相邻的节点不需要遍历寻找就可以被计算出来。

生成结构网格的方法主要有：贴体坐标法和块结构化网格。

非结构网格和结构网格对应，是指内部节点没有毗邻单元。

其主要的方法有阵面、Delaunay三角剖分、四叉树和八叉树法。

本文主要针对于二维平面上的Delaunay三角剖分。

2第二章问题分析我们假设有这样有这样一个场景，当你拿着新买的碟子的边在盛食物，刚加了两个饼突然碟子沿中间裂开了。

我们知道这是碟子的质量问题，捡起来重买一个就可以。

但是在现实生活中我们遇到的好多工程，如桥梁、高楼等，他们对质量要求很高，在建筑之前必须对整个建筑做相应的应力分析，否则后果不堪设想。

我们以这个盘子为例，若有工程师愿意对此盘子做有限元分析，我们首先来看看它的单元怎么生成。

首先我们看第一步：问题及求解域。

很明显这个问题中我们是要研究整个盘子的应力分析，而其求解域就是整个盘子，我们可以简化一下，它就是一个数学圆，或者里面又套了一个小圆，其为底部的边界。

第二步：网格剖分。

这一步是关键，也是难点。

我们遇到的首个问题是，怎么把实物，此处为碟子，描述到计算机中，变成可用于计算的对象？我们有两种方法。

第一种：我们可以现实中丈量实物的长宽等尺寸；第二种：对实物进行拍照，然后对其图像处理，进行像素提取、边界提取、纹理处理等操作。

显然第一种简单，第二种很麻烦。

但是工程中，我们所面对的都是庞大有复杂的对象，或是测量成本很大，或是根本不能进行现实测量的，如山体卫星图等。

所以第二种方法虽然复杂，但是可以推广且尤为重要，故本文重点用第二种方法来完成此实物的网格剖分。

图 2.1所以我们首先要进行边界提取。

边界提取目前有成熟的算法，对于本文中的实物，我们采用有向边缘法，其算法理论在算法设计环节中重点介绍。

3其次我们需要对边界进行多边形逼近，所幸的是，我们得用计算机到的边界其实已经离散化了，所以这一步我们可以省略。

最后我们的重点到了，就是对上述的边界，或是点的集合进行三角网格剖分。

我们采用已经成熟的空洞算法，又名增量算法，来生成有限点集的Delaunay三角形。

下面是Delaunay三角形剖分的理论基础。

45 第三章 理论介绍3.1 Delaunay 三角形定义及性质在数学中我们有对偶图的概念。

即图二以图一的面为顶点，图一中共边的两面对应的图二中顶点为相邻的两点，则图二为图一的对偶图。

图 3.1而Delaunay 三角形是Voronoi 图的对偶图。

所以为定义Delaunay 图，我们首先定义Voronoi 图，这又需要用到点的领域的概念。

定义1、点的领域。

点的领域是指在平面的有限区域内，给定一些有限的点集}p ,p ,p ,{p n 321 ，对每个点n i p i ,2,1,=，其领域是指这样的点的集合},2,1),)(,(min ))(,(|),{(n j P y x d P y x d y x V j i i ===。

即每个点的领域中的点到其他给定点的距离大于或等于到该点的距离。

定义2、Voronoi 图。

在平面上给定有限的点集，每个点和相邻点之间线段的垂直平分线彼此相截围成的图是Voronoi 图。

所以我们看到，Voronoi 图实际上就是点集的领域图。

定义3、Delaunay 三角形。

Voronoi 图的对偶图就是Delaunay 三角形图。

在上述中我们会发现有些问题存在。

当所给的点集中存在四点共圆，且圆中无其他给定点，则Voronoi 图中会有四线共点，对应的Delaunay 三角形中会出现四边形。

如右上图。

但是不影响，此时我们链接四边形的任意一条对角线分成两个三角形即可。

下文中我们会看到这么处理不会影响Delaunay 三角形的性质，符合其空外接圆准则。

图 3.2图 3.3性质1、空外接圆准则。

所给的点集生成Delaunay三角形中，每个三角形都不会以其他给定点为内点，最多以其为边界点。

1定理1、给定点集所生成的Delaunay三角剖分中，每个单元都符合上述准则。

定理2、所给定点集的某一三角剖分，若每个单元都符合Delaunay三角准则，则这个三角剖分一定是Delaunay三角剖分。

3.2 Delaunay三角形的优点性质2、最小角最大准则。

即在给定点集的所有三角剖分中，我们会有Delaunay三角形中的最小角比其他剖分中最小角大。

这个准则和Delaunay三角形的空外接圆准则是一致的，但这条准则也让我们看到Delaunay三角剖分均匀的这一优点。

最小角最大准则，可以有效的降低1杨钦.限定Delaunay三角网格剖分技术6剖分中不规范的三角形出现的次数，有利于有限元分析时单元光滑、均匀，符合实际要求，从而提高逼近的精度。

所以，这也足以说服我们剖分时采用Delaunay 三角剖分。

性质3、局部重连（加点）。

当在点集中增加点时，只需做局部的改动，就能保证整体是Delaunay三角形。

为了说明这种局部改动，可以参考下图。

图 3.4性质4、局部重连（删点）。

当点集Delaunay三角化之后需要删去一个点，只需删去以此点为顶点的三角形，剩余的仍然是剩余顶点的Delaunay三角形。

以上的定理深刻的描述了Delaunay三角形的本质，也提供了构造Dlaunay 三角剖分的方法论。

我们在构造Delaunay三角剖分时，也就不走定义的那条路，即先构造Voronoi图，然后根据Delaunay三角形式Voronoi图的对偶图，从而得到Delaunay图，转而寻找符合空外接圆准则的三角剖分图，这样就保证了最终的确是Voronoi图。

而这样操作容易操作，算法的时间复杂度也会降低。

下面我们介绍构造点集Delaunay三角剖分的算法。

3.3 二维点集Delaunay三角剖分算法生成二维点集的Delaunay三角剖分算法经典的有两种：空洞算法和换边法。

783.3.1 空洞算法（增量算法）空洞算法的思想源于数学归纳法。

我们都知道数学归纳法的三大步骤：首先证明初始时命题成立，其次证明当目前命题成立时，下一步命题也成立，最后得到最终结果。

空洞算法也是这样：初始时，当给定4个点时可以生成2个Delaunay 三角形，然后增加一个点时，也可以生成Delaunay 三角形，最后一直加点，知道所有点集中的点被加入。

我们已经看到，增加点时由局部重连性质可以保证第二步确实成立。

那么问题只剩开始的四个点能够生成三角剖分吗？其实只要四点不共线就能保证，如下图。

图 3.53.3.2 换边法（局部变换法）换边法的原理仍然是空外接圆准则。

他要求先对所给点集做初始三角剖分，然后逐一检验每个三角形是否符合空外接圆准则，否则就四边形进行换边。

如下图：图 3.6则问题是，换掉对角线能保证局部是Delaunay 三角形吗？ 为此，我们需要以下引理。

9 引理1、圆上一条弦对某点所张开的角度，若点在圆内，则大于对应的圆周角；若在圆上，则等于对应的圆周角；若在圆外，则小于对应的圆周角。

回到原问题，我们假设ADC ∆的外接圆ABC O 包含点B ，则由上面引理，有 DAC DBC ∠>∠,ACD ABD ∠>∠。

那么在ABD ∆的外接圆ABD O 中，点C 肯定在其外，否则ABD ACD ∠≥∠， 同理，点A 在BCD ∆的外接圆BCD O 外，否则CBD CAD ∠≥∠。

从而换边之后符合Delaunay 三角剖分的空外接圆准则。

第四章算法设计本节中着重对边界提取和空洞算法详细介绍和进行分析。

4.1 边界提取。

我们的主要工作是确定图中满足灰度值要求象素群的边缘线或轮廓线。

由m*n个点组成的二值点阵图形在上一步我们已经实现。

这阶段主要来将其轮廓提取出来。

一般的，我们面对的二值图如下图所示：图 4.1若定义黑色和白色的公共边所连接构成的线为边缘线，这时我们可以看到边缘线将图片的轮廓完全描绘出来。