非凸邊界上自然單元法形函數(shù)計(jì)算方法研究

張 欣

0 引言

有限單元法在計(jì)算力學(xué)及其相關(guān)工程領(lǐng)域得到了廣泛的應(yīng)用，但由于其結(jié)點(diǎn)之間必須連接成相關(guān)的單元，因此在處理諸如裂紋擴(kuò)展、大變形等問題時因網(wǎng)格的變形和扭曲而影響解的精度甚至造成求解困難。近年來興起的無網(wǎng)格方法的主要優(yōu)點(diǎn)是不必將所研究的區(qū)域劃分為網(wǎng)格，但由于很多無網(wǎng)格方法的近似函數(shù)不具備插值性質(zhì)，因而有準(zhǔn)確施加本質(zhì)邊界條件的困難。

較晚出現(xiàn)的自然單元法(NEM)[1]，在凸區(qū)域的邊界上其形函數(shù)可以滿足 δ函數(shù)性質(zhì)并具有線性插值性，在非凸邊界上則需要采用α-shape方法[2]或約束自然單元法(C-NEM)[3]來計(jì)算形函數(shù)，但這兩種方法在實(shí)施時存在諸多限制和不便之處。

本文就非凸邊界上自然單元法形函數(shù)插值性能及其計(jì)算方法進(jìn)行了研究和探討，在揭示了造成非凸邊界上形函數(shù)不在邊界結(jié)點(diǎn)間線性變化的原因及其本質(zhì)的基礎(chǔ)上，提出了一種新的非凸邊界上自然單元法形函數(shù)計(jì)算方法，實(shí)現(xiàn)了形函數(shù)在邊界結(jié)點(diǎn)間的線性變化，對結(jié)點(diǎn)布置和區(qū)域邊界的凹凸程度并無限制，對各種形式的凸或非凸邊界計(jì)算方法是統(tǒng)一的。

1 自然單元法與自然鄰點(diǎn)插值

自然單元法是以自然鄰點(diǎn)插值作為試函數(shù)的求解偏微分方程的數(shù)值方法，Sukumar等在文獻(xiàn)[4]中將自然單元法應(yīng)用于二維彈性力學(xué)問題的研究。

Sibson[5]利用Voronoi圖和二階Voronoi單胞的概念定義了點(diǎn)x的自然鄰點(diǎn)坐標(biāo)。對圖1的二維情形，點(diǎn)x對于結(jié)點(diǎn)I的形函數(shù)為:

其中，AI為計(jì)算點(diǎn)x與其自然鄰點(diǎn)I對應(yīng)的二階Voronoi單胞的面積。如圖1所示，點(diǎn)x與其自然鄰點(diǎn)4對應(yīng)的二階Voronoi單胞是四邊形abfe所圍成的區(qū)域。

Non-Sibson[6]插值是較晚提出的一種基于Voronoi圖的插值方法，對于圖1的情形，點(diǎn)x對于結(jié)點(diǎn)I的non-Sibson插值形函數(shù)定義為:

其中，sJ(x)為計(jì)算點(diǎn)x與結(jié)點(diǎn)J相關(guān)的Voronoi邊的邊長; hJ(x)為計(jì)算點(diǎn)x到結(jié)點(diǎn)J的距離的一半;n為點(diǎn)x的自然鄰點(diǎn)個數(shù)。

在凸區(qū)域的邊界上，由于計(jì)算點(diǎn)與邊界結(jié)點(diǎn)相應(yīng)的二階Voronoi單胞的面積趨于無窮，因此在凸區(qū)域的邊界上計(jì)算點(diǎn)x的形函數(shù)僅當(dāng)其自然鄰點(diǎn)為邊界結(jié)點(diǎn)時才不為 0，且為計(jì)算點(diǎn)坐標(biāo)的線形函數(shù)。但是在非凸區(qū)域的邊界上，由于邊界結(jié)點(diǎn)的Voronoi單胞面積(三維時為體積)不再是無限的，因而插值函數(shù)在邊界節(jié)點(diǎn)間不再是線形變化的。

