廣義有限單元法研究進展

吳繼華，王志非 （同濟大學土木工程學院，上海 200092）

1 概 述

傳統的有限單元法（FEM）理論成熟，并且有許多商業軟件可用，是目前在工程上應用最廣泛的數值計算方法。但傳統FEM對單元劃分依賴度高，在處理大變形或裂紋擴展等問題時很難劃分出或自動生成合理的單元以滿足計算精度要求。無網格法克服了FEM對網格的依賴性，在求解涉及網格畸變、網格移動等問題時有明顯的優勢。在無網格法的研究中，單位分解法的出現使得在近似函數中可以加入能夠反映待求解問題特性的函數，從而也就催生了廣義有限單元法（GFEM）。GFEM采用分片任意高階多項式甚至是任意函數的線性組合作為逼近空間，因此GFEM可以通過提高廣義結點函數的階數和引入特定的局部逼近函數，在不增加結點個數的前提下，有效地提高總體近似解的精度。近年來GFEM的理論基礎和應用都得到了迅猛發展。

2 廣義有限單元法的研究歷史及現狀

GFEM的概念在國外最早由Babuska I等人提出[1]，但只有單位分解法[2]出現后才開始了系統的GFEM 研究，T·Strouboulisd等將FEM的插值函數作為單位分解法的單元分解函數，可見GFEM是FEM和無網格法中的單位分解法的綜合[3]。Babuska I、C·A·Duarte、Barros FB等人發表了一系列論文[4~7]，詳細闡述了GFEM的基本原理、收斂性和誤差估計等問題，并將GFEM應用在三維結構力學、動態裂紋發展、結構非線性等領域，使GFEM的理論基礎和應用都得到迅速發展。我國學者梁國平等吸收了流形方法的思想，將每個結點只有一個位移未知量的拉格朗日型插值空間推廣為每個結點有任意多個廣義位移的函數展開式，也提出了廣義結點有限元法[8]，并從數學上論證了它的可行性。欒茂田[9][10]等人對上述方法進行了深入研究，給出了平面問題一階和二階GFEM的具體形式。程堯舜、辛婭云[11]對欒茂田等人的方法進行了改進，將坐標原點移到考察結點上，給出了插值函數的新形式，使得廣義自由度的意義得到明確，邊界條件也變得容易處理。程堯舜、李鳳嶺[12]詳細闡述了廣義自由度置零法、罰函數法等處理GFEM位移邊界條件的方法，并提出了一種新的處理GFEM位移邊界條件的方法，提高了計算精度。

GFEM相比較傳統的FEM而言，具有其自身的優越性，但是GFEM每個結點有多個自由度，增加了計算量。程堯舜、李剛[13]針對這一缺陷提出了一種新型的GFEM，稱之為半彌散單元法（SDEM）。SDEM利用考察結點的鄰近結點的函數值來表示結點的廣義自由度，使得每個結點的自由度數和FEM相同，在系統自由度數不增加的前提下可以大大提高計算精度。S·Rajendran[14]結合FEM中形函數和無網格法中的單位分解法也提出了一種新的四邊形單元，這種方法和SDEM有一定的相似之處，也能使得自由度數降低，但是和SDEM相比計算量更大，四邊形單元在單元劃分時適應性也比三角形單元低。黃書彪[15]通過引進更高階的插值函數對SDEM改進，提高了算法精度，并且提出了新的處理位移邊界條件的方法，通過改變已知位移邊界結點的鄰近點，使得在SDEM中可以使用FEM相同的方法來處理位移邊界條件。B·R·Zhang、S·Rajendran[16][17]將新提出的四邊形單元應用于非線性問題和動力問題，得到了較好的結果。黃亮[18]將SDEM應用于求解二維穩態熱傳導問題中，數值算例的計算結果表明，在相同網格劃分下SDEM的精度遠高于FEM，體現了SDEM在計算中的優勢。

3 半彌散單元法的基本原理

設單元內的真實場函數是一個二次完全多項式：

式中：a、b、c、d、e、f均為常數。那么結點i 處的場函數值和場函數的一階偏導數值如下：

下面分別用GFEM和無網格法的思想構造局部近似場函數。首先考查GFEM中插值函數的特性。如圖1所示，三結點三角形單元的3個結點編碼按逆時針依次為i、j、m，為了減小條件數，將坐標原點移到考察結點上，則未知場的二階單元插值函數表示為如下形式：

式中：x、y 為坐標；Fh為表示單元局部場函數；di1、di2、di3為廣義自由度；Li為傳統有限單元法中三結點三角形單元的插值函數（見圖1）。

由式(4)可知，GFEM中每個結點的自由度為FEM的3倍，如果廣義自由度不做處理，計算量將大大增加。考察單元插值函數與真實場函數的聯系，將結點i的坐標代入式(4)可得到結點i處場函數值Fh=di1，因此di1可確定為：

圖1 三角形單元

