自然單元法與有限元法和無單元法的比較

琚宏昌,劉平偉,張勝利

(廣西工學院鹿山學院,廣西柳州 545616)

其中

自然單元法(natural element method,NEM)是一種求解偏微分方程的新型數(shù)值方法。Braun等[1]于1995年首先將基于Sibson插值的自然單元法應用于求解高度不規(guī)則網(wǎng)格的偏微分方程,并指出這種數(shù)值方法可用于有限元法的應用領域,從此,自然單元法得到了極大的關注。

Sukumar等[2-3]將自然單元法應用于二維彈性力學問題的研究,詳細討論了自然單元法在裂紋體及材料的不連續(xù)體中的應用。Cueto等[4-6]研究了在自然單元法中施加本質邊界條件的方法,給出了其在固體力學和流體力學中應用的算例。Gonzalez等[7]對自然單元法數(shù)值積分方案和自然單元法中的體積閉鎖問題進行了研究。Yvonnet等[8]基于約束Delaunay結構對自然單元法做了擴展,解決了非凸邊界上不滿足線性插值的問題。Cho等[9]對二維接近不可壓縮物體的大變形問題進行了分析。Calvo等[10]采用自然單元法分析了大應變超彈性問題。Alfaro等[11-12]采用基于α-shape的自然單元法對剛塑性模型的三維鋁擠壓成型過程進行了分析,并對自然單元法的計算精度特別是計算效率問題作了較為深入的研究。

國內(nèi)學者在自然單元法研究和應用方面也做了很多工作。朱懷球[13]對Sibson插值基函數(shù)的性質進行了研究,給出了基函數(shù)的一階導數(shù)數(shù)學表達式及其數(shù)學性質。朱合華等[14]將自然單元法應用于求解彈塑性問題。王兆清等[15]對自然單元法的發(fā)展概況及進展進行了綜述。盧波等[16-18]對自然單元法數(shù)值積分方案進行了研究,對有限元法、無單元法及自然單元法進行了比較。余天堂[19]用自然單元法對加錨體進行了分析計算。李武等[20]分析了自然單元法數(shù)值積分產(chǎn)生誤差的各種可能的原因,并提出使用蒙特-卡羅方法以提高數(shù)值積分的精度和效率。戴斌等[21]對自然單元法三維算法進行了初步研究。曾祥勇等[22]應用自然單元法分析了Winkler地基上薄板彎曲問題。蔡永昌等[23]采用局部Petrov-Galerkin法建立系統(tǒng)平衡方程以減少自然單元法的計算時間,提高該方法的計算效率。江濤等[24-26]應用自然單元法研究了裂紋與材料邊界等問題。

本文對自然單元法、有限元法和無單元法的形函數(shù)性質進行了比較,并對自然單元法中尚待改進的問題進行了初步探討。自然單元法兼具無網(wǎng)格方法和有限元法的優(yōu)點,又克服了二者的不足,是一種發(fā)展前景廣闊的求解偏微分方程的數(shù)值方法。

1 自然單元法形函數(shù)

自然單元法雖然在構造其形函數(shù)的過程中繼承了無網(wǎng)格的思想,但其插值方法與無網(wǎng)格伽遼金方法及再生核質點方法等有所不同。自然單元法采用自然鄰節(jié)點插值構造近似函數(shù)和試函數(shù)。

1.1 自然鄰節(jié)點插值

自然鄰節(jié)點插值的概念是建立在Voronoi圖基礎上的。圖1(a)所示為平面7節(jié)點Voronoi劃分及待插值點x的二次Voronoi結構。自然鄰節(jié)點根據(jù)空圓規(guī)則確定,即若Delaunay三角形的外接圓包含待插值點x,則該三角形的3個頂點即為待插值點的自然鄰節(jié)點。標號為1～4的4個節(jié)點構成了待插值點x的自然鄰節(jié)點,待插值點 x與自然鄰節(jié)點一起形成的Voronoi圖稱為待插值點 x的二次Voronoi結構,如圖1(b)所示。

