（水工结构工程专业论文）自然单元法在水工结构数值分析中的应用.pdf

摘要 自然单元法是近期出现的一种很有希望的无单元方法，它基于弗洛诺伊( v o r o n o i ) 图和德劳内( d e l a u n a y ) _ 三角形几何结构，采用希布逊( s i b s o n ) 或非希布逊( n o n - - s i b s o n ) 插值方法全域构造插值函数，既具有无网格方法和经典有限元方法的优点又克服了两者 的一些缺点，它可以像有限元法一样直接施加本质边界条件，不存在基于移动最小二乘 拟合的无网格方法不能直接施加本质边界条件的难题，而且自然单元法是无网格方法， 可以方便地处理有限元方法较难处理的一些问题，例如移动边界和大变形等问题，因此 受到很多学者的关注。 本文的主要工作如下： 1 本文较为详细的阐述了自然单元法的基本理论：自然邻点插值的几何基础德劳 内( d e l a u n a y ) - - 角形和弗洛诺伊( v o r o n o i ) 图，自然邻点希布逊( s i b s o n ) 插值和非希布逊 o n s i b s o n ) 插值的构造以及它们的性质，一维和二维自然邻点插值的特点等。 2 对于实现自然单元法所必需的背景积分网格形成即离散点的三角化、形函数的 计算、刚度矩阵的集成也作了必要的介绍，并给出了编程实现的流程图，此外，还简单 介绍了适用于本文程序的弹性理论。 3 根据上述理论，本文编制了弹性自然单元法c h 程序s g n e m l 0 ，利用本程序 对几个水工结构物进行分析，并和有限元分析结果或者解析解进行比较分析，结果表明， 自然单元法在水工结构数值分析上的应用是可行的。 关键词：自然单元法，c + + ，v o r o n o i 图，d d a u n a y 三角化，重力坝，地下厂房 a b s t r a c t t h en a t u r a le l e m e n tm e t h o d ( n e m ) i sak i n do fe l e m e n tf r e em e t h o d ( e f m ) n a t u r a l n e i g h b o ri n t e r p o l m i o no fn e m r e l i e so nc o n c e p t ss u c h 船v o r o n o id i a g r a m sa n dd e l a u n a y t e s s e l l a t i o n si nc o m p u t a t i o n a lg e o m e t r y , a n da d o p t ss i b s o no rn o n - s i b s o n i a ni n t e r p o l a t i o n m e t h o dt oc o n s t r u c ti n t e r p o l a n ti nt h ew h o l ea r e a n e mn o to n l yh o l d st h ea d v a n t a g e so f t h e c l a s s i c a lf e ma n de f mb u ta l s oo v e r c o m e st h es h o r t c o m i n go ft h e m ，i tc a r ld i r e c t l yi m p o s e s t h ee s s e n t i a lb o u n d a r yc o n d i t i o n sj u s tl i k ef e md o e s ，h e n c ea v o i dt h ed i f f i c u l t yo fi m p o s i t i o n o fe s s e n t i a lb o u n d a r yc o n d i t i o ni nu s i n go fe f mb a s e do nm o v i n gl e a s ts q u a r em e t h o d ，a n d b e c a u s en e mi so n eo fe f m i tc a nd e a l 、) l ，i t hs o m ep r o b l e m ss u c ha sm o v i n gb o u n d a r y p r o b l e ma n dl a r g ed e f o r m a t i o nt h a tf e mi sd i f f i c u l tt oh a n d l e t h e r e f o r e ，m a n ys c h o l a r sa r e p a y i n gm o r ea n dm o r ea t t e n t i o no ni t t h em a i nw o r k so ft h et h e s i sa r el i s t e da sf o l l o w s ： 1 t h et h e s i si n t r o d u c e si nd e t a i lt h eb a s i ct h e o r i e so fn e mi n c l u d i n gg e o m e t r i c f o u d a t i o no fi n t e r p o l a t i o n ：v o r o n o id i a g r a ma n dd e l a u n a yt r i a n g u l a t i o n ，t h ec o n s t r u c t i o no f n a t u r a ln e i g h b o u rs i b s o ni n t e r p o l a n ta n dn o n s i b s o n i a ni n t e r p o l a n ta n dt h ep r o p e r t i e so f t h e m ， t h ec h a r a c t e r i s t i c so f n a t u r a ln e i g h o u ri n t e r p o l a t i o ni no n ea n dt w od i m e n s i o n 2 t h et h e s i sr e p r e s e n t sn e c e s s a r yc o n t e n t ss u c ha st r i a n g u l a t i o no fn o d a ld i s c r e t i z a t i o n u s e da sb a c k g r o u n di n t e g r a t i o nm e s h ，n a t u r a ln e i g h b o u rs h a p ef u n c t i o nc a l c u l a t i o na n d a s s e m b l a g eo fs t i f f n e s sm a t r i xw h i c ha l la r ee s s e n t i a lt oi m p l e m e n t a t i o no fn e ma n do f f e r f l o w c h a r to fp r o g r a m m i n gr e a l i z a t i o nr i g h ta f t e re a c hr e p r e s e n t a t i o n f u r t h e r m o r e ，t h e r e l a v a n te l a s t i ct h e o r yi si n t r o d u c e dt o o 3 a c c o r d i n gt o t h ea b o v et h e o r y , t h ee l a s t i cn a t u r a le l e m e n ta n a l y s i sp r o g r a m s g n e m l 0i nc + + i sd e v e l o p e d a n dt h er e s u l t so faf e we x a m p l e su s i n gt h i sp r o g r a ma r e c o m p a r e dw i t l lf e m o ra n a l y t i c a ls o l u t i o n f r o mt h e s ew ec a nc o n c l u d et h a tt h en e mi sv a l i d t ot h en u m e r i c a la n a l y s i so fh y d m u l i cs t r u c t u r e k e y w o r d s ： n a t u r a le l e m e n tm e t h o d ，c + + ，v o r o n o id i a g r a m ，d e l a u n a yt r i a n g u l a t i o n ， g r a v i t yd a m ，u n d e r g r o u n dp l a n t i i 独立完成与诚信声明 本人的学位论文是在导师指导下独立撰写并完成的， 学位论文没有剽窃、抄袭等违反学术道德、学术规范的侵 权行为，否则本人愿意承担由此产生的一切法律责任和法 律后果，特此郑重声明。 学位论文作者( 签字) ：寸闺摒 0 6 年月f 7 日 第一章绪论 第一章绪论 1 1 常见的数值算法 1 1 1 有限差分和有限元法 有限差分方法( f d m ) 是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方 法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以泰勒 级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从 而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代 数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。 有限元方法是解决工程和数学物理问题的数值方法。可用有限元方法解决的典型问 题包括结构分析、热传导、流体流动、质量传输和电磁电位。用有限元方法求解一个问 题本质是要求解联立代数方程组，而不是直接求解微分方程。有限元的数值解给出连续 体中多个离散点的未知量的近似值。因此模拟物体的过程是将一个物体划分成由小的物 体或单元( 有限元) 组成的等价系统、这些单元通常与两个或更多的单位( 节点) 相互连接， 或与边界线或表面相互连接，这个过程叫做离散化。在有限元方法中，代替一次求解整 个物体，建立每一个有限单元的方程，并组合这些方程得出整个物体的解答。简言之， 结构问题的求解通常是指确定每个节点的位移和构成承载结构的每个单元内的应力，自 然单元法的解题过程和有限元法大体是一致的。对于有限元方法，其基本思路和解题步 骤亩归纳为： 1 建立积分方程，根据变分原理或方程余量与权函数正交化原理，建立与微分方程初 边值问题等价的积分表达式，这是有限元法的出发点。 2 区域单元剖分，根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相 互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作，这部分 工作量比较大，除了给计算单元和节点进行编号和确定相互之间的关系之外，还要 表示节点的位置坐标，同时还需要列出自然边界和本质边界的节点序号和相应的边 界值。 华北水利水电学院碗士论文 3 确定单元基函数，根据单元中节点数目及对近似解精度的要求，选择满足- 二定插值 条件的插值函数作为单元基函数。有限元方法中的基函数是在单元中选取的，由于 各单元具有规则的几何形状，在选取基函数时可遵循一定的法则。 4 单元分析：将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近；再 将近似函数代入积分方程，并对单元区域进行积分，可获得含有待定系数( 即单元中 各节点的参数值) 的代数方程组，称为单元有限元方程。 5 总体合成：在得出单元有限元方程之后，将区域中所有单元有限元方程按一定法则 进行累加，形成总体有限元方程。 6 。边界条件的处理：一般边界条件有三种形式，分为本质边界条件( e s s e n t i a lb o u n d a r y c o n d i t i o n ，狄里克雷( d i r i c h l e t ) 边界条件1 、自然边界条件( n a t u r a lb o u n d a r y c o n d i t i o n ，黎曼( r i e m a n n ) 边界条件) 、混合边界条件( m i x i n gb o u n d a r yc o n d i t i o n 柯 西( c a u c h y ) 边界条件) 。对于自然边界条件， 一般在积分表达式中可自动得到满足。 对于本质边界条件和混合边界条件，需按一定法则对总体有限元方程进行修正满足。 7 解有限元方程：根据边界条件修正的总体有限元方程组，是含所有待定未知量的方 程组，采用适当的数值计算方法求解，可求得各节点的函数值。 有限差分法直观，理论成熟，精度可选，易于编程，易于并行，f d m 的并行是当 前和将来应用的一个不错的方向。但是不规则区域处理繁琐，虽然网格生成可以使f d m 应用于不规则区域，但是对区域的连续性等要求较严。 有限元方法比有限差分优越的方面主要在于能够适应不规则区域。对于椭圆型方 程，如果区域规则，传统有限差分和有限元都能求解，在求解效率，收敛快慢、内存需 要等方面有限差分有优势，而有限元方法适合处理复杂区域，并且精度可选，缺憾在于 内存和计算量较大。 1 1 2 无网格法的发展现状 无网格法”刮又称为无单元法，继有限差分法，边界元法，有限元方法后，得到了 广泛的研究和重视，获到了巨大的发展。自从有限元法及电子计算机的广泛应用，计算 固体力学的发展是空前的，但计算固体力学仍面l 临着许多难以处理的问题，例如，在仿 真挤压和模压加工过程时，必须处理非常大的网格变形：在铸造的模拟计算中，固体和 液体之间的界面的变化是很重要的；在仿真材料失效过程时，需要模拟裂纹在任意方向 上的扩展：在高等材料的发展中，需要能跟踪相边界的增长和大量微观开裂等。有限元 2 第一章绪论 在处理这些问题时，就遇到了困难，其主要原因是网格的存在妨碍了处理与原始网格不 一致的不连续界面。在处理不连续界面时，给予网格的计算方法用的主要策略是在每一 计算步中重构网格，以保证网格线仍然与变形后的界面一致，这样不仅要碰到数值上的 困难，而且使计算精度下降和计算程序复杂，计算费用大大增加，有时还要因为变形太 大，使网格畸变而计算中断，而无网格法是部分或彻底取消网格的数值计算方法。 从2 0 世纪7 0 年代开始，已经提出了许多不用单元和网格的方法，有的以近似函数 命名，有的以离散化的方法命名，尚没有统一的名称。国际上用“e l e m e n tf r e em e t h o d ” 或“m e s h l e s sm e t h o d ”统称那些只取节点而不必生成单元网格的各种算法。历史上曾提 出过五花八门的无单元和无网格算法，名称也极不统一，如s a t l u r i 总结了几种主要算 法：光滑粒子水动学法( s p h - s m o o t h i n gp a r t i c l eh y d r o d y n a m i cm e t h o d ) 、有限点法 ( f p m f i n i t ep a r t i c l em e t h o d ) 、单位分解法( p u m p a r t i t i o no fu n i t y l 、无单元g a l e r k i n 法 ( e f g m e l e m e n tf r e eg a l e r k i nm e t h o d ) 、扩散单元法( d e m - - d i f f u s ee l e m e n tm e t h o d ) 、局 部p 无网格法( m l p g ) 、再生核质点法( r k p m - r e p r o d u c i n gk e r n e lp a r t i c l em e t h o d ) 、局部 边界积分法( l b i e - - l o c mb o u n d a r yi n t e g r a t i o nm e t h o d ) 。另外还有无网格局部伽辽金法 ( m e s h l e s sl o c a lp e 订o v - g a l e r k i nm e t h o d ) 、小波伽辽金法( w a v e l e t g a l e r k i nm e t h o d ) 、多尺 度再生核质点法( m u l t i s c a l er e p r o d u c i n gk e r n e lp a r t i c l em e t h o d ) 、小波质点法( w a v e l e t p a r t i c l em e t h o d ) 、移动最小二乘再生核法( m o v i n gl e a s t - s q u a r er e p r o d u c i n gk e r n e l m e t h o o ；多象限法( m u l t i q u a d r i e sm e t h o d ) 、h p 云团法( h p c l o u d sm e t h o d ) 、h p 无网格云 团法( h p m e s h l e s sc l o u d sm e t h o d ) 、无单元流形法( m a n i f o l dm e t h o d ) 、自然单元法( n a t u r a l e l e m e mm e t h o d ) 、有限覆盖无单元法( f i n i t e c o v e re l e m e n t f r e em e t h o d ) 。 在1 9 9 6 年，应用力学和工程计算机方法( c o m p u t e r m e t h o d s i n a p p l i e d m e c h a n i c sa n d e n g i n e e r i n g ) 杂志出版了一次会议的论文专辑，其中的第一篇论文把“无网格法 ( m e s h l e s s m e t h o d ) 作为标题来评述这些方法，这里也以无网格法作为统一的名称。 国内对无单元法的研究始于1 9 9 5 年，由于是消化和吸收已有的成果，所以发展很 快，其中清华大学的周维垣对无单元法进行了基本理论阐述，并结合数值流形法进行了 断裂力学的应用研究，此外还有陆明万和张雄等都为此作出了贡献，正是由于这些人的 引入及发展，使得国内很多人熟悉了这一数值计算的分支，并有越来越多的人投入到这 一领域里面来。 华北水利水电学院硕士论文 1 2 自然单元法的研究现状与进展 自然单元法是新近发展起来的无网格方法c a - v , b r a u n 和s a m b r i d g e 首先提出的一种 求解偏微分方程的算法，后来s u k u m a r 对二维空间中的自然单元法在解决弹性力学椭圆 型边界值问题的应用进行了比较深入的研究。自然单元法是一种多变量数据插值方法， 它是s i b s o n 提出的依赖于计算几何学的插值概念。最初自然邻点插值是用于数据插值和 地质现象模型，较之s h e p a r d 插值、移动最小二乘法，受到的重视还是很少的。自然单 元法利用需求解点自然邻接点的二次v o r o n o i 结构来构造近似位移函数，在局部的 d e l a u n a y 三角形予域上采用标准的g a l e r k i n 方法建立整体求解的平衡控制方程。同有限 元法比较，自然单元法的计算精度与同样网格划分的三节点或四节点单元有限元计算结 果相当，自然相邻插值函数比其他无网格法插值函数的计算速度快，同样由于自然邻点 插值在凸域的边界上的相邻点之间是严格线性的，所以边界的处理也相当简单。自然单 元法与其他数值方法的最根本区别在于其插值格式的不同，将自然邻点插值用于 g a l e r k i n 过程，就得到基于v o r o n o i 结构的自然单元g a l e r k i n 法，自然邻点插值有自然 邻点s i b s o n 插值和l a p l a c e 插值( n o n - s i b s o n 插值) 两种，l a p l a c e 插值比s i b s o n 插值在计 算上要简单的多，并且不论对凸的或非凸的区域都能精确旆加本质边界条件，以l a p l a c e 插值为试函数的自然单元法在数值实施上比以s i b s o n 插值为试函数的自然单元法简单。 自然单元法的背景积分网格为d e l a u n a y 三角形( 或四面体) ，不同于有限元法里面 的d e l a t m a y 网格需要对三角形的形状进行限制，自然单元法所采用的d e l a u n a y 背景积 分网格单元的形状、尺寸、角度可以是任意的。二维空间的有限元方法中，四边形的边 界逼近不如三角形单元精确，但三角形单元的精度往往难以达到要求。自然单元法代替 三角形单元有限元可以实现四边形单元精度和三角形单元边界逼近的结合。三维空间 中，目前六面体单元的自适应加密还没有较好的方法，在应用中多采用的是四面体单元。 在实际的工程计算中，采用自然单元法这样的无网格法可以完全采用软件自动进行 d e l a u n a y 算法的背景积分网格划分，大大提高了计算效率。将白适应算法与自然单元法 结合，可以解决岩土工程的复杂问题，如土体加筋的强度与变形、边坡问题等。一般来 讲，完全能够胜任复杂岩土及地下工程的数值分析任务。为了比较准确地描述岩土及地 下工程问题，其数值计算模型一般均取的比较大，相应的单元和结点数目较多，计算时 间也较长，因此，岩土及地下工程数值分析通用软件中需要解决的主要问题应该是在保 证求解精度和求解效率的条件下，把用户的前后处理工作减少到最少。但是由于采用标 4 第一章绪论 准的g a l e r k i n 方法系统平衡方程，需要在背景三角形积分网格里采用3 个以上的积分点 完成平衡方程的数值积分，加上计算形函数和寻找自然邻结点的时间消耗总的来说目前 自然单元法仍比相当精度的4 结点有限元法多一些。 自1 9 9 5 年b r a u n 和s a m b r i d g e 提出自然单元法概念以来，以v o r o n o i 图为几何基础 的数值方法在国外得到极大的关注。许多学者为此展开研究工作，并在不同的领域内得 到了一定的应用。s a m b r i d g e 和b r a u n 等将基于s i b s o n 插值的自然单元法应用于求解基 于高度不规则网格的偏微分方程，并应用于地球物理的研究，在他们的文献中，研究了 自然邻点插值在三维空间的数值计算方法，并指出了自然邻点插值可应用于有限元法、 地形数据处理等领域；s u k u m a r 等在一系列文章中把s i b s o n 插值格式的自然单元法应用 于二维弹性力学问题的研究，通过标准的g a l e r k i n 过程得到离散的线性方程组，给出了 形函数及其导数的计算方法和丰富的算例，并对在材料不连续模型和非凸体域( 存在裂 纹) 中，自然单元法的应用方式作了详细的讨论。在另一篇文章呻1 中，通过将s i b s o n 自 然邻点坐标嵌入b e m s t e i n b e z i e r 多项式，构造出具有c 1 的自然邻点插值，通过该插值 格式插值节点的位移及其导数，进而应用于求解四阶偏微分方程。还有篇文章凹1 中利 用l a p l a c e 插值格式的g a l e r k i n 过程，求解线弹性力学问题的椭圆型偏微分方程；基于 l a p l a c e 插值的自然单元法可以和线性有限元法实现无缝耦合，进而发挥各自的优势， 方便求解一些复杂问题。蔡永昌、朱合华、王建华基于自然邻点近似位移函数提出了 一种用于求解弹性力学平面问题的无网格局部p e t r o v g a l e r k i n 方法，该方法用加权残数 法推导控制方程，采用的是s i b s o n 插值格式，得到的系统矩阵是带状稀疏的。朱怀球、 吴江航等对s i b s o n 插值基函数的性质进行了研究，给出了基函数的一阶导数的一种数学 表达式及其数学性质，并将其应用于计算流体力学的研究中。c u e t o 等研究了自然单元 法中本质边界条件的施加方法，对自然单元g a l e r k i n 法的研究进展作了综述，描述了自 然单元g a l e r k i n 法的基本理论和数值计算方法，给出了在固体力学和流体力学中应用的 算例。最近，s u k u m a r 将v o r o n o i 单元和自然邻点插值的思想应用于求解在任意网格上 的扩散方程的有限差分方法，提出了v o r o n o i 单元有限差分法( v o r o n o ic e l lf i n i t e d i f f e r e n c em e t h o d ) 。王兆清与冯伟就当前自然单元法的发展概况及进展进行了综述。同 济大学的戴斌在他的硕士论文阳3 中从l a s s e r r e 凸多面体体积公式出发推导- 5 - - 维自然 邻结点坐标及其导数的算法，该算法通过直接计算二次v o r o n o i 胞元的体积求得自然邻 结点坐标及其导数，从理论上说它可以用于任意维空间中自然邻结点及其导数的计算。 华北水利水电学院硕士论文 作者并对自然单元法、l p 自然单元法等算法开发了应用软件“f n e m v l 0 ”。徐斌在他的 硕士论文阳3 中将无网格自然单元法用于弹塑性分析中，并编制了基于自然元法的弹塑性 数值分析程序。 1 3 本文的主要内容 自然单元法是新兴的无网格的数值算法，但由于发展时间较短，还有许多方面有待 发展和完善。根据国内外有关学者在有限元和自然单元法方面的研究成果，参照相关的 编程方法，做了以下的具体工作： 1 、对无网格自然单元法( s g n e m l o ) 所涉及的基本理论进行了较为详细的阐述。 2 、对于程序核心部分给出了详细的介绍，并附有计算程序的总体框架以及流程图。 3 、利用已编制的程序就几个算例加以应用，并和传统有限元方法的计算结果进行比较， 从而验证自然单元算法应用的可行性。 6 第二章自然单元法的基本理论 第二章自然单元法的基本理论 2 1v o r o n o i 图和d e l a u n a y 三角形 自然单元法是基于v o r o n o i 图和d e l a u n a y ( 德劳内) 三角形几何结构的数值计算方 法，这两个几何概念是形函数构造的几何基础1 0 - 1 2 。 2 1 1 概述 弗洛诺伊( v o m n o i ) 和它的对偶德劳内镶嵌( d e l a u n a yt e s s e l a t i o n ) 是一种定义不规 则点集最为基础和有效的几何作图。v o m n o i 图，又叫泰森多边形( t h i e s s e np o l y g o n ) 或狄 利克莱i 虱( d i r i c h l e td i a g r a m ) ，它是由一组由连接两邻点直线的垂直平分线组成的连续多 边形组成( 图2 - 1 中的a b c d ) 。n 个在平面上不同的点( 本图中共( 1 ，2 ，3 ，4 ，5 ) 5 个点) ，按照 最邻近原则划分平面，每个点与它的最近邻区域相关联。d e l a u n a y 三角形是由与相邻 v o r , o n o i 多边形共享一条边的相关点连接而成的三角形( 如三角形1 ：的边1 5 共享 v o r o n o i 多边形的动边) ，d d a u n a y 三角形的外接圆圆心是与三角形相关的v o m n o i 多边 形的一个顶点( 图中的空圆所示) ，d e l a u n a y 三角形是v o r o n o i 图的偶图。d e l a u n a y 三 角形产生的基本准则：任何一个d e l a u n a y 三角形的外接圆的内部不能包含其他任何点 ( d e l a u n a y1 9 3 4 ) 。l a w s o n ( 1 9 7 2 ) 提出了最大化最小角原则，每两个相邻的三角形构成凸 四边形的对角线，在相互交换后，六个内角的最小角不再增大。 j，一一、 图2 - 1v o r o n o i 图和d e h m 科三角形及夕卜接圃关系 f i g 2 - 1r e l a t i o n s h i p b e t w e e n v o r o n o i d i a g r a m a n d d e l a t m a y i r i a n g l e w i t h c i r c t m a e i r c l e 华北水利水电学院硕士论文 对给定的初始点集，有多种三角网剖分方式，其中d e l a u n a y 三角网具有以下特征： 1 不论从区域何处开始构网，最终都将得到一致的结果，d e l a u n a y 三角网是唯一的； 2 三角网的外边界构成了点集p 的凸多边形“外壳( c o n v e xh u l l ) ”； 3 点集中没有任何点在三角形的外接圆内部，反之，如果一个三角网满足此条件，那 么它就是d e l a u n a y 三角网，又可以表述为：在d e l a u n a y 三角形网中任一三角形的 外接圆范围内不会有其它点存在并与其通视，即空圆特性( e m p t yc i r c u r n c i r e l e c r i t e r i o n ) ： 4 在构网时，总是选择最邻近的点形成三角形并且不与约束线段相交，不会出现两个 三角形的边互相分割的现象； 5 形成的三角形网总是具有最优的形状特征，如果任意两个相邻三角形形成的凸四边 形的对角线可以互换的话，那么两个三角形6 个内角中最小的角度不会变大；如果 将三角网中的每个三角形的最小角进行升序排列，则d e l a u n a y 三角网的排列得到的 数值最大，从这个意义上讲，d e l a u n a y 三角网是“最接近于规则化”的三角网，这 个性质在离散点构网时很重要。 2 1 2 一阶v o r o n o i 图 为了说明在自然单元法中用到的v o r o n o i 图的概念，我们现在只考虑二维欧几里德 空间且，事实上这一理论在一般的n 维空间都适合。 假设在二维空间r 中有由m 个互异点构成的点集= ，n 2 ，, l m ) ，而集合的 一阶v o r o n o i 图指的是将平面分割成小区域i ( 这个小区域是封闭的、凸的，或者是无界 的) ，每个区域巧与点卿对应，使得在巧中的所有的点与 1 ( 最近邻接点) 的距离较与 其它点n ( 这里的n j 属于，但，不等于，) 的距离近，巧是与g l 距离更近点的聚集区域， 这里面就将霉称为点 1 的v o r o n o i 胞元，v o r o n o i 多边形巧的数学定义为： 乃= x r 2 ：d ( x ，x ，) 1 ) 的v o r o n o i 图来。对于自然邻点坐标，本文特别感兴趣的是七= 2 的 情形，即二阶v o r o n o i 图。 节点集合n 的二阶v o r o n o i 图是将平面分割成二阶v o r o n o i 胞元巧，每一个乃与 一个自然邻点对( 门，协) 对应，这样是与n i 最近点的集中区域，而n ，是第二最近 邻点，需要强调的是胞元只有当且仅当”，和是邻接点时才非空a 二阶v o r o n o i 胞 元数学定义如下： 巧= ( x r 2 ：d ( x ，x ) d ( x ，x j ) d ( x ，x 岸) ，v 足，j ( 2 2 ) 二阶v o r o n o i 胞元在自然邻点插值形函数的构造中具有重要的作用。 2 2 自然邻点插值 本节将介绍自然邻接插值的函数的构造和性质。自然单元法和其它无单元法和有限 元法的根本差别在于其形函数的构造方法，自然单元法中构造形函数的方法有l a p l a c e 插值( 又称n o n - s i b s o n 插值) 和s i b s o n 插值两种，本文中采用的是s i b s o n 插值，但也介 绍了l a p l a c e 插值的构造方法。以后如不特殊说明，所有的互异点集v 代表的都是凸域。 2 2 1s i b s o n 插值 自然邻点坐标的构造过程即是对引入图形当中的任意点x 的邻接关系加以量化的过 程，如图2 - 4 ，s i b s o n 采用二阶v o r o n o i 胞元的概念，引入了自然邻接点和自然邻接坐标， 通过定义自然邻接得出了新的一种度量邻接程度的方法，在图( a ) 中x 被放到了集合的 v o r o n o i 图中，将点x 连同点集合一起做v o r o n o i 图，则x 的自然邻接点就是与它在 新的三角形中若边的点，达到同一结果的直接方法是利用空圆准则，即如果x 位于三角 形d t ( n ，r t 。，n 。) 中，那么，吼是它的邻接点。自点x 和它的邻接点做中垂线，这 样就得到了二阶v o r o n o i 胞元( 封闭的多边形a b e d ) 。可以看到x 有四个自然邻接点，即 节点( 1 - 4 ) 。 令r ( x ) 为一阶v o r o n o i 胞元t ( x ) 的勒贝格测度( l e b e s g u em e a s u r e ，一维为长度， 1 0 第二章自然单元法的基本理论 二维为面积，三维为体积) ，令畸( ，= 1 , 4 ) 为点x 和点，的v o r o n o i 胞元重叠区域( 即二阶 v o r o n o i 胞元) 乃的测度，在二维情况下，测度为面积，因此我们可以得到这样的表达式 图2 _ 4 ( a ) 点x 乖l 原来的 “o r o n o i 图， f i 舀2 - 4 ( a ) 0 i j 西r l a i i d r 商盘蚪n a n d 蝎 点x 的一阶和二阶v o r o n o i 单元 l s t - o r d e r z l2 r i d - o r d e r v o r o n o ic e l l s 址j o u t x 彳( x ) = _ i r ( x ) 和4 ( x ) ；_ ( x ) ，我们可以定义点x 关于自然邻点，的自然邻点坐标，它是 一个比值，是这两个点的v o r o n o i 胞元的重叠的区域与x 的v o r o n o i 胞元面积的相除后 得到的。 力( x ) = 4 ( x ) a ( x ) ( 2 3 ) 在这里，从1 到月( n 为点x 的自然邻接点个数) ，并有： 彳( x ) = ：，4 ( x ) ( 2 4 ) 这四个区域如图2 - 4 ( b ) ，它们是二阶胞元，但它们的整体( 封闭的多边形a b c d ) 删 - - 阶的v o m n o i 胞元，参照图可以得到形函数萌( x ) 的表达式为： 舭) = 等 ( 2 5 ) 自然邻点坐标的导数通过求导等式( 2 5 ) 来得到： 舭) = 坐学u ：1 ，2 ) ( 2 6 ) 若考虑对于矢量函数u ( x ) ：qcr 2 斗r 2 的插值方案，可用作出形式如下的等式：+ u “( x ) = 办( x ) u ， ( 2 7 ) 华北水利水电学院硕士论文 这里u ，是竹个自然邻点处的节点位移矢量，而办( x ) 是与每个节点对应的形函数，在自 然邻点插值中，形函数c a x ) 被看成点x 在平面上的自然邻点坐标( n a t u r a ln e i g h b o r c o o r d i n a t e s ) ，这里值得注意的是由于形函数办( x ) 具有紧支特性，方程( 2 - 7 ) 是局部插值方 n o n - s i b s o n 插值的严格的定义：令= 如。，”：， 是d 维空间r 4 中一个不同 的点构成的点集，用正表示点x ，的v o r o n o i 胞元： 乃= x r 4 ：d ( x ，x ，) 0 办( x ) = 0 旧( x ) = 0 根据空圆准则，容易知对于形函数办( x ) ，要想在x 点处有非零的贡献，点x 必须位于有 ，点作为一个顶点的德劳内三角形外接圆内部，容易推得旃( x ) 的紧支域是凸壳c 日( ) 与所有通过，点德劳内三角形外接圆并集的交集。如图2 - 6 ( a ) ，一个单位正方形被离散 成2 5 ( 5 x s ) n 等距节点，a 点的紧支域如图2 - 6 ( b ) 所示，节点a 位于中心处，在这里 九( 扎) 呈现出整体最大值1 ，办( a ) 紧支域就是包含节点a 的所有的德劳内三角形( 此例 船 但 一 聊弛 x x x 华北水利水电学院硕士论文 中共8 个) 外接圆的并集。从图2 - 6 ( b ) 嗍，力( x ) 的紧支域的表面好比是橡胶皮被 拉伸来满足节点的数值。 图2 - 6 自然元法形函数的紧支域( a ) 节寺棚档和( b ) 节点a 的形函数以( x ) f 培2 - 6 s u p p o a f o r n e ms h a p e f m 蜊d o n ：( a ) n o d a l g r i d a n 1 嘞s h a p e f 硼d o l l 九( x ) f o r n o d e a 比较自然单元法插值函数的紧支域和其它一些广泛应用的面逼近方案的紧支域，发 现更多的自然单元法插值的本质和内在优越点。如谢巴德插值( s h e p a r di n t e r p o l a n t ) 和移 动最小二乘法插值( m o v i n gl e a s ts q u a r e sa p p r o x i m a t i o n s ) 都是基于与距离有关的权函 数，权函数一般来说都是各向同性的( 二维为圆，三维为球体) ，并且在确定半径的圆或 球内是非负的，与到点x 的距离成单调递减的关系，这些方法的基本原理是认为在x 点 处，与x 距离近的节点较那些较远的的节点提供了更大的权重。自然邻点插值采取了完 全不同的视角，在点x 的权重不是由在各个方向上同样的长度来决定的，而是通过合适 的空间维度的相应的勒贝格测度决定的，这样就允许了各向异性的紧支特性。在无单元 法中，大多的紧支尺寸都是基于距离的权重的，由于在高密度的节点区内，点x 的贡献 趋向于不成比例的偏大，所以对不规则排列的节点的处理是很困难的。在自然邻点插值 中，由于构造上的原因，节点分布和密度在给节点分配在点x 处的权重时就已经考虑在 内了。自然单元法中节点之间几何关系被称为空间邻接性( s p a t i a l a d j a c e n c y ) ，图2 7 表 现了空间邻接的特点。 邻接点的个数甩是位置x 及节点密度的函数，在d 维空间中，自然邻点的个数下限 为d + 1 个。图2 8 为一个均匀分布5 5 的栅格的凸壳内栉的变化，任何一个点至少有四 个邻接点，n 有两个值4 或6 ，而6 个邻接点的情况将在透镜形状的区域内获得。 第二章自然单元法的基本理论 图2 - 7 点b , c , d 晃羔a 的自然邻撬董，点e p 却不是 f 嘻2 - 7 v o r o n o i 删g h b o t a x b ，c ，a n d d a r e m i g h b o u r s o f a ，b u t e m d f a m n o t 图2 - 8 期缏4 栅格当中自然啕辘最的变化 f 嘻粥曲n o f 嘞删r 画g h t m f b r a 唧妇鲥d 2 2 3 5 光滑性 自然单元法形函数有很好的光滑性，自然邻接形函数在节点处是零阶可导且连续的 ( c o ) ，而在紧支域其它点处为无穷阶可导且连续( c 。，见形函数的计算式) 。如图2 3 所 示，力( 盖) 是x 的连续函数，当x 从任何方向趋近于盖，( 即：x j x ，) 时，旃( x ) 的连续性都 能够得到保证。除了节点外，所有点处的形函数力( x ) 的可导性是显而易见的，由于办( x ) 在紧支域内除了节点外的任何点都是光滑变化的，并且导数是无穷阶的。这一事实在第 四章推求自然单元法的形函数时仍然可以看到，点x 位于节点，的自然邻点外接圆上， 即x 规。此时只有一个圆内接三角形对形函数计算中的面积计算做出贡献，由于 c ，( x ) = v ( j = 1 ，3 ) ，子三角的面积和它的导数都为零。所彩。黼。7 ”“亡的耳黼7 一 华北水利水电学院硕士论文 阶或更高阶) 在x 点处都为零，因此说办( x ) 除在，点外是无穷阶可导且连续的。 2 2 4 一维自然邻点插值 如图2 - 9 ，长为上的一维杆件，将之离散为m 个不等距的点，根据v o r o n o i 图的定 义，v o r o n o i 图的顶点是任意两个相邻节点的中点，图中实圆代表节点，空圆代表v o r o n o i 图的顶点，在开集( o ，上) 上所有点都有两个自然邻点，而在边界上的点只有一个自然邻 123 m - 1m k l 一 l a ) 1 卜学一卜2 t = 工二i 二二 二 一号 n l ( 0 )鼍n i + 1 ( 1 ) 峭 ( b l 图2 _ 9 一维自然单元法形函数( a ) 物理空间，和( b ) 参考空间 f i g 2 - g n e m s h a p e f l l l l g t i o d s i n o n e d i m e n s i o n ：( a ) p h y s i c a ls p a c