非凸邊界上自然單元法形函數計算方法研究

張 欣

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

2 非凸邊界上自然單元法形函數計算的邊界結點法

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

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

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

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

3 算例

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

4 結語

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

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