楊建, 謝偉, 張志偉

(1.西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安 710072; 2.液體火箭發(fā)動(dòng)機(jī)技術(shù)國(guó)防科技重點(diǎn)實(shí)驗(yàn)室, 陜西 西安 710100)

有限單元法(FEM)是一種通過(guò)將無(wú)限自由度問(wèn)題離散為有限自由度，從而得到近似結(jié)果的高效數(shù)值算法，該方法自問(wèn)世以來(lái)就在很多領(lǐng)域內(nèi)得到了廣泛應(yīng)用。目前對(duì)于三維有限元問(wèn)題，四面體和六面體單元具有原理簡(jiǎn)單、網(wǎng)格生成算法成熟等優(yōu)點(diǎn)，被絕大多數(shù)商業(yè)軟件所采用。然而在面對(duì)大型復(fù)雜結(jié)構(gòu)時(shí)，四面體網(wǎng)格計(jì)算精度差，六面體網(wǎng)格難以劃分等缺點(diǎn)也逐漸突出。

近年來(lái)，多面體單元由于其單元面的形狀和數(shù)量都具有任意性，為網(wǎng)格生成提供了巨大的靈活性等特點(diǎn)而受到廣泛關(guān)注[1-9]。然而因?yàn)殡y以構(gòu)造滿足單元協(xié)調(diào)性的多項(xiàng)式形式的形函數(shù)以及網(wǎng)格生成算法不夠成熟，與常規(guī)有限元相比其發(fā)展仍然處于起步階段。

劉桂榮教授和Nguyen-Thoi創(chuàng)立了一系列的光滑有限元法(S-FEM)[10]，由最初的光滑子單元域有限元法(CS-FEM)[11]，之后推廣到光滑節(jié)點(diǎn)域有限元法(NS-FEM)[12]、光滑面域有限元法(FS-FEM)[13]和光滑邊域有限元法(ES-FEM)[14]。這些方法都具有一個(gè)共同的特點(diǎn)，即不需要對(duì)形函數(shù)求導(dǎo)，降低了形函數(shù)必須連續(xù)的要求，這一特點(diǎn)可以有效地解決使用多面體網(wǎng)格遇到的問(wèn)題。

在以上這些S-FEM算法中，二維光滑邊域有限元法被證明具有最佳的穩(wěn)定性和精度[14]。所以本文在二維ES-FEM基礎(chǔ)上，建立了基于多面體單元的三維光滑邊域有限元法，劃分了三維ES-FEM在多面體單元中的的光滑域，推導(dǎo)了幾何矩陣和剛度矩陣。通過(guò)Matlab軟件編制的計(jì)算程序計(jì)算了彈性力學(xué)典型三維算例，并將結(jié)果與常規(guī)有限元結(jié)果等進(jìn)行了對(duì)比研究，證實(shí)了基于多面體單元的三維ES-FEM方法的可行性與可靠性。

1 ES-FEM-Poly的基本理論

在多面體網(wǎng)格中,ES-FEM-Poly以多面體單元的邊為基準(zhǔn),在與該邊相鄰的n個(gè)單元中劃分n個(gè)光滑子域,n個(gè)光滑子域組成了一個(gè)以該邊為基準(zhǔn)的光滑域。本文研究的多面體為包括四面體,六面體在內(nèi)的任意面體。圖1給出了以五面體單元為例的光滑子域的劃分方法,光滑子域由基準(zhǔn)邊的2個(gè)端點(diǎn)B、D,鄰接單元面的形心F、G,鄰接單元的體心H組成的五結(jié)點(diǎn)六面體。在多面體網(wǎng)格中,一個(gè)邊的鄰接單元數(shù)量不定,所以組成一個(gè)光滑域的光滑子域數(shù)量也不定,模型中邊的數(shù)量等于光滑域數(shù)量。

圖1 五面體單元光滑子域劃分方法

(1)