圖1 Sibson和Non-sibson插值構造示意圖

自然單元法形函數(shù)通過采用自然鄰點插值思想的Sibson或Non-sibson插值構造,如圖1(b)和(c)所示。

采用Sibson插值,自然單元法形函數(shù)公式為(參見圖1(b))式中:Ai(x)為待插值點x的二次Voronoi結構與初始Voronoi結構中i節(jié)點所對應的Voronoi多邊形相互重疊部分的面積,例如在圖1(b)中,A3(x)為四邊形cdef的面積,A(x)為四邊形 abcd的面積;n為x的自然鄰節(jié)點總數(shù)。

采用Non-sibson插值,自然單元法形函數(shù)公式為(參見圖1(c))式中:sj(x)為與節(jié)點j關聯(lián)的Voronoi邊的長度;hj(x)為插值點 x到節(jié)點j的Voronoi邊的垂距。

利用自然鄰節(jié)點坐標 Φi(x),Sibson和Nonsibson的位移插值格式為式中:uh(x)是任一點x的位移;ui為節(jié)點i發(fā)生的位移。

Sibson插值需要構造二階Voronoi圖,而Nonsibson插值利用Voronoi單胞的邊長和點到Voronoi邊的距離構造插值基函數(shù),使Non-sibson插值比Sibson插值的計算大為簡化。Sibson插值在凸區(qū)域的邊界是線性精確的,但是對于凹區(qū)域的邊界,插值是不精確的,而Non-sibson插值則沒有這個限制,因此采用Non-sibson插值,可以準確施加本質邊界條件。

由自然鄰節(jié)點插值的定義可知,節(jié)點的影響域為共享該節(jié)點的Delaunay三角形外接圓的并集,如圖2(a)所示。節(jié)點 A處分別采用Sibson和Nonsibson插值方法構造的形函數(shù)如圖2(b)和(c)所示。

圖2 Sibson和Non-sibson形函數(shù)

1.2 自然單元法形函數(shù)的性質

文獻[26]證明,自然單元法形函數(shù)滿足正定性、插值性、單位分解條件和線性一致性:

此外,自然單元法形函數(shù)還具有在邊界上滿足線性插值,以及除了在節(jié)點處C0連續(xù)外,其余區(qū)域具有C∞光滑性的特點。

1.3 自然單元法基本原理

以二維彈性力學為例,說明自然單元法的基本原理。在邊界 Γ包圍的區(qū)域Ω,其平衡方程和邊界條件可表示為

方程(5)的變分形式為

式中:t為面力向量;ε為應變矩陣;H10為一維Sobolev空間。

將式(3)代入式(6),并考慮到檢驗函數(shù)節(jié)點變量變分的任意性,可得離散系統(tǒng)節(jié)點變量的線性方程為

其中

式中:K,D,f分別為剛度矩陣、節(jié)點位移矩陣和節(jié)點力矩陣;E為彈性矩陣;Bi為自然單元法形函數(shù)的導數(shù)矩陣。

2 自然單元法、有限元法和無單元法形函數(shù)的統(tǒng)一性

從自然單元法基本原理可以看出,在構造自然單元法形函數(shù)的過程中不涉及矩陣的運算,而且形函數(shù)的構造簡單,其導數(shù)的計算也相對簡單,滿足正定性、插值性、單位分解條件和線性一致性。

2.1 有限元法形函數(shù)的性質

有限元法是將整個求解域離散為有限個單元,然后在單元的基礎上采用分片多項式插值構造插值函數(shù)。有限元法的形函數(shù)具有如下性質:①在節(jié)點上插值函數(shù)的值Ni(xj)=δij;②在單元中任一點各插值函數(shù)之和等于1,即。由此看出,有限元法形函數(shù)構成了單位分解。

2.2 無單元法形函數(shù)的特性

