自然单元法与有限元法和无单元法的比较

2012-07-14琚宏昌刘平伟张胜利

水利水电科技进展 2012年2期

琚宏昌,刘平伟,张胜利

(广西工学院鹿山学院,广西柳州 545616)

其中

自然单元法(natural element method,NEM)是一种求解偏微分方程的新型数值方法。Braun等[1]于1995年首先将基于Sibson插值的自然单元法应用于求解高度不规则网格的偏微分方程,并指出这种数值方法可用于有限元法的应用领域,从此,自然单元法得到了极大的关注。

Sukumar等[2-3]将自然单元法应用于二维弹性力学问题的研究,详细讨论了自然单元法在裂纹体及材料的不连续体中的应用。Cueto等[4-6]研究了在自然单元法中施加本质边界条件的方法,给出了其在固体力学和流体力学中应用的算例。Gonzalez等[7]对自然单元法数值积分方案和自然单元法中的体积闭锁问题进行了研究。Yvonnet等[8]基于约束Delaunay结构对自然单元法做了扩展,解决了非凸边界上不满足线性插值的问题。Cho等[9]对二维接近不可压缩物体的大变形问题进行了分析。Calvo等[10]采用自然单元法分析了大应变超弹性问题。Alfaro等[11-12]采用基于α-shape的自然单元法对刚塑性模型的三维铝挤压成型过程进行了分析,并对自然单元法的计算精度特别是计算效率问题作了较为深入的研究。

国内学者在自然单元法研究和应用方面也做了很多工作。朱怀球[13]对Sibson插值基函数的性质进行了研究,给出了基函数的一阶导数数学表达式及其数学性质。朱合华等[14]将自然单元法应用于求解弹塑性问题。王兆清等[15]对自然单元法的发展概况及进展进行了综述。卢波等[16-18]对自然单元法数值积分方案进行了研究,对有限元法、无单元法及自然单元法进行了比较。余天堂[19]用自然单元法对加锚体进行了分析计算。李武等[20]分析了自然单元法数值积分产生误差的各种可能的原因,并提出使用蒙特-卡罗方法以提高数值积分的精度和效率。戴斌等[21]对自然单元法三维算法进行了初步研究。曾祥勇等[22]应用自然单元法分析了Winkler地基上薄板弯曲问题。蔡永昌等[23]采用局部Petrov-Galerkin法建立系统平衡方程以减少自然单元法的计算时间,提高该方法的计算效率。江涛等[24-26]应用自然单元法研究了裂纹与材料边界等问题。

本文对自然单元法、有限元法和无单元法的形函数性质进行了比较,并对自然单元法中尚待改进的问题进行了初步探讨。自然单元法兼具无网格方法和有限元法的优点,又克服了二者的不足,是一种发展前景广阔的求解偏微分方程的数值方法。

1 自然单元法形函数

自然单元法虽然在构造其形函数的过程中继承了无网格的思想,但其插值方法与无网格伽辽金方法及再生核质点方法等有所不同。自然单元法采用自然邻节点插值构造近似函数和试函数。

1.1 自然邻节点插值

自然邻节点插值的概念是建立在Voronoi图基础上的。图1(a)所示为平面7节点Voronoi划分及待插值点x的二次Voronoi结构。自然邻节点根据空圆规则确定,即若Delaunay三角形的外接圆包含待插值点x,则该三角形的3个顶点即为待插值点的自然邻节点。标号为1～4的4个节点构成了待插值点x的自然邻节点,待插值点 x与自然邻节点一起形成的Voronoi图称为待插值点 x的二次Voronoi结构,如图1(b)所示。

图1 Sibson和Non-sibson插值构造示意图

自然单元法形函数通过采用自然邻点插值思想的Sibson或Non-sibson插值构造,如图1(b)和(c)所示。