ES-FEM-Poly的形函數(shù)并不是關(guān)于相關(guān)結(jié)點(diǎn)位移的連續(xù)函數(shù),不需要導(dǎo)數(shù),針對(duì)這一特點(diǎn)使用線性插值法[15],直接計(jì)算相關(guān)結(jié)點(diǎn)在積分面上高斯積分點(diǎn)處的形函數(shù)值,將高斯積分點(diǎn)取為每個(gè)面的形心。以圖1所示的五面體為例,表1給出了對(duì)應(yīng)的形函數(shù)值表。表中第i行第j列的值代表第i個(gè)結(jié)點(diǎn)在第j個(gè)相關(guān)結(jié)點(diǎn)處的形函數(shù)值。

表1 形函數(shù)值表

對(duì)于一個(gè)擁有n個(gè)結(jié)點(diǎn)的多面體單元,每個(gè)結(jié)點(diǎn)在體心處的形函數(shù)值是1/n,在其他位置計(jì)算形函數(shù)的方法與五面體單元相同。

(2)

(3)

(4)

將(1)式代入(4)式中,可得

(5)

ES-FEM-Poly的總體平衡方程可由光滑Garlerkin弱化形式[13]得到

(6)

式中：b為體積力向量；t為邊界牽引力向量；D為彈性系數(shù)矩陣。

將(1)式和(5)式代入(6)式并簡(jiǎn)化可得

(7)

(8)

接下來(lái),將每個(gè)光滑域的剛度矩陣按照常規(guī)有限元方法組裝成整體剛度矩陣,然后求解平衡方程。

從以上推導(dǎo)過(guò)程可以看出,ES-FEM-Poly對(duì)光滑域的邊界進(jìn)行積分而無(wú)需對(duì)光滑域積分,從而無(wú)需像常規(guī)有限元那樣對(duì)坐標(biāo)進(jìn)行映射,它只要求多面體單元任意邊的2個(gè)結(jié)點(diǎn)與單元形心所組成的三角形面必須位于單元內(nèi)部,因此大大降低了對(duì)網(wǎng)格質(zhì)量的要求,即使是單元內(nèi)兩相鄰面內(nèi)角大于180°的畸形網(wǎng)格也同樣能滿足要求。

2 算例分析

本文利用Matlab軟件編寫(xiě)了ES-FEM-Poly算法程序,計(jì)算了空心球模型和懸臂梁模型,將計(jì)算結(jié)果與常規(guī)有限元四面體單元(FEM-T4)和六面體單元(FEM-H6)所得結(jié)果通過(guò)應(yīng)力相對(duì)誤差和能量相對(duì)誤差進(jìn)行比較。

應(yīng)力相對(duì)誤差

(9)

能量相對(duì)誤差

(10)

由于本文所用網(wǎng)格類型不同,且形狀不規(guī)則,引入單元特征邊長(zhǎng)來(lái)表征網(wǎng)格疏密程度。

(11)

式中：VΩ,Ne分別為整個(gè)模型的體積和單元個(gè)數(shù)。

2.1 空心球模型

空心球內(nèi)外半徑分別為a=5 m,b=10 m,彈性模量E=72 GPa,泊松比μ=0.33,因其對(duì)稱性,仿真計(jì)算時(shí)只取1/8模型,在其內(nèi)表面施加P=100 MPa的壓力,在其3個(gè)對(duì)稱面上分別施加x,y,z方向的位移約束,如圖2所示。

圖2 1/8空心球模型示意圖 圖3 球模型參考應(yīng)力云圖

選取1 016 379個(gè)六面體單元網(wǎng)格計(jì)算所得結(jié)果為參考解,圖3給出了參考應(yīng)力云圖,圖4給出了不同數(shù)量多面體單元對(duì)應(yīng)的應(yīng)力云圖。從圖中可以看出,應(yīng)力分布隨著單元數(shù)量的增加越來(lái)越接近參考解的應(yīng)力分布。

圖4 不同數(shù)量多面體單元的應(yīng)力云圖

圖5給出了不同數(shù)值方法下應(yīng)力相對(duì)誤差隨單邊特征邊長(zhǎng)的收斂曲線。從圖中可以看出,3條曲線比較接近,表明不同方法應(yīng)力誤差相差不大,但ES-FEM-Poly曲線的斜率明顯大于FEM-H6和FEM-T4的曲線斜率,表明：隨單元數(shù)量增加,lg(h)減小,ES-FEM-Poly計(jì)算結(jié)果更趨近參考解,表明其有更好的收斂性。當(dāng)h<0.5 m，即lg(h)<-0.3時(shí),ES-FEM-Poly的精度高于FEM-T4;當(dāng)h<0.38 m，即lg(h)<-0.42時(shí),ES-FEM-Poly的精度高于FEM-H6。