將式(5)和式(6)代入式(4)，并利用FEM中三結點三角形單元的行函數Li的性質可得：

由式(7)可知，單元插值函數可以精確的表達二階真實場函數。其中參數di1即為該結點的場函數值，di2,di3分別該結點場函數對x、y的偏導數值的一半，即：

其次，從無網格法的思想出發，構造結點i附近的局部近似場函數，取二階表達式如下：

將結點i的坐標和場函數值代入上式可得：

在結點i附近選取5個鄰近點，將鄰近點坐標和場函數值代入式(9)可得一組線性方程組：

式中，xij、yij是第j個鄰近點的位置坐標；Fij為第j個鄰近點函數值，j為鄰近點編號，j取1、2、3、4、5。參數中已由式(10)確定，另外5個參數可以由線性方程組(11)求得。因此參數已經可以完全由結點的場函數值確定。

考察結點局部近似場函數與真實場函數的聯系，結點i處場函數值和一階導數值為：

由結點i局部近似場函數真實場函數求得的函數值和場函數的導數值應該和由求得的值相同，因此：

比較(8)和(13)兩式，可以建立兩種局部近似場函數中廣義自由度的聯系：

根據上式，GFEM的廣義自由度可以由無網格法中廣義自由度來表示，而無網格法中廣義自由度可以由結點的場函數值來表達。因此，GFEM中的結點廣義自由度可以由結點鄰近點的場函數值解出。和一般的GFEM不同，SDEM正是利用這一聯系，在求解之前已經把結點廣義的廣義自由度用結點鄰近點的場函數值表示出來，使得自由度大大減少。

4 結 論

GFEM就其本質來講是在FEM的結點上增加新的自由度，再匹配以可以表現問題特性的插值函數，這就使得在遇到復雜問題的時候GFEM可以不增加結點而只需要增加自由度以達到足夠的精度，避免了傳統有限單元法的網格重構過程。相比FEM，廣義結點有限元具有其自身的優越性，但是每個結點有多個自由度，這增加了數據的儲存量，增加了計算量。SDEM的提出給廣義有限單元法的研究帶來了新的思想，利用考察結點的鄰近結點的函數值來表示結點的廣義自由度，使得每個結點的自由度數和FEM相同，在系統自由度數不增加的前提下可以大大提高計算精度。

[1]Babushka I,Osborn.Generalized finite element methods:their performan-ce and their relation to mixed methods[J].SIAM Journal of Numerical Analysis,1983(3).

[2]Babuska I,Melenk JM.The partition of unity method[J].Int.J.Numer.Meth,Engng,1997(40).

[3]Strouboulis T,Babuska I,Copps K.Thedesign and analysisof the Genera-lized Finite Element Method[J].Comput.Methods Appl.Mech.Engrg.2000(181).

[4]Duarte CA,Babuska I,Oden JT.Generalized finite element methods for three-dimensional structural mechanics problems[J].Computer&Struc-tures,2000(77).

[5]Duarte CA,Hamzeh ON,Liszka TJ,et al.A generalized finite element method for thesimulation of three-dimensional dynamic crack propagation[J].Comput.Methods Appl.Mech.Engrg,2001(15-17).

[6]Strouboulis T,Babuska I,Hidajat R.The generalized finite element meth-od for Helmholtzequation:Theory,computation,and open problems[J].Comput.Methods Appl.Mech.Engrg,2006(37-40).

[7]Barros FB,Proenca SPB,de Barcellos CS.Generalized finiteelement met-hod in structural nonlinear analysis:ap-adaptivestrategy[J].Computational Mechanics,2004(2).

[8]梁國平,何江衡.廣義有限元方法——一類新的逼近空間[J].力學進展,1995(4).

[9]欒茂田,田榮,楊慶.廣義結點有限元法[J].計算力學學報,2000(2).

[10]田榮,欒茂田,楊慶.高階形式廣義結點有限元法及其應用[J].大連理工大學學報,2000(4).

[11]辛婭云.改進的廣義結點有限元法[D].上海:同濟大學,2006.

[12]李鳳嶺.廣義結點有限元法中位移邊界條件實施研究[D].上海：同濟大學,2009.

[13]李剛.半彌散單元法[D].上海：同濟大學，2007.

[14]S·Rajendran,B·R·Zhang.A“FE-meshfree”QUAD4 element based on partition of unity[J]Comput.Methods Appl.Mech.Engrg,2007(197).

[15]黃書彪.改進的半彌散單元法[D].上海：同濟大學，2008.

[16]B·R·Zhang,S·Rajendran.“FE-Meshfree”QUAD4 element for freevibrationanalysis[J]Comput.Methods Appl.Mech.Engrg.2008(197).

[17]S·Rajendran,B·R·Zhang.A partition of unity-based‘FE-meshfree’QUAD4 element for geometric non-linear analysis[J]Int.J.Numer.Meth.Engng,2010(82).

[18]黃亮.半彌散單元法在維穩熱傳問題中的應用[D].上海：同濟大學，2013.