采用Sibson插值,自然单元法形函数公式为(参见图1(b))式中:Ai(x)为待插值点x的二次Voronoi结构与初始Voronoi结构中i节点所对应的Voronoi多边形相互重叠部分的面积,例如在图1(b)中,A3(x)为四边形cdef的面积,A(x)为四边形 abcd的面积;n为x的自然邻节点总数。

采用Non-sibson插值,自然单元法形函数公式为(参见图1(c))式中:sj(x)为与节点j关联的Voronoi边的长度;hj(x)为插值点 x到节点j的Voronoi边的垂距。

利用自然邻节点坐标 Φi(x),Sibson和Nonsibson的位移插值格式为式中:uh(x)是任一点x的位移;ui为节点i发生的位移。

Sibson插值需要构造二阶Voronoi图,而Nonsibson插值利用Voronoi单胞的边长和点到Voronoi边的距离构造插值基函数,使Non-sibson插值比Sibson插值的计算大为简化。Sibson插值在凸区域的边界是线性精确的,但是对于凹区域的边界,插值是不精确的,而Non-sibson插值则没有这个限制,因此采用Non-sibson插值,可以准确施加本质边界条件。

由自然邻节点插值的定义可知,节点的影响域为共享该节点的Delaunay三角形外接圆的并集,如图2(a)所示。节点 A处分别采用Sibson和Nonsibson插值方法构造的形函数如图2(b)和(c)所示。

图2 Sibson和Non-sibson形函数

1.2 自然单元法形函数的性质

文献[26]证明,自然单元法形函数满足正定性、插值性、单位分解条件和线性一致性:

此外,自然单元法形函数还具有在边界上满足线性插值,以及除了在节点处C0连续外,其余区域具有C∞光滑性的特点。

1.3 自然单元法基本原理

以二维弹性力学为例,说明自然单元法的基本原理。在边界 Γ包围的区域Ω,其平衡方程和边界条件可表示为

方程(5)的变分形式为

式中:t为面力向量;ε为应变矩阵;H10为一维Sobolev空间。

将式(3)代入式(6),并考虑到检验函数节点变量变分的任意性,可得离散系统节点变量的线性方程为

其中

式中:K,D,f分别为刚度矩阵、节点位移矩阵和节点力矩阵;E为弹性矩阵;Bi为自然单元法形函数的导数矩阵。

2 自然单元法、有限元法和无单元法形函数的统一性

从自然单元法基本原理可以看出,在构造自然单元法形函数的过程中不涉及矩阵的运算,而且形函数的构造简单,其导数的计算也相对简单,满足正定性、插值性、单位分解条件和线性一致性。

2.1 有限元法形函数的性质

有限元法是将整个求解域离散为有限个单元,然后在单元的基础上采用分片多项式插值构造插值函数。有限元法的形函数具有如下性质:①在节点上插值函数的值Ni(xj)=δij;②在单元中任一点各插值函数之和等于1,即。由此看出,有限元法形函数构成了单位分解。

2.2 无单元法形函数的特性

无单元法主要依靠形函数逼近来实现。无单元法按形函数逼近方式的不同,可分为3类:移动最小二乘逼近法(MLS)、再生核近似法(RKM)和单位分解法(PU)[16]。

对于平面问题,MLS插值可以精确再生多项式基线性组合函数,得到满足连续性条件的方程:RKM和MLS的形函数具有相近的性质,两者均构成了单位分解。

2.3 3种形函数统一于单位分解

有限元法利用单元离散求解域,单元之间共享节点,因而形成了有限元自身的覆盖。有限元形函数构造简单,但依赖于对整个求解域的网格剖分,网格质量直接影响计算精度。

无单元法和自然单元法通过节点的影响域实现对求解域的覆盖。无单元法形函数构造只需要节点的位置信息,构造过程中必须保证在每个节点的影响域内包含“足够多同时又尽量少”的节点。“足够多”是为了保证矩阵可逆,而“尽量少”是为了保证局部逼近的特性。节点影响域的大小是人为给定的,影响域的大小直接影响计算精度,而且形函数的构造复杂,计算量大。自然单元法形函数的构造仅取决于节点的位置信息,节点的影响域是由求解域内节点的Voronoi结构所限定的相邻关系决定的,形函数的构造形式简单,计算量小。