Cueto E等[2]使用建立點(diǎn)集α-shape的方法來構(gòu)建非凸區(qū)域的模型，以實(shí)現(xiàn)在非凸區(qū)域邊界上位移場是線性變化的，并稱這種方法為α-NEM。在計(jì)算力學(xué)中往往需要在場變量變化劇烈的地方加大布點(diǎn)密度，這種情況下點(diǎn)集的α-shape不一定能重構(gòu)區(qū)域的真實(shí)形狀，如此又引出了隨密度變化的α-shape，即α值隨局部的點(diǎn)的分布密度而變化。因此布點(diǎn)密度的變化不僅是計(jì)算精度的需要，也是模擬計(jì)算區(qū)域幾何模型的需要。

Yvonnet J等[3]通過引入可見性準(zhǔn)則來建立點(diǎn)集的約束Voronoi圖，并以此作為形函數(shù)計(jì)算的依據(jù)，實(shí)現(xiàn)了近似函數(shù)在非凸邊界上的線性插值，這種方法特別適用于具有裂紋邊界的非凸區(qū)域，稱為約束自然單元法(C-NEM)。

非凸邊界上自然單元法形函數(shù)計(jì)算的α-shape方法和C-NEM方法各有其優(yōu)點(diǎn)和不足。α-shape方法需要對結(jié)點(diǎn)布置加以控制，且對高度非凸邊界如裂紋尖端等情形是不適用的。C-NEM方法通過定義點(diǎn)集的約束Voronoi圖來限制相應(yīng)的鄰點(diǎn)關(guān)系，比較適合于具有裂紋或材料邊界的情形，但對于復(fù)雜邊界其計(jì)算復(fù)雜性和計(jì)算量都較大。本文提出一種計(jì)算非凸邊界上自然單元法形函數(shù)的統(tǒng)一方法，在明確定義邊界結(jié)點(diǎn)的基礎(chǔ)上，通過簡單的局部判定來限制相關(guān)的鄰點(diǎn)關(guān)系，從而實(shí)現(xiàn)在非凸邊界上形函數(shù)的線性變化和近似函數(shù)在邊界上的線性插值。相對于 α-NEM和C-NEM，將本文方法稱為邊界結(jié)點(diǎn)法(B-NEM)。

2 非凸邊界上自然單元法形函數(shù)計(jì)算的邊界結(jié)點(diǎn)法

影響形函數(shù)計(jì)算的因素是計(jì)算點(diǎn)的自然鄰點(diǎn)和相關(guān)結(jié)點(diǎn)的自然鄰點(diǎn)。為了使形函數(shù)在邊界結(jié)點(diǎn)間是線性變化的，邊界上的計(jì)算點(diǎn)x的Voronoi單胞在邊界線的外側(cè)應(yīng)該是無限的。由于非凸邊界上的計(jì)算點(diǎn)實(shí)際上是位于原點(diǎn)集凸殼的內(nèi)部，因而其Voronoi單胞總是有限的，而從另一個角度來看，引起這一現(xiàn)象的根本原因是非凸邊界上的計(jì)算點(diǎn)一定有位于其所在邊界線段外側(cè)的自然鄰點(diǎn)。

邊界結(jié)點(diǎn)法的基本思想是通過邊界結(jié)點(diǎn)確定邊界線段，且邊界線段是有向的，任一點(diǎn)x位于邊界線段ij的內(nèi)側(cè)是指按照i，j，x的循環(huán)次序所定義的三角形的有向面積為正值。通過限制位于邊界線段外側(cè)的結(jié)點(diǎn)成為計(jì)算點(diǎn)的自然鄰點(diǎn)，可以實(shí)現(xiàn)邊界上計(jì)算點(diǎn)的Voronoi單胞在邊界線段的外側(cè)是無限的，事實(shí)上通過局部邊界線段來限制相關(guān)點(diǎn)對的鄰點(diǎn)關(guān)系可以滿足可見性準(zhǔn)則的要求。