無單元法主要依靠形函數(shù)逼近來實現(xiàn)。無單元法按形函數(shù)逼近方式的不同,可分為3類:移動最小二乘逼近法(MLS)、再生核近似法(RKM)和單位分解法(PU)[16]。

對于平面問題,MLS插值可以精確再生多項式基線性組合函數(shù),得到滿足連續(xù)性條件的方程:RKM和MLS的形函數(shù)具有相近的性質,兩者均構成了單位分解。

2.3 3種形函數(shù)統(tǒng)一于單位分解

有限元法利用單元離散求解域,單元之間共享節(jié)點,因而形成了有限元自身的覆蓋。有限元形函數(shù)構造簡單,但依賴于對整個求解域的網(wǎng)格剖分,網(wǎng)格質量直接影響計算精度。

無單元法和自然單元法通過節(jié)點的影響域實現(xiàn)對求解域的覆蓋。無單元法形函數(shù)構造只需要節(jié)點的位置信息,構造過程中必須保證在每個節(jié)點的影響域內(nèi)包含“足夠多同時又盡量少”的節(jié)點?！白銐蚨唷笔菫榱吮ＷC矩陣可逆,而“盡量少”是為了保證局部逼近的特性。節(jié)點影響域的大小是人為給定的,影響域的大小直接影響計算精度,而且形函數(shù)的構造復雜,計算量大。自然單元法形函數(shù)的構造僅取決于節(jié)點的位置信息,節(jié)點的影響域是由求解域內(nèi)節(jié)點的Voronoi結構所限定的相鄰關系決定的,形函數(shù)的構造形式簡單,計算量小。

雖然構造方式不同,但有限元法、無單元法及自然單元法形函數(shù)均構成了單位分解。另外自然單元法和有限元法形函數(shù)是插值函數(shù),而無單元法形函數(shù)是逼近函數(shù)。有限元法、無單元法及自然單元法的形函數(shù)如圖3所示。

圖3 規(guī)則分布節(jié)點上有限元法、無單元法及自然單元法生成的形函數(shù)

3 自然單元法存在的問題

3.1 弱形式的積分

以Delaunay三角形作為背景積分網(wǎng)格,在單位正方形區(qū)域邊界上施加線性位移場,對標準Galerkin方法的控制方程弱形式數(shù)值積分進行位移分片檢驗,計算結果顯示位移和應變相對誤差范數(shù)分別可達到10-4和10-3,沒有通過小片檢驗[2]。在其進行的平衡小片檢驗中位移和能量相對誤差范數(shù)也分別可達到10-4和10-3,同樣未能通過小片檢驗。自然單元法弱形式積分的這種誤差,是由于自然單元法形函數(shù)的支撐域不能被分解為三角形,且自然單元法形函數(shù)不是多項式引起的。為了通過分片測試,在每個Ddaunay三角形內(nèi)需要采用至少3個積分點,或采用高階高斯積分來減少積分誤差。但這種方法不僅對提高計算精度不明顯,而且增加了計算的復雜性,造成求解效率略低于相當精度的四節(jié)點有限元。因此自然單元法弱形式的積分方案值得進一步研究。

3.2 節(jié)點應力恢復和邊界應力精度的改善

自然單元法形函數(shù)為非多項式形式,在節(jié)點處C0連續(xù)而其余區(qū)域C∞光滑,其在節(jié)點處的導數(shù)不存在,這意味著自然單元法無法通過位移直接求得節(jié)點的應變及應力值。此外,其形函數(shù)在邊界上相鄰兩節(jié)點間滿足線性插值,這使得應力及應變在邊界上為分段常數(shù),應力及應變值在邊界上的節(jié)點處不連續(xù)(跳躍)。因此,需要根據(jù)自然單元法算得的位移場恢復節(jié)點上的應力及應變值,并消除其在邊界上的不連續(xù)(跳躍)現(xiàn)象,構造全域光滑的應力及應變場[16]。