虽然构造方式不同,但有限元法、无单元法及自然单元法形函数均构成了单位分解。另外自然单元法和有限元法形函数是插值函数,而无单元法形函数是逼近函数。有限元法、无单元法及自然单元法的形函数如图3所示。

图3 规则分布节点上有限元法、无单元法及自然单元法生成的形函数

3 自然单元法存在的问题

3.1 弱形式的积分

以Delaunay三角形作为背景积分网格,在单位正方形区域边界上施加线性位移场,对标准Galerkin方法的控制方程弱形式数值积分进行位移分片检验,计算结果显示位移和应变相对误差范数分别可达到10-4和10-3,没有通过小片检验[2]。在其进行的平衡小片检验中位移和能量相对误差范数也分别可达到10-4和10-3,同样未能通过小片检验。自然单元法弱形式积分的这种误差,是由于自然单元法形函数的支撑域不能被分解为三角形,且自然单元法形函数不是多项式引起的。为了通过分片测试,在每个Ddaunay三角形内需要采用至少3个积分点,或采用高阶高斯积分来减少积分误差。但这种方法不仅对提高计算精度不明显,而且增加了计算的复杂性,造成求解效率略低于相当精度的四节点有限元。因此自然单元法弱形式的积分方案值得进一步研究。

3.2 节点应力恢复和边界应力精度的改善

自然单元法形函数为非多项式形式,在节点处C0连续而其余区域C∞光滑,其在节点处的导数不存在,这意味着自然单元法无法通过位移直接求得节点的应变及应力值。此外,其形函数在边界上相邻两节点间满足线性插值,这使得应力及应变在边界上为分段常数,应力及应变值在边界上的节点处不连续(跳跃)。因此,需要根据自然单元法算得的位移场恢复节点上的应力及应变值,并消除其在边界上的不连续(跳跃)现象,构造全域光滑的应力及应变场[16]。

文献[18]对基于Non-sibson插值的自然单元法的应力恢复和误差估计进行了研究。该文采用自然单元法求得的节点位移通过MLS拟合构造新的位移场,以此位移场计算节点处的应变和应力。再利用节点的恢复应力,采用自然邻点插值方法得到全域光滑的应力场。

文献[8]对基于稳定协调节点积分方案的自然单元法(nodal-NEM)的应力恢复和误差估计进行了研究,通过nodal-NEM方法可以方便地得到节点处的光滑应变。以节点处的光滑应变作为相应Voronoi单胞内的近似应变场,则Voronoi单胞内的近似应变场和应力场为常数,恢复应力场则通过节点处的应力由Sibson插值方法得到[26]。

4 结语

本文基于自然单元法、有限元法和无单元法形函数的构造方法和特点,进行了3种方法的比较,指出了3种形函数的单位分解共性以及自然单元法有待改进的问题。自然单元法是一种基于点集的Voronoi图和Delaunay三角化几何结构,以自然邻节点插值为试函数的一种新型求解偏微分方程的数值方法,其形函数满足正定性、插值性、单位分解条件和线性一致性,所构造的近似函数在区域边界上具有线性插值性,因而可以像有限元法一样方便地准确施加本质边界条件。由于自然单元法是无网格方法,可以方便地处理诸如复杂几何形状和裂纹扩展等问题,因此自然单元法既具有有限元法和无单元法的优点,又克服了两者的一些不足。自然单元法形函数的计算不涉及矩阵求逆,也不涉及权函数及相关参数的选择问题,可以方便实施。自然单元法兼具无单元法和有限元法的优点,是一种很有发展前途的数值方法。

[1]BRAUN J,SAMBRIDGE M.A numerical method for solving partial differential equations on highly irregular evolving grids[J].Nature,1995,376:655-660.