圖5 各算法的應(yīng)力相對(duì)誤差收斂曲線

圖6給出了不同數(shù)值方法下能量相對(duì)誤差隨單邊特征邊長(zhǎng)的收斂曲線。從圖中可以看出,ES-FEM-Poly和FEM-H6的精度相似,且高于FEM-T4;但ES-FEM-Poly曲線的斜率大于FEM-H6和FEM-T4的曲線斜率,表明ES-FEM-Poly的收斂性要比FEM-H6和FEM-T4好。

圖6 各算法的能量相對(duì)誤差收斂曲線

2.2 梁模型

如圖7所示,長(zhǎng)和寬為5 m,高為60 m的三維懸臂梁模型一端固支,另一端受到方向沿Y軸負(fù)方向的10 MPa均布剪力。彈性模量E=72 GPa,泊松比μ=0.33。

圖7 梁模型示意圖 圖8 梁模型參考應(yīng)力云圖

選取2 976 750個(gè)六面體單元網(wǎng)格計(jì)算所得結(jié)果為參考解,圖8給出了參考應(yīng)力云圖,圖9給出了不同數(shù)量多面體單元對(duì)應(yīng)的應(yīng)力云圖。從圖中可以看出,應(yīng)力分布隨著單元數(shù)量的增加越來(lái)越接近參考解的應(yīng)力分布。

圖9 不同數(shù)量多面體單元的應(yīng)力云圖

圖10給出了梁模型應(yīng)力相對(duì)誤差隨單邊特征邊長(zhǎng)的收斂曲線。從圖中可以看出,ES-FEM-Poly的曲線明顯低于FEM-T4和FEM-H6,表明ES-FEM-Poly的計(jì)算精度顯著高于FEM-T4和FEM-H6。當(dāng)h=0.5 m，即lg(h)=-0.3時(shí),ES-FEM-Poly的應(yīng)力相對(duì)誤差是FEM-H6的1/4,是FEM-T4的1/6。

圖10 各算法的應(yīng)力相對(duì)誤差收斂曲線

圖11給出了梁模型能量相對(duì)誤差隨單邊特征邊長(zhǎng)的收斂曲線。從圖中可以看出,使用多面體單元的ES-FEM-Poly的計(jì)算精度與FEM-H6相近,但顯著高于FEM-T4。當(dāng)h=0.5 m，即lg(h)=-0.3時(shí),ES-FEM-Poly的能量相對(duì)誤差是FEM-T4的1/24。

圖11 各算法的能量相對(duì)誤差收斂曲線

3 結(jié) 論

本文建立了基于多面體網(wǎng)格的光滑邊域有限元法,介紹了多面體單元的光滑域劃分方法、形函數(shù)的構(gòu)造以及幾何矩陣和剛度矩陣的推導(dǎo),利用Matlab軟件編程計(jì)算了空心球模型和梁模型,并將結(jié)果與常規(guī)有限元結(jié)果進(jìn)行對(duì)比研究,可以看出：①ES-FEM-Poly曲線的斜率明顯大于FEM-H6和FEM-T4的曲線斜率,表明隨著單元數(shù)量增加,ES-FEM-Poly計(jì)算結(jié)果更趨近參考解,有更好的收斂性；②ES-FEM-Poly曲線明顯低于FEM-H6和FEM-T4曲線,表明在單元特征邊長(zhǎng)一定時(shí),ES-FEM-Poly計(jì)算結(jié)果更接近參考解,精度更高。

綜上所述，基于多面體網(wǎng)格的三維光滑邊域有限元法,相比于常規(guī)有限元方法具有更好的精度和收斂性,可以用該方法對(duì)結(jié)構(gòu)復(fù)雜的三維彈性力學(xué)問(wèn)題進(jìn)行分析。