文獻[18]對基于Non-sibson插值的自然單元法的應力恢復和誤差估計進行了研究。該文采用自然單元法求得的節(jié)點位移通過MLS擬合構造新的位移場,以此位移場計算節(jié)點處的應變和應力。再利用節(jié)點的恢復應力,采用自然鄰點插值方法得到全域光滑的應力場。

文獻[8]對基于穩(wěn)定協(xié)調(diào)節(jié)點積分方案的自然單元法(nodal-NEM)的應力恢復和誤差估計進行了研究,通過nodal-NEM方法可以方便地得到節(jié)點處的光滑應變。以節(jié)點處的光滑應變作為相應Voronoi單胞內(nèi)的近似應變場,則Voronoi單胞內(nèi)的近似應變場和應力場為常數(shù),恢復應力場則通過節(jié)點處的應力由Sibson插值方法得到[26]。

4 結 語

本文基于自然單元法、有限元法和無單元法形函數(shù)的構造方法和特點,進行了3種方法的比較,指出了3種形函數(shù)的單位分解共性以及自然單元法有待改進的問題。自然單元法是一種基于點集的Voronoi圖和Delaunay三角化幾何結構,以自然鄰節(jié)點插值為試函數(shù)的一種新型求解偏微分方程的數(shù)值方法,其形函數(shù)滿足正定性、插值性、單位分解條件和線性一致性,所構造的近似函數(shù)在區(qū)域邊界上具有線性插值性,因而可以像有限元法一樣方便地準確施加本質邊界條件。由于自然單元法是無網(wǎng)格方法,可以方便地處理諸如復雜幾何形狀和裂紋擴展等問題,因此自然單元法既具有有限元法和無單元法的優(yōu)點,又克服了兩者的一些不足。自然單元法形函數(shù)的計算不涉及矩陣求逆,也不涉及權函數(shù)及相關參數(shù)的選擇問題,可以方便實施。自然單元法兼具無單元法和有限元法的優(yōu)點,是一種很有發(fā)展前途的數(shù)值方法。

[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網(wǎng)格生成與C∞插值基函數(shù)的新軟件開發(fā)及其在計算流體力學中的應用研究[D].北京:北京大學,2000.

[14]朱合華,楊寶紅,蔡永昌,等.無網(wǎng)格自然單元法在彈塑性分析中的應用[J].巖土力學,2004,25(4):671-674.

[15]王兆清,馮偉.自然單元法研究進展[J].力學進展,2004,34(4):437-445.

[16]盧波,葛修潤,孔祥禮.有限元法、無單元法及自然單元法之比較研究[J].巖石力學與工程學報,2005,24(5):780-786.

[17]沈振中,陳小虎,吳越建.求解堤壩滲流場的罰函數(shù)無單元法[J].河海大學學報:自然科學版,2008,36(1):44-48.

[18]張鑫,沈振中,徐力群.含孔巖體破裂過程的無單元法數(shù)值模擬[J].河海大學學報:自然科學版,2008,36(5):722-726.

[19]余天堂.加錨體自然單元法分析與程序設計[J].計算力學學報,2007,24(5):698-703.

[20]李武,姚文娟,朱合華,等.蒙特卡羅數(shù)值積分在自然單元法中的應用[J].巖土工程學報,2008,30(5):698-704.

[21]戴斌,王建華.自然單元法原理與三維算法實現(xiàn)[J].上海交通大學學報,2004,38(7):1223-1228.

[22]曾祥勇,朱愛軍,鄧安福.自然單元法在Winkler地基薄板計算中的應用[J].計算力學學報,2008,25(4):844-848.

[23]蔡永昌,朱合華,王建華.基于Voronoi結構的無網(wǎng)格局部Petrov-Galerkin方法[J].力學學報,2003,35(2):187-193.

[24]江濤,章青.線性強化材料彈塑性分析的自然單元法[J].力學季刊,2010,31(2):288-296.

[25]江濤,章青.自然單元法計算裂紋與材料邊界問題[J].應用力學學報,2009,26(4):690-694.

[26]江濤.固體力學自然單元法的若干問題研究[D].南京:河海大學,2009.