吳繼華,王志非 (同濟大學(xué)土木工程學(xué)院,上海 200092)
傳統(tǒng)的有限單元法(FEM)理論成熟,并且有許多商業(yè)軟件可用,是目前在工程上應(yīng)用最廣泛的數(shù)值計算方法。但傳統(tǒng)FEM對單元劃分依賴度高,在處理大變形或裂紋擴展等問題時很難劃分出或自動生成合理的單元以滿足計算精度要求。無網(wǎng)格法克服了FEM對網(wǎng)格的依賴性,在求解涉及網(wǎng)格畸變、網(wǎng)格移動等問題時有明顯的優(yōu)勢。在無網(wǎng)格法的研究中,單位分解法的出現(xiàn)使得在近似函數(shù)中可以加入能夠反映待求解問題特性的函數(shù),從而也就催生了廣義有限單元法(GFEM)。GFEM采用分片任意高階多項式甚至是任意函數(shù)的線性組合作為逼近空間,因此GFEM可以通過提高廣義結(jié)點函數(shù)的階數(shù)和引入特定的局部逼近函數(shù),在不增加結(jié)點個數(shù)的前提下,有效地提高總體近似解的精度。近年來GFEM的理論基礎(chǔ)和應(yīng)用都得到了迅猛發(fā)展。
GFEM的概念在國外最早由Babuska I等人提出[1],但只有單位分解法[2]出現(xiàn)后才開始了系統(tǒng)的GFEM 研究,T·Strouboulisd等將FEM的插值函數(shù)作為單位分解法的單元分解函數(shù),可見GFEM是FEM和無網(wǎng)格法中的單位分解法的綜合[3]。Babuska I、C·A·Duarte、Barros FB等人發(fā)表了一系列論文[4~7],詳細闡述了GFEM的基本原理、收斂性和誤差估計等問題,并將GFEM應(yīng)用在三維結(jié)構(gòu)力學(xué)、動態(tài)裂紋發(fā)展、結(jié)構(gòu)非線性等領(lǐng)域,使GFEM的理論基礎(chǔ)和應(yīng)用都得到迅速發(fā)展。我國學(xué)者梁國平等吸收了流形方法的思想,將每個結(jié)點只有一個位移未知量的拉格朗日型插值空間推廣為每個結(jié)點有任意多個廣義位移的函數(shù)展開式,也提出了廣義結(jié)點有限元法[8],并從數(shù)學(xué)上論證了它的可行性。欒茂田[9][10]等人對上述方法進行了深入研究,給出了平面問題一階和二階GFEM的具體形式。程堯舜、辛婭云[11]對欒茂田等人的方法進行了改進,將坐標原點移到考察結(jié)點上,給出了插值函數(shù)的新形式,使得廣義自由度的意義得到明確,邊界條件也變得容易處理。程堯舜、李鳳嶺[12]詳細闡述了廣義自由度置零法、罰函數(shù)法等處理GFEM位移邊界條件的方法,并提出了一種新的處理GFEM位移邊界條件的方法,提高了計算精度。
GFEM相比較傳統(tǒng)的FEM而言,具有其自身的優(yōu)越性,但是GFEM每個結(jié)點有多個自由度,增加了計算量。程堯舜、李剛[13]針對這一缺陷提出了一種新型的GFEM,稱之為半彌散單元法(SDEM)。SDEM利用考察結(jié)點的鄰近結(jié)點的函數(shù)值來表示結(jié)點的廣義自由度,使得每個結(jié)點的自由度數(shù)和FEM相同,在系統(tǒng)自由度數(shù)不增加的前提下可以大大提高計算精度。S·Rajendran[14]結(jié)合FEM中形函數(shù)和無網(wǎng)格法中的單位分解法也提出了一種新的四邊形單元,這種方法和SDEM有一定的相似之處,也能使得自由度數(shù)降低,但是和SDEM相比計算量更大,四邊形單元在單元劃分時適應(yīng)性也比三角形單元低。黃書彪[15]通過引進更高階的插值函數(shù)對SDEM改進,提高了算法精度,并且提出了新的處理位移邊界條件的方法,通過改變已知位移邊界結(jié)點的鄰近點,使得在SDEM中可以使用FEM相同的方法來處理位移邊界條件。B·R·Zhang、S·Rajendran[16][17]將新提出的四邊形單元應(yīng)用于非線性問題和動力問題,得到了較好的結(jié)果。黃亮[18]將SDEM應(yīng)用于求解二維穩(wěn)態(tài)熱傳導(dǎo)問題中,數(shù)值算例的計算結(jié)果表明,在相同網(wǎng)格劃分下SDEM的精度遠高于FEM,體現(xiàn)了SDEM在計算中的優(yōu)勢。
設(shè)單元內(nèi)的真實場函數(shù)是一個二次完全多項式:
式中:a、b、c、d、e、f均為常數(shù)。那么結(jié)點i 處的場函數(shù)值和場函數(shù)的一階偏導(dǎo)數(shù)值如下:
下面分別用GFEM和無網(wǎng)格法的思想構(gòu)造局部近似場函數(shù)。首先考查GFEM中插值函數(shù)的特性。如圖1所示,三結(jié)點三角形單元的3個結(jié)點編碼按逆時針依次為i、j、m,為了減小條件數(shù),將坐標原點移到考察結(jié)點上,則未知場的二階單元插值函數(shù)表示為如下形式:
式中:x、y 為坐標;Fh為表示單元局部場函數(shù);di1、di2、di3為廣義自由度;Li為傳統(tǒng)有限單元法中三結(jié)點三角形單元的插值函數(shù)(見圖1)。
由式(4)可知,GFEM中每個結(jié)點的自由度為FEM的3倍,如果廣義自由度不做處理,計算量將大大增加。考察單元插值函數(shù)與真實場函數(shù)的聯(lián)系,將結(jié)點i的坐標代入式(4)可得到結(jié)點i處場函數(shù)值Fh=di1,因此di1可確定為:
圖1 三角形單元
將式(5)和式(6)代入式(4),并利用FEM中三結(jié)點三角形單元的行函數(shù)Li的性質(zhì)可得:
由式(7)可知,單元插值函數(shù)可以精確的表達二階真實場函數(shù)。其中參數(shù)di1即為該結(jié)點的場函數(shù)值,di2,di3分別該結(jié)點場函數(shù)對x、y的偏導(dǎo)數(shù)值的一半,即:
其次,從無網(wǎng)格法的思想出發(fā),構(gòu)造結(jié)點i附近的局部近似場函數(shù),取二階表達式如下:
將結(jié)點i的坐標和場函數(shù)值代入上式可得:
在結(jié)點i附近選取5個鄰近點,將鄰近點坐標和場函數(shù)值代入式(9)可得一組線性方程組:
式中,xij、yij是第j個鄰近點的位置坐標;Fij為第j個鄰近點函數(shù)值,j為鄰近點編號,j取1、2、3、4、5。參數(shù)中已由式(10)確定,另外5個參數(shù)可以由線性方程組(11)求得。因此參數(shù)已經(jīng)可以完全由結(jié)點的場函數(shù)值確定。
考察結(jié)點局部近似場函數(shù)與真實場函數(shù)的聯(lián)系,結(jié)點i處場函數(shù)值和一階導(dǎo)數(shù)值為:
由結(jié)點i局部近似場函數(shù)真實場函數(shù)求得的函數(shù)值和場函數(shù)的導(dǎo)數(shù)值應(yīng)該和由求得的值相同,因此:
比較(8)和(13)兩式,可以建立兩種局部近似場函數(shù)中廣義自由度的聯(lián)系:
根據(jù)上式,GFEM的廣義自由度可以由無網(wǎng)格法中廣義自由度來表示,而無網(wǎng)格法中廣義自由度可以由結(jié)點的場函數(shù)值來表達。因此,GFEM中的結(jié)點廣義自由度可以由結(jié)點鄰近點的場函數(shù)值解出。和一般的GFEM不同,SDEM正是利用這一聯(lián)系,在求解之前已經(jīng)把結(jié)點廣義的廣義自由度用結(jié)點鄰近點的場函數(shù)值表示出來,使得自由度大大減少。
GFEM就其本質(zhì)來講是在FEM的結(jié)點上增加新的自由度,再匹配以可以表現(xiàn)問題特性的插值函數(shù),這就使得在遇到復(fù)雜問題的時候GFEM可以不增加結(jié)點而只需要增加自由度以達到足夠的精度,避免了傳統(tǒng)有限單元法的網(wǎng)格重構(gòu)過程。相比FEM,廣義結(jié)點有限元具有其自身的優(yōu)越性,但是每個結(jié)點有多個自由度,這增加了數(shù)據(jù)的儲存量,增加了計算量。SDEM的提出給廣義有限單元法的研究帶來了新的思想,利用考察結(jié)點的鄰近結(jié)點的函數(shù)值來表示結(jié)點的廣義自由度,使得每個結(jié)點的自由度數(shù)和FEM相同,在系統(tǒng)自由度數(shù)不增加的前提下可以大大提高計算精度。
[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].力學(xué)進展,1995(4).
[9]欒茂田,田榮,楊慶.廣義結(jié)點有限元法[J].計算力學(xué)學(xué)報,2000(2).
[10]田榮,欒茂田,楊慶.高階形式廣義結(jié)點有限元法及其應(yīng)用[J].大連理工大學(xué)學(xué)報,2000(4).
[11]辛婭云.改進的廣義結(jié)點有限元法[D].上海:同濟大學(xué),2006.
[12]李鳳嶺.廣義結(jié)點有限元法中位移邊界條件實施研究[D].上海:同濟大學(xué),2009.
[13]李剛.半彌散單元法[D].上海:同濟大學(xué),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].上海:同濟大學(xué),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]黃亮.半彌散單元法在維穩(wěn)熱傳問題中的應(yīng)用[D].上海:同濟大學(xué),2013.