賈 程，陳卉卉

(鹽城工學(xué)院土木工程學(xué)院，江蘇鹽城224051)

0 引言

目前，有限單元法已經(jīng)被廣泛地應(yīng)用于解決各種工程問(wèn)題，但該方法仍然存在一些缺陷.傳統(tǒng)的有限單元法通常采用協(xié)調(diào)模型，使得剛度矩陣“過(guò)剛”，能導(dǎo)致計(jì)算所得位移和應(yīng)力的精度下降，尤其對(duì)線性三角形單元而言更為明顯;能使得等參單元對(duì)網(wǎng)格畸變敏感.學(xué)者們提出了許多提高有限元精度的方法，如非協(xié)調(diào)模式法［1-2］、四邊形面積坐標(biāo)法［3］、無(wú)網(wǎng)格法［4-5］等.近些年來(lái)，單位分解法［5-6］引起了研究者的注意.單位分解有限元法較普通有限元的優(yōu)點(diǎn)是不需增加額外的節(jié)點(diǎn)，就能構(gòu)造高階的場(chǎng)函數(shù)，這為創(chuàng)建高效的單元格式提供了一種方法.另外，它能夠自由地選擇局部近似函數(shù)，這使其便于分析復(fù)雜問(wèn)題.

筆者采用線性三角形“單元”形函數(shù)作單位分解函數(shù)，利用最小二乘法建立局部近似場(chǎng)函數(shù)進(jìn)而構(gòu)造有限元無(wú)網(wǎng)格耦合三角形單元.其形函數(shù)具有Kronecker delta性質(zhì).通過(guò)分析二維固體的自由振動(dòng)和強(qiáng)迫振動(dòng)響應(yīng)，表明該單元沒(méi)有虛假的零能模式，計(jì)算結(jié)果的精度優(yōu)于線性三角形單元和線性四邊形等參元.

1 有限元無(wú)網(wǎng)格耦合三角形單元公式

1.1 有限元無(wú)網(wǎng)格耦合三角形單元的形函數(shù)

對(duì)于二維線彈性問(wèn)題，采用三角形網(wǎng)格離散問(wèn)題域，如圖1所示.單元內(nèi)任一點(diǎn)的位移可表示成:

其中:N'= [N1，N2，N3]是傳統(tǒng)的線性三角形單元的形函數(shù)矩陣［7］;ue={u1(x，y)u2(x，y)u3(x，y)}T是局部節(jié)點(diǎn)位移函數(shù).在該節(jié)點(diǎn)處等于其位移值，即 ui(xi，yi)=ui，i=1，2，3.局部節(jié)點(diǎn)位移函數(shù)ui(x，y)利用Rajendran等［6］使用的最小二乘點(diǎn)插值無(wú)網(wǎng)格法(LSPIM)，由i點(diǎn)的支持域內(nèi)的節(jié)點(diǎn)值擬合得到:

其中，Φi是LSPIM法的關(guān)于節(jié)點(diǎn)i的形函數(shù)矩陣，ui是節(jié)點(diǎn)i支持域內(nèi)節(jié)點(diǎn)的位移參數(shù)向量，N為節(jié)點(diǎn)i支持域內(nèi)的所有節(jié)點(diǎn)數(shù).一個(gè)單元所有節(jié)點(diǎn)支持域的合集構(gòu)成了一個(gè)單元的支持域Ω.

節(jié)點(diǎn)支持域和單元支持域的定義如圖1所示:

將式(2)帶入式(1)，得

從方程(3)中，可以得到該單元的形函數(shù)矩陣，設(shè)單元支持域Ω內(nèi)的節(jié)點(diǎn)數(shù)為M:

圖1 支持域的定義Fig.1 Definition of the support.

1.2 局部節(jié)點(diǎn)位移函數(shù)的形函數(shù)Φi

局部節(jié)點(diǎn)位移可以寫成下列形式:

為方便起見(jiàn)，該單元的第一個(gè)節(jié)點(diǎn)記為節(jié)點(diǎn)i.