[2]SUKUMAR N,MORAN B,BELYTSCHKO T.The natural elements method in solid mechanics[J].International Journal for Numerical Methods in Engineering,1998,43:839-887.

[3]SUKUMAR N,MORAN B,SEMENOV A Y,et al.Natural neighbour Galerkin methods[J].International Journal for Numerical Methods in Engineering,2001,50:1-27.

[4]CUETO E,CALVO B,DOBLARE M.Modellingthreedimensional piece-wise homogeneous domains using the αshape-based natural elementsmethod[J].International Journal for Numerical Methods in Engineering,2002,54:871-897.

[5]CUETO E,CEGONINO J,CALVO B,et al.On the imposition of essential boundary conditions in natural neighbour Galerkin methods[J].Communications in Numerical Methods in Engineering,2003,19:361-376.

[6]CUETO E,SUKUMAR N,CALVO B,et al.Overview and recenet advances in natural neighbour Galerkin methodes[J].International Journal for Numerical Methods in Engineering,2003(4):307-384.

[7]GONZALEZ D,CUETO E,DOBLARE M.Volumetric locking in natural neighbour Galerkin methods[J].International Journal for Numerical Methods in Engineering,2004,61:611-632.

[8]YVONNET J,COFFIGNAL G,RYCKELYNCK D,et al.A simple error indicator for meshfree methods based on natural neighbors[J].Computers and Structures,2006,84:1301-1312.

[9]CHO J R,LEE H W.2-D large deformation analysis of nearly incompressible body by natural elementmethod[J].Computers&Structes,2006,84:293-304.

[10]CALVO B,MARTINEZ M A,Doblare M.On solving large strain hyperelastic problems with the natural element method[J].International Journal for Numerical Methods in Engineering,2005,62:159-185.

[11]ALFAR O I,BEL D,CUETO E,et al.Three-dimensional simulation of aluminium extrusion by the α-shape based Natural element method[J].Computer Methods in Applied Mechanics and Engineering,2006,195:4269-4286.

[12]ALFARO I,YVONNET J,CHINESTA F,et al.A Study on the performance of natural neighbour-based Galerkin methods[J].International Journal for Numerical Methods in Engineering,2007,71:1436-1465.

[13]朱怀球.Voronoi cells网格生成与C∞插值基函数的新软件开发及其在计算流体力学中的应用研究[D].北京:北京大学,2000.

[14]朱合华,杨宝红,蔡永昌,等.无网格自然单元法在弹塑性分析中的应用[J].岩土力学,2004,25(4):671-674.

[15]王兆清,冯伟.自然单元法研究进展[J].力学进展,2004,34(4):437-445.

[16]卢波,葛修润,孔祥礼.有限元法、无单元法及自然单元法之比较研究[J].岩石力学与工程学报,2005,24(5):780-786.

[17]沈振中,陈小虎,吴越建.求解堤坝渗流场的罚函数无单元法[J].河海大学学报:自然科学版,2008,36(1):44-48.

[18]张鑫,沈振中,徐力群.含孔岩体破裂过程的无单元法数值模拟[J].河海大学学报:自然科学版,2008,36(5):722-726.

[19]余天堂.加锚体自然单元法分析与程序设计[J].计算力学学报,2007,24(5):698-703.

[20]李武,姚文娟,朱合华,等.蒙特卡罗数值积分在自然单元法中的应用[J].岩土工程学报,2008,30(5):698-704.

[21]戴斌,王建华.自然单元法原理与三维算法实现[J].上海交通大学学报,2004,38(7):1223-1228.

[22]曾祥勇,朱爱军,邓安福.自然单元法在Winkler地基薄板计算中的应用[J].计算力学学报,2008,25(4):844-848.

[23]蔡永昌,朱合华,王建华.基于Voronoi结构的无网格局部Petrov-Galerkin方法[J].力学学报,2003,35(2):187-193.

[24]江涛,章青.线性强化材料弹塑性分析的自然单元法[J].力学季刊,2010,31(2):288-296.