丁勝勇 邵國建

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

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

丁勝勇 邵國建

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

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

關(guān)鍵詞:多邊形單元;多邊形有限元;Wachspress插值;Gauss積分

有限單元法是固體力學(xué)中解決邊值問題最常用的一種數(shù)值方法.在二維問題中，傳統(tǒng)有限元方法常將單元?jiǎng)澐殖沙?yīng)變?nèi)Y(jié)點(diǎn)三角形單元和雙線性四結(jié)點(diǎn)四邊形單元.多邊形單元作為一種新的單元形式，對(duì)于一些復(fù)雜幾何形狀的區(qū)域，具有區(qū)域適應(yīng)性強(qiáng)和前處理方便等優(yōu)點(diǎn).在對(duì)多晶體材料的數(shù)值分析中，使用多邊形單元能有效降低分析單元數(shù)量，大大減少了計(jì)算所需時(shí)間［1］.在模擬裂紋擴(kuò)展的問題中，應(yīng)用多邊形有限元只需將裂紋穿過的單元分成多邊形單元，克服了傳統(tǒng)有限元法在裂紋區(qū)域需要不斷重新劃分網(wǎng)格的缺點(diǎn)［2］.

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

1 Wachspress多邊形有限元法

1.1 Wachspress插值形函數(shù)

由于多邊形單元的邊數(shù)是任意的，因此構(gòu)造滿足單元位移協(xié)調(diào)性的多項(xiàng)式形式的位移插值函數(shù)變得異常困難.1975年，Wachspress利用映射幾何技術(shù)成功地在凸多邊形上構(gòu)造出具有有理函數(shù)形式的多邊形單元形函數(shù)，但因?yàn)槠錁?gòu)造過程繁瑣而未能得到廣泛的應(yīng)用.對(duì)于具有n條邊的多邊形，其表達(dá)式為

式中，ωi為權(quán)函數(shù)，i=1，2，…，n.

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

圖1 多邊形單元

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

在有限元單元?jiǎng)哦染仃囉?jì)算中，需要計(jì)算形函數(shù)對(duì)x，y的偏導(dǎo)數(shù)，而形函數(shù)偏導(dǎo)數(shù)的計(jì)算關(guān)鍵在于對(duì)其權(quán)函數(shù)偏導(dǎo)數(shù)的計(jì)算.對(duì)于式(4)有

1.2 積分方案

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

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

圖2 多邊形單元積分方案

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

1.2.2 三角形中的Gauss積分

針對(duì)平面三角形區(qū)域積分，對(duì)于每一個(gè)三角形區(qū)域，式(6)的積分計(jì)算式可表示為

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

式(8)可運(yùn)用Gauss積分公式，計(jì)算得到數(shù)值積分結(jié)果.

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

式中，Cn，i為 Gauss 積分中的權(quán)函數(shù);rn，i為 Gauss 積分中的積分點(diǎn)坐標(biāo).

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

圖3 三角形中Gauss積分點(diǎn)分布

本文中所提的多邊形有限元法積分點(diǎn)數(shù)為原多邊形單元經(jīng)三角形化后在每個(gè)三角形中的數(shù)值積分點(diǎn)個(gè)數(shù).

2 數(shù)值算例

2.1 小片試驗(yàn)

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

圖4 單向拉伸方板

小片的網(wǎng)格劃分如圖5所示，對(duì)5種不同的網(wǎng)格采用1.1節(jié)中ω(Ⅲ)型Wachspress形函數(shù)進(jìn)行小片試驗(yàn)，積分方案采用Gauss積分.表1為Wachspress形函數(shù)在5種網(wǎng)格中采用Gauss積分在不同積分點(diǎn)數(shù)時(shí)位移計(jì)算結(jié)果的相對(duì)誤差.

圖5 小片試驗(yàn)的5種網(wǎng)格

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

表1 Gauss積分相對(duì)誤差

2.2 懸臂梁

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

圖6 懸臂梁計(jì)算模型

圖7 懸臂梁模型計(jì)算網(wǎng)格

圖8 A點(diǎn)和B點(diǎn)y方向位移值

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

2.3 混合型網(wǎng)格

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

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

圖9 混合網(wǎng)格模型

為疏密網(wǎng)格中間過渡的單元，不但可以使分析變得方便，而且能有效降低過渡區(qū)單元數(shù)量，提高計(jì)算效率.

3 結(jié)語

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

［1］范亞玲，張遠(yuǎn)高，陸明萬.二維任意多邊形有限單元［J］.力學(xué)學(xué)報(bào)，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］唐旭海，鄭超，張建海.多邊形有限元法模擬裂紋擴(kuò)展［C］//第17屆全國結(jié)構(gòu)工程學(xué)術(shù)會(huì)議論文集(第Ⅰ冊).武漢，2008:446－450.

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

［4］王兆清.多邊形有限元研究進(jìn)展［J］.力學(xué)進(jìn)展，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］蔡永昌，郭盛勇，楊健，等.多邊形單元的構(gòu)造新方法及實(shí)現(xiàn)［J］.同濟(jì)大學(xué)學(xué)報(bào):自然科學(xué)版，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］王兆清，李淑萍.多邊形有限單元形函數(shù)的比較研究［J］.應(yīng)用力學(xué)學(xué)報(bào)，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

中圖分類號(hào):TP391

A

1001－0505(2013)01-0216-05

doi:10.3969/j.issn.1001 －0505.2013.01.039

收稿日期:2012-04-21.

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

基金項(xiàng)目:國家自然科學(xué)基金資助項(xiàng)目(50978083).

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