由于傳統(tǒng)的最小二乘近似a=(PTP)-1PTui使得節(jié)點(diǎn)i的位移近似值不等于該點(diǎn)的位移值，即ui≠p(xi，yi)a，導(dǎo)致位移條件施加比較困難.故為了使節(jié)點(diǎn)i的位移近似值等于該點(diǎn)的位移值，利用方程(5)中的第一個(gè)方程解出a1，再?gòu)姆匠?5)其余的方程中消去a1，可得

則由最小二乘法得

節(jié)點(diǎn)i鄰域內(nèi)位移又可以表示為

則

并將式(7)、(9)帶入式(8)可得ui(x，y)表達(dá)式

ui(x，y)表達(dá)式(10)可簡(jiǎn)寫成:

此處，Φi為局部節(jié)點(diǎn)位移函數(shù)的形函數(shù)，可寫成:

其中1為所有元素均為1的(N-1)行列向量，利用式(4)進(jìn)而可以求出單元的形函數(shù)矩陣Ψ.

從式(10)可以看出，對(duì)于(xi，yi)點(diǎn)存在ui(xi，yi)=ui，再根據(jù)式(4)可得 Ψ = ［1，0，0|0，…，0］.同理對(duì)于 i=2，3 點(diǎn)，存在.故該有限元無(wú)網(wǎng)格耦合三角形單元的形函數(shù)具有Kronecker delta性質(zhì)，能夠像普通的有限元一樣直接施加位移邊界條件.

1.3 自由振動(dòng)和強(qiáng)迫振動(dòng)分析

由哈密頓原理，可得系統(tǒng)的運(yùn)動(dòng)方程為

為有限元無(wú)網(wǎng)格耦合三角形單元的剛度矩陣，位移應(yīng)變矩陣B的維數(shù)是3×2M，M為單元支持域Ω內(nèi)的節(jié)點(diǎn)數(shù).

M為質(zhì)量矩陣:

其中，Q是(4)式中元素組成的2×2M矩陣.

C為阻尼矩陣，為簡(jiǎn)便起見(jiàn)，取瑞利阻尼:

其中α，β是瑞利阻尼系數(shù).

荷載列陣F:

文中運(yùn)動(dòng)方程(13)的求解采用Newmark方法［7-8］，求解參數(shù)取 α =0.25，δ=0.5.

若阻尼和荷載項(xiàng)為零，則得自由振動(dòng)方程:

2 數(shù)值算例

2.1 錐形梁的自由振動(dòng)

如圖2(a)所示錐形懸臂梁，梁厚t=0.05 m，彈性模量E=200 GPa，泊松比 μ=0.3，ρ=8 000 kg/m3，為了與四邊形等參元比較，采用6×4形式劃分網(wǎng)格，如圖2(b)所示，四邊形采用畸變的網(wǎng)格.前六階頻率計(jì)算結(jié)果列于表1.

圖2 錐形梁的幾何模型和網(wǎng)絡(luò)Fig.2 The geometry and mesh of the taper beam

從表1可以看出，有限元無(wú)網(wǎng)格耦合三角形單元沒(méi)有出現(xiàn)虛假的零能模式，并且其結(jié)果的誤差遠(yuǎn)小于線性三角形單元(FEM-T3)和四邊形等參元(FEM-Q4).

表1 錐形梁的固有頻率Tab.1 The natural frequencies for the taper beam

2.2 懸臂梁的強(qiáng)迫振動(dòng)

如圖3(a)所示矩形懸臂梁，彈性模量E=1Pa，μ =0.3，ρ=1.0 kg/m3，瑞利阻尼系數(shù) α =0.005，β=0.272.采用10×4劃分網(wǎng)格，如圖3(b)所示.考慮兩種荷載工況.