如圖 2所示，需要限制邊界線段 12外側(cè)的結(jié)點(diǎn) 3和結(jié)點(diǎn)4成為計(jì)算點(diǎn)x的鄰點(diǎn)，為了保證形函數(shù)計(jì)算的連續(xù)性，結(jié)點(diǎn) 3，4也不應(yīng)成為位于△125外接圓內(nèi)且位于邊界線段 12內(nèi)側(cè)的計(jì)算點(diǎn)的鄰點(diǎn)。注意到結(jié)點(diǎn) 3，4位于邊界線段 12的外側(cè)與計(jì)算點(diǎn)x位于邊界線段 31的外側(cè)是等效的，因此可以通過判斷計(jì)算點(diǎn)x位于邊界線段 31的內(nèi)側(cè)還是外側(cè)來決定結(jié)點(diǎn) 3，4能否成為其鄰點(diǎn)，換句話說邊界線段 12外側(cè)的結(jié)點(diǎn) 3，4不能以越過邊界的方式影響計(jì)算點(diǎn)x的形函數(shù)計(jì)算。

另一個需要處理的問題是計(jì)算區(qū)域的三角化應(yīng)該能準(zhǔn)確反映區(qū)域的邊界形狀，這可以通過允許空外接圓包含相應(yīng)邊界線段外側(cè)的結(jié)點(diǎn)來實(shí)現(xiàn)。圖 3中△123的外接圓包含位于邊界線段外側(cè)的結(jié)點(diǎn) 4，因此△123是有效的，但是應(yīng)該禁止△541成為有效的三角形，因?yàn)榻Y(jié)點(diǎn) 1位于邊界線段 45的外側(cè)。

3 算例

考慮圖 4情形，區(qū)域包括了兩類非凸邊界情形，一類是在裂紋分析中可能出現(xiàn)的高度非凸的邊界，另一類是一般的非凸區(qū)域，為了使邊界顯示的更加清楚，在裂紋邊界處有意地使裂紋兩側(cè)的邊界進(jìn)行了一定程度的分離，事實(shí)上裂紋兩側(cè)的邊界結(jié)點(diǎn)位置也可以是重合的，但需要進(jìn)行不同的結(jié)點(diǎn)編號。兩個不同結(jié)點(diǎn)1，13的形函數(shù)計(jì)算結(jié)果見圖 5，限于篇幅未示出形函數(shù)的導(dǎo)數(shù)的計(jì)算結(jié)果。

4 結(jié)語

為實(shí)現(xiàn)非凸邊界上自然單元法形函數(shù)在邊界結(jié)點(diǎn)間的線性變化，本文提出通過邊界結(jié)點(diǎn)確定相應(yīng)邊界線段的內(nèi)外側(cè)，以簡單的局部判斷避免在三角化時出現(xiàn)跨越邊界的三角形，在形函數(shù)計(jì)算時避免結(jié)點(diǎn)以跨越邊界的方式影響形函數(shù)的計(jì)算，實(shí)現(xiàn)了形函數(shù)在邊界結(jié)點(diǎn)間的線性變化。通過對多種形式的邊界進(jìn)行實(shí)際計(jì)算，結(jié)果顯示無論對凸邊界、一般非凸邊界、裂紋邊界、材料邊界以及內(nèi)部結(jié)點(diǎn)比較接近于邊界時的情形，計(jì)算所得形函數(shù)在邊界線段上是線性變化的，且采用的方法或判斷準(zhǔn)則是一致的。因此本文方法具有方便實(shí)用、適用面廣和計(jì)算方法統(tǒng)一等特點(diǎn)，為自然單元法的進(jìn)一步應(yīng)用奠定了堅(jiān)實(shí)的基礎(chǔ)。

[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] Cueto E，DoblaréM，Gracia L.Im posing essential boundary conditions in the natural neighbour Galerkin method by means of density-scaledα-shapes[J].International Journal for Numerical Methods in Engineering，2000(49):519-546.

[3] Yvonnet J，Ryckelynck D，Lorong P，et al..A new extension of the natural element method for non-convex and discontinuous problems:the constrained natural element method(C-NEM) [J].International Journal for Numerical Methods in Engineering，2004(60):1451-1474.

[4] 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.

[5] Sibson R.A Vector Identity for the Dirichlet Tesselation[J]. Mathematical Proceedings of the Cambridge Philosophical Society，1980(87):151-155.

[6] Belikov V V，Ivanov V D，Kontorovich V K，et al..The non-Sibsonian interpolation.A newmethod of interpolation of the values ofa function on an arbitrary set of points[J].Computational Mathematics and Mathematical Physics，1997(37):9-15.