Wachspress型多邊形有限元法積分方案

丁勝勇 邵國建

(河海大學工程力學系，南京 210098)

Wachspress型多邊形有限元法積分方案

丁勝勇 邵國建

(河海大學工程力學系，南京 210098)

摘 要:針對3種基于Wachspress插值的多邊形單元形函數，討論其各自的特點，選取其中最合理的形函數構造形式進行其偏導數公式推導，建立多邊形有限元的數值列式.對現行多邊形有限元法的積分方案進行探討，將傳統有限元法常用的Gauss積分運用到多邊形有限元法積分方案中，結合數值算例驗證了該積分方案的合理性和可行性，并給出了在不同類型網格中Gauss積分點數選取的建議.驗證在混合型網格中多邊形有限元法的合理性，對多邊形單元在疏密網格中作為一種過渡單元使用的可行性進行初步探討，為解決有限元計算中疏密網格的過渡問題提供了一種新的思路.

關鍵詞:多邊形單元;多邊形有限元;Wachspress插值;Gauss積分

有限單元法是固體力學中解決邊值問題最常用的一種數值方法.在二維問題中，傳統有限元方法常將單元劃分成常應變三結點三角形單元和雙線性四結點四邊形單元.多邊形單元作為一種新的單元形式，對于一些復雜幾何形狀的區域，具有區域適應性強和前處理方便等優點.在對多晶體材料的數值分析中，使用多邊形單元能有效降低分析單元數量，大大減少了計算所需時間［1］.在模擬裂紋擴展的問題中，應用多邊形有限元只需將裂紋穿過的單元分成多邊形單元，克服了傳統有限元法在裂紋區域需要不斷重新劃分網格的缺點［2］.

1975年，Wachspress首先提出了具有有理函數形式的凸多邊形單元形函數，稱為Wachspress插值法［3］，但因為其構造方法過于復雜而沒有得到重視.從20世紀90年代，工程數值計算和計算機圖形學等領域學者又對多邊形單元形函數構造方法做了大量的研究，取得了許多重要的成果［4－5］.文獻［1，6－8］提出了一些改進后的Wachspress型插值多邊形形函數.由于多邊形單元形函數是非多項式形式，以現有的數值積分方法，其數值積分結果的精度在理論上無法達到傳統有限元的計算精度，如何在現有理論框架下高效地獲得滿意的數值積分結果是一個值得進一步探討的問題.傳統有限元法常使用Gauss積分，并且根據其插值形函數的次數選取積分點數.為了將Gauss積分引入到多邊形有限元法，并探求多邊形有限元法積分點數的選取是否存在類似的規律，本文對3種基于Wachspress插值的多邊形形函數各自的特點進行了分析，并選取其中最合理的形式進行程序研制.對現有多邊形有限元法的積分方案進行進一步探討，結合算例對Gauss積分中積分點數的選擇進行數值計算對比，給出在不同類型網格中計算效率最佳的積分點數的選取建議.最后，對多邊形單元在疏密網格中作為一種過渡單元使用的可行性進行初步探索.

1 Wachspress多邊形有限元法

1.1 Wachspress插值形函數

由于多邊形單元的邊數是任意的，因此構造滿足單元位移協調性的多項式形式的位移插值函數變得異常困難.1975年，Wachspress利用映射幾何技術成功地在凸多邊形上構造出具有有理函數形式的多邊形單元形函數，但因為其構造過程繁瑣而未能得到廣泛的應用.對于具有n條邊的多邊形，其表達式為

式中，ωi為權函數，i=1，2，…，n.

Wachspress型插值多邊形形函數經過近十幾年不斷改進，已經成功地運用到多邊形有限單元法中.目前比較常見的Wachspress多邊形形函數有3種，圖1中的多邊形單元權函數表達式分別為

圖1 多邊形單元

式(2)～(4)可以在式(1)中通過數學變換獲得，但其顯式表達有一定區別.相對而言，ω(Ⅰ)型和ω(Ⅱ)型權函數的表達式較為簡單，且計算程序也便于編寫;ω(Ⅲ)型權函數的計算量會隨著單元邊數的增加而增加，程序編寫也最為繁瑣.但 ω(Ⅰ)，ω(Ⅱ)型權函數并不嚴格滿足Kronecker性質，在本點無意義，因此，本文采用ω(Ⅲ)型權函數進行分析計算.

在有限元單元勁度矩陣計算中，需要計算形函數對x，y的偏導數，而形函數偏導數的計算關鍵在于對其權函數偏導數的計算.對于式(4)有

1.2 積分方案

1.2.1 多邊形有限元法積分方案

多邊形單元由于邊數不定，無法直接在整個單元上進行積分，現行處理該積分的主要方法是將多邊形單元三角形化［4－5，9］，如圖 2 所示.然后在每個三角形上進行積分，其數學公式為

圖2 多邊形單元積分方案