(1)如圖4(a)所示，在A點(diǎn)受一集中簡(jiǎn)諧荷載F(t)=cos(ωt)，ω =0.3，計(jì)算時(shí)間取t=100 s，時(shí)間步長(zhǎng)取Δt=1 s.圖5繪出了考慮阻尼時(shí)，線性三角形單元(FEM-T3)、四邊形等參元(FEMQ4)、八節(jié)點(diǎn)等參元(FEM-Q8)、本文耦合單元以及理論解［9］的A點(diǎn)豎向位移的動(dòng)力響應(yīng).從圖上可以看出，本文耦合單元和八節(jié)點(diǎn)等參元的振幅都非常接近理論解，精度比四邊形等參元高，遠(yuǎn)優(yōu)于線性三角形單元.圖6為無(wú)阻尼時(shí)，4種單元和理論解［9］的動(dòng)力響應(yīng)，同樣可以看出，本文耦合單元的結(jié)果更接近理論解和八節(jié)點(diǎn)等參元，比四邊形等參元和線性三角形單元精確.

上例自由振動(dòng)和本例的結(jié)果說(shuō)明筆者所提的耦合單元能夠應(yīng)用到動(dòng)力響應(yīng)分析中，并且具有較高的計(jì)算精度.

(2)如圖4(b)所示，t=0 s在A點(diǎn)受一突加集中荷載F=1 N，持續(xù)時(shí)間t=2 000 s，時(shí)間步長(zhǎng)取Δt=1 s.利用筆者所提的耦合單元，計(jì)算A點(diǎn)豎向位移的動(dòng)力響應(yīng)，其結(jié)果與理論解［9］一并繪于圖7.

從圖7同樣可以看出，該耦合單元精度較高.在沒(méi)有阻尼時(shí)，A點(diǎn)在做振幅是接近常數(shù)的受迫振動(dòng).考慮阻尼時(shí)，由于阻尼的作用，位移的振動(dòng)響應(yīng)隨著時(shí)間的延長(zhǎng)而迅速衰減，直至消失.

3 結(jié)論

將有限元無(wú)網(wǎng)格耦合三角形單元應(yīng)用于分析二維固體的自由振動(dòng)和強(qiáng)迫振動(dòng)響應(yīng).通過(guò)典型的數(shù)值算例計(jì)算表明，該單元不存在虛假的零能模式，計(jì)算響應(yīng)的振幅接近八節(jié)點(diǎn)等參元，其精度優(yōu)于線性三角形單元和線性四邊形等參元.該單元還可以應(yīng)用到板、殼等結(jié)構(gòu)的動(dòng)力分析中.

［1］ WILSON EL，TAYLORR L，DOHERTY WP，et al.Incompatible Displacement Modes in Numerical and Computer Models in Structural Mechanics［M］.New York:Academic Press，1973:43-57.

［2］ 秦力一，許德剛，周愛(ài)民.基于切應(yīng)力條件的廣義協(xié)調(diào)等參元［J］.鄭州大學(xué)學(xué)報(bào):工學(xué)版，2005，26(2):92-94.

［3］ 龍馭球，李聚軒，龍志飛，等.四邊形單元面積坐標(biāo)理論［J］.工程力學(xué)，1997，14(3):1-11.

［4］ BELYTSCHKO，LU Y Y，GU L.Element free Galerkin methods［J］.International Journal for Numerical Methods in Engineering，1994，37(2):229–256.

［5］ BABUSˇI，MELENK JM.The partition of unity method［J］.International Journal for Numerical Methods in Engineering，1997，40(4):727–758.

［6］ RAJENDRAN S，ZHANG B R.A “FE-meshfree”QUAD4 element based on partition of unity［J］.Computer Methods in Applied Mechanics and Engineering，2007，197(1-4):128 –147.

［7］ 陳國(guó)榮.有限單元法原理及應(yīng)用［M］.北京:科學(xué)出版社，2009:310-319.

［8］ ABBASSIAN F，DAWSWELL D J，KNOWLES N C.Free Vibration Benchmarks［M］.Glasgow:National Agency for Finite Element Methods ＆ Standards，1987:72-77.

［9］ 克拉夫，彭津.結(jié)構(gòu)動(dòng)力學(xué)［M］.北京:高等教育出版社，2006:308-314.