琚宏昌,劉平偉,張勝利

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

其中

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

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

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

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

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

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

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

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

圖1 Sibson和Non-sibson插值構(gòu)造示意圖

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

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

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

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

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

由自然鄰節(jié)點(diǎn)插值的定義可知,節(jié)點(diǎn)的影響域?yàn)楣蚕碓摴?jié)點(diǎn)的Delaunay三角形外接圓的并集,如圖2(a)所示。節(jié)點(diǎn) A處分別采用Sibson和Nonsibson插值方法構(gòu)造的形函數(shù)如圖2(b)和(c)所示。

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

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

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

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

1.3 自然單元法基本原理

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

方程(5)的變分形式為

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

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

其中

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

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

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

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

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

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

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

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

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

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

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

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

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

3 自然單元法存在的問題

3.1 弱形式的積分

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

3.2 節(jié)點(diǎn)應(yīng)力恢復(fù)和邊界應(yīng)力精度的改善

自然單元法形函數(shù)為非多項(xiàng)式形式,在節(jié)點(diǎn)處C0連續(xù)而其余區(qū)域C∞光滑,其在節(jié)點(diǎn)處的導(dǎo)數(shù)不存在,這意味著自然單元法無法通過位移直接求得節(jié)點(diǎn)的應(yīng)變及應(yīng)力值。此外,其形函數(shù)在邊界上相鄰兩節(jié)點(diǎn)間滿足線性插值,這使得應(yīng)力及應(yīng)變?cè)谶吔缟蠟榉侄纬?shù),應(yīng)力及應(yīng)變值在邊界上的節(jié)點(diǎn)處不連續(xù)(跳躍)。因此,需要根據(jù)自然單元法算得的位移場(chǎng)恢復(fù)節(jié)點(diǎn)上的應(yīng)力及應(yīng)變值,并消除其在邊界上的不連續(xù)(跳躍)現(xiàn)象,構(gòu)造全域光滑的應(yīng)力及應(yīng)變場(chǎng)[16]。

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

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

4 結(jié) 語

本文基于自然單元法、有限元法和無單元法形函數(shù)的構(gòu)造方法和特點(diǎn),進(jìn)行了3種方法的比較,指出了3種形函數(shù)的單位分解共性以及自然單元法有待改進(jìn)的問題。自然單元法是一種基于點(diǎn)集的Voronoi圖和Delaunay三角化幾何結(jié)構(gòu),以自然鄰節(jié)點(diǎn)插值為試函數(shù)的一種新型求解偏微分方程的數(shù)值方法,其形函數(shù)滿足正定性、插值性、單位分解條件和線性一致性,所構(gòu)造的近似函數(shù)在區(qū)域邊界上具有線性插值性,因而可以像有限元法一樣方便地準(zhǔn)確施加本質(zhì)邊界條件。由于自然單元法是無網(wǎng)格方法,可以方便地處理諸如復(fù)雜幾何形狀和裂紋擴(kuò)展等問題,因此自然單元法既具有有限元法和無單元法的優(yōu)點(diǎn),又克服了兩者的一些不足。自然單元法形函數(shù)的計(jì)算不涉及矩陣求逆,也不涉及權(quán)函數(shù)及相關(guān)參數(shù)的選擇問題,可以方便實(shí)施。自然單元法兼具無單元法和有限元法的優(yōu)點(diǎn),是一種很有發(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ā)及其在計(jì)算流體力學(xué)中的應(yīng)用研究[D].北京:北京大學(xué),2000.

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

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

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

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

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

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

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

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

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

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

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

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

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