現有數值積分方法在處理多項式函數積分時可獲得精確值，而對非多項式函數積分理論上并不完善.對于式(6)，由于被積函數f(x，y)中包含多邊形單元形函數表達式，而多邊形單元的形函數又為非多項式函數形式，因此，在理論上多邊形有限元計算結果的精度要受到影響.

1.2.2 三角形中的Gauss積分

針對平面三角形區域積分，對于每一個三角形區域，式(6)的積分計算式可表示為

假設積分區間是［a(x)，b(x)］，首先需要將［a(x)，b(x)］變換成［－1，1］，然后再應用 Gauss求積.下面給出最簡單的線性變換公式:

式(8)可運用Gauss積分公式，計算得到數值積分結果.

為了將Gauss積分應用到多邊形單元中，在式(7)的基礎上通過積分的換元法對積分域進行變換，可得到Gauss積分式為

式中，Cn，i為 Gauss 積分中的權函數;rn，i為 Gauss 積分中的積分點坐標.

式(9)中積分函數關于局部坐標L1，L2對稱，因此在L1和L2上應選擇相同積分點數，但該三角形中的Gauss積分點分布并不對稱，因此取直角邊為單位長度的等腰直角三角形，其不同Gauss積分點數分布如圖3所示.

圖3 三角形中Gauss積分點分布

本文中所提的多邊形有限元法積分點數為原多邊形單元經三角形化后在每個三角形中的數值積分點個數.

2 數值算例

2.1 小片試驗

考慮受單位分布力作用下的單向拉伸方板，板的邊長為單位1，其邊界條件如圖4所示，上部受均布拉力q=1 Pa.材料常數楊氏模量E=100 kPa，泊松比v=0.3，假設為平面應力狀態.

圖4 單向拉伸方板

小片的網格劃分如圖5所示，對5種不同的網格采用1.1節中ω(Ⅲ)型Wachspress形函數進行小片試驗，積分方案采用Gauss積分.表1為Wachspress形函數在5種網格中采用Gauss積分在不同積分點數時位移計算結果的相對誤差.

圖5 小片試驗的5種網格

由表1可以看出，在多邊形單元中運用Gauss積分的計算結果能夠通過小片試驗，因此該積分方案是可行的.在圖5(a)三角形網格和圖5(b)矩形網格中，多邊形有限元積分點數量對計算結果精度基本沒有影響，其計算結果精度與傳統有限元相當.Wachspress形函數在三角形單元上等價于三角形面積坐標，在矩形單元上等價于雙線性多項式形函數［10］，這里數值計算結果的精度得到了很好的體現.在圖5(c)、(d)、(e)的四邊形、五邊形和六邊形網格中，當積分點數為1×1時，計算結果相對誤差超過1%，誤差較大，而通過增加積分點數，可使數值計算結果精度得到很快提高并趨于穩定.

表1 Gauss積分相對誤差

2.2 懸臂梁

對如圖6所示的懸臂梁問題，取其長度L=8 m，高度h=1 m，單位厚度.x=0處為固定位移約束，在懸臂梁的右端施加P=1 kN的集中力作用.材料的彈性模量E=2 GPa，泊松比v=0.3.按照平面應力問題計算，利用本文Wachspress多邊形單元程序對圖7的網格分別進行計算.其中，矩形、四邊形、五邊形和六邊形網格單元數分別為256，256，240，266 個.取A(8，0)，B(2，0.5)兩點y方向位移值進行驗證，計算結果如圖8所示.在矩形、四邊形、五邊形和六邊形網格中，計算結果誤差都在1%以下.由此可見，使用本文中的積分方案，多邊形有限元計算結果可以滿足工程計算的要求.

圖6 懸臂梁計算模型

圖7 懸臂梁模型計算網格

圖8 A點和B點y方向位移值

理論上Gauss積分方法在積分點數足夠的情況下，針對多項式函數的積分是精確的.而多邊形單元形函數為非多項式形式，因此，多邊形有限元法的數值積分結果精度在理論上無法達到傳統有限元法的計算精度.圖8中在積分點數增加后計算結果未在一個固定值上穩定，而是一直在理論解值周圍作微小波動，這也進一步證明了這點.根據小片試驗和懸臂梁算例的計算結果，考慮到計算效率和計算結果精度2個方面，建議在多邊形有限元法使用Gauss積分時，在矩形網格中，每個三角化區域積分點數采用1×1，在任意四邊形、五邊形和六邊形網格中，每個三角化區域積分點數采用2×2.

2.3 混合型網格

本算例針對懸臂梁問題，網格劃分如圖9所示.采用以上給出的積分點數選取建議，對圖9中混合網格模型進行計算，考察A(8，0)和B(2，0.5)兩點在y方向的位移，計算結果相對誤差都在1%以下，這表明多邊形有限元使用在混合網格中也是合理的.

鑒于多邊形單元邊數可以是任意的，若在有限元分析中作為一種輔助單元使用，例如在圖9中作

圖9 混合網格模型

為疏密網格中間過渡的單元，不但可以使分析變得方便，而且能有效降低過渡區單元數量，提高計算效率.

3 結語

近年來，多邊形單元以其獨特的優點得到了越來越多學者的關注，也取得了很多重要的成果.由于多邊形單元形函數為非多項式，運用現行積分方法無法獲得其精確值，通過理論分析其積分點數的選取存在很大困難.本文將Gauss積分運用到多邊形有限元積分方案中，通過對不同積分點數下數值算例的結果進行對比分析，驗證了多邊形有限元在工程計算中的可行性，給出了在不同類型網格中Gauss積分點的選取建議，為多邊形有限元法的推廣運用提供了合理的依據.在疏密網格中，多邊形單元可用作過渡單元，為解決疏密網格的過渡問題提供了一種新的思路.

［1］范亞玲，張遠高，陸明萬.二維任意多邊形有限單元［J］.力學學報，1995，27(6):742－746.

Fan Yaling，Zhang Yuangao，Lu Mingwan.Arbitrarily polygonal 2-D finite elements［J］.Acta Mechanica Sinica，1995，27(6):742 －746.(in Chinese)

［2］唐旭海，鄭超，張建海.多邊形有限元法模擬裂紋擴展［C］//第17屆全國結構工程學術會議論文集(第Ⅰ冊).武漢，2008:446－450.

［3］Wachspress E L.A rational finite element basis［M］.New York:Academic Press，Inc，1975.

［4］王兆清.多邊形有限元研究進展［J］.力學進展，2006，36(3):344 －353.Wang Zhaoqing.Advances in polygonal finite element method［J］.Advances in Mechanics，2006，36(3):344－353.(in Chinese)

［5］Sukumar N，Malsch E A.Recent advances in the construction of polygonal finite element interpolants［J］.Archives of Computational Methods in Engineering，2006，13(1):129－163.

［6］Meyer M，Lee H，Barr A，et al.Generalized barycentric coordinates on irregular polygons［J］.Journal of Graphics Tools，2002，7(1):13 －22.

［7］Dasgupta G，Asce M.Interpolants within convex polygons:Wachspress'shape functions［J］.Journal of Aerospace Engineering，2003，16(1):1－8.

［8］Malsch E A，Dasgupta G.Interpolations for temperature distributions:a method for all non-concave polygons［J］.International Journal of Solids and Structures，2004，41(8):2165 －2188.

［9］蔡永昌，郭盛勇，楊健，等.多邊形單元的構造新方法及實現［J］.同濟大學學報:自然科學版，2009，37(7):883－887.

Cai Yongchang，Guo Shengyong，Yang Jian，et al.A new construction and implementation of polygonal elements［J］.Journal of Tongji University:Natural Science，2009，37(7):883 －887.(in Chinese)

［10］王兆清，李淑萍.多邊形有限單元形函數的比較研究［J］.應用力學學報，2007，24(4):604 －608.

Wang Zhaoqing，Li Shuping.Some remarks on shape functions of polygonal finite element［J］.Chinese Journal of Applied Mechanics，2007，24(4):604 －608.(in Chinese)

Integration scheme of Wachspress interpolation polygonal finite element method

Ding Shengyong Shao Guojian
(Department of Engineering Mechanics，Hohai University，Nanjing 210098，China)

Abstract:According to three kinds of Wachspress interpolation based shape function of polygonal element，the characteristics of their structural forms are discussed.The formula of partial derivative is derived by applying the best reasonable structural form and the polygonal FEM(finite element method)formulations are established.The current integration scheme in the polygonal FEM is studied.Based on this，the Gaussian integration，which is widely used in traditional FEM，is applied to the integration scheme of polygonal FEM.Numerical examples are presented to demonstrate the rationality and effectiveness of the proposed method.And the suggestion of selecting the number of the Gaussian integration points in different grids is given.The rationality of the polygonal FEM in mixed mesh is tested，and the feasibility of using the polygonal element as transition element in the coarse mesh and fine mesh is preliminarily discussed.This study provides a new idea to solve the transition problem of the coarse mesh and fine mesh in the FEM calculation.

Key words:polygonal element;polygonal finite element;Wachspress interpolation;Gauss integration

中圖分類號:TP391

A

1001－0505(2013)01-0216-05

doi:10.3969/j.issn.1001 －0505.2013.01.039

收稿日期:2012-04-21.

丁勝勇(1987—)，男，博士生;邵國建(聯系人)，男，博士，教授，博士生導師，gjshao@hhu.edu.cn.

基金項目:國家自然科學基金資助項目(50978083).

引文格式:丁勝勇，邵國建.Wachspress型多邊形有限元法積分方案［J］.東南大學學報:自然科學版，2013，43(1):216－220.［doi:10.3969/j.issn.1001 －0505.2013.01.039］