張 欣

0 引言

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

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

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

1 自然單元法與自然鄰點插值

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 算例

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

4 結(jié)語

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