US－FE－LSPIM四邊形單元性能進(jìn)一步研究★

賈 程

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

0 引言

有限單元法已經(jīng)成為工程領(lǐng)域中一個重要的分析工具。但是，傳統(tǒng)的等參元在網(wǎng)格畸變時精度下降給人們帶來很大的麻煩。在計算大變形問題時，由于精度下降還需要重新劃分網(wǎng)格。為提高有限單元法的精度，尋找更先進(jìn)的方法，學(xué)者們做了大量的研究工作。提出了QM6單元、增強(qiáng)假設(shè)應(yīng)變法(EAS)、四邊形面積坐標(biāo)法和無網(wǎng)格法(EFG)[1]等。

最近，RAJENDRAN 等[2]提出的 FE－LSPIM 四邊形單元屬于雜交單元。該單元結(jié)合有限單元法和無網(wǎng)格法的優(yōu)點(diǎn)，具有高精度和Kronecker－delta性質(zhì)。但是，該單元施加單元整體邊界需滿足的位移不方便，需要用罰函數(shù)或拉格朗日乘子法。文獻(xiàn)[3]改進(jìn)了它在單元整體邊界位移施加中的不足，提出US－FE－LSPIM四邊形單元。本文進(jìn)一步研究US－FE－LSPIM四邊形單元在網(wǎng)格畸變、剪切自鎖和體積自鎖時的性能。

1 US-FE-LSPIM四邊形單元

1.1 FE－LSPIM四邊形單元的形函數(shù)

單元內(nèi)位移u寫成:

其中，N′=[N1N2N3N4];Ni為四邊形等參元形函數(shù)，i=1，2，3，4。

ui(x，y)為節(jié)點(diǎn)位移函數(shù)，i=1，2，3，4，在該節(jié)點(diǎn)處等于其位移值，即:ui(xi，yi)=ui，i=1，2，3，4，位移函數(shù) ui(x，y)由 i點(diǎn)的支持域內(nèi)的節(jié)點(diǎn)值運(yùn)用最小二乘點(diǎn)插值法(LSPIM)得到:

其中，Φi=[Φi1Φi2Φi3… Φi

N]，i=1，2，3，4。

Φ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ù)。對于該單元內(nèi)的其他節(jié)點(diǎn)，N可能不同。一個單元所有節(jié)點(diǎn)支持域的合集構(gòu)成了一個單元的支持域Ω。

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

設(shè)單元支持域Ω內(nèi)的節(jié)點(diǎn)數(shù)為M，則1≤N≤M，方程(3)就可以寫成單元支持域的形式:

其中，Φ為相應(yīng)于單元支持域的LSPIM單元形函數(shù)矩陣;u為單元支持域內(nèi)所有節(jié)點(diǎn)的x方向的位移向量。

圖1 支持域的定義

將式(4)代入式(1)得:

從式(5)中，得到該單元的形函數(shù)矩陣:

Φ中的元素Φi可以由i點(diǎn)的支持域內(nèi)的節(jié)點(diǎn)值運(yùn)用最小二乘點(diǎn)插值法(LSPIM)得到:

其中，ˉQ，ˉq定義見文獻(xiàn)[2];1為所有元素均為1的列向量。則利用式(6)進(jìn)而可以求出單元的形函數(shù)矩陣ψ。

1.2 US－FE－LSPIM 四邊形單元

線彈性體在平衡狀態(tài)下的虛功方程為:

其中，b為體力;t為面力;δu為虛位移場;δε為相應(yīng)的虛應(yīng)變場;σ為真實(shí)的應(yīng)力場。

為了重構(gòu)一個完整的二次位移場，形函數(shù)需滿足以下方程:

FE－LSPIM四邊形單元形函數(shù)ψ能夠滿足式(9)所有的方程[2]，而常規(guī)的四邊形等參元形函數(shù)只能滿足低階的條件，即p，q取0或1。

由于方程(8)中的虛位移可以是任意的，只要它滿足本質(zhì)邊界條件和域內(nèi)連續(xù)性條件。常規(guī)的四邊形等參元形函數(shù)插值的位移場δˉun由于能夠滿足單元間C0階連續(xù)和單元內(nèi)C1階連續(xù)的要求，故是個合適的選擇。

在傳統(tǒng)的等參元中，σ由等參元形函數(shù)插值的應(yīng)力ˉσ代替。但是由于不能滿足式(9)中的高階條件，導(dǎo)致等參元在網(wǎng)格畸變時精度下降。為了提高在單元二次位移域下的計算精度，選擇FE－LSPIM插值的應(yīng)力場作為試函數(shù)，即

將以上的虛位移和應(yīng)力函數(shù)代入式(8)，整理可得:

2 數(shù)值算例

2.1 劃分為五個畸變網(wǎng)格的懸臂梁

如圖2所示，懸臂梁被劃分成五個畸變的網(wǎng)格，彈性模量E=1 500，泊松比υ=0.25，厚度t=1。考慮兩種荷載工況:a.彎矩 M作用下的純彎曲;b.橫向剪力作用下的線性彎曲。A點(diǎn)的豎直位移列于表1。

圖2 劃分為五個畸變網(wǎng)格的懸臂梁

表1 懸臂梁A點(diǎn)的豎直位移

從表1可以看出US－FE－LSPIM四邊形單元在畸變的網(wǎng)格下具有較高的精度。

2.2 剪切自鎖問題

為研究剪切自鎖問題，考慮如圖3所示懸臂梁。彈性模量E=107，泊松比 υ =0.25，厚度 t=0.1。彎矩 M=h，長寬比 L/h 分別取3，30，300和3 000，使用2×2高斯積分。梁端的豎直位移的計算結(jié)果列于表2。

圖3 懸臂梁剪切自鎖問題

表2 懸臂梁端部位移

從表2中可以看出，在彎矩作用下，隨著長寬比的增大，四邊形等參元Q4的精度急劇下降，表現(xiàn)出剪切鎖定現(xiàn)象。QM6單元和US－FE－LSPIM單元在長寬比3～300時，精度變化不大;在長寬比為3 000時，QM6單元精度下降較大，而 US－FE－LSPIM單元精度略有下降，顯示了較強(qiáng)的抗剪切自鎖性能。

2.3 體積自鎖問題

如圖4所示的懸臂梁，受端部彎矩作用，考慮平面應(yīng)變問題。彈性模量 E=1 500，泊松比 υ 分別取 0.25，0.49，0.499，0.499 9。e分別取0和3。A點(diǎn)的豎直位移計算結(jié)果列于表3。

圖4 懸臂梁體積自鎖問題

表3 懸臂梁體積自鎖問題A點(diǎn)的豎直位移

從表3中可以看出，Q4四邊形單元在e等于0和3時都表現(xiàn)出明顯的體積自鎖現(xiàn)象。和QM6單元一樣，US－FE－LSPIM單元無體積自鎖現(xiàn)象。但在e=3時，US－FE－LSPIM單元的精度高于QM6單元。

3 結(jié)語

本文進(jìn)一步研究了US－FE－LSPIM單元的靜力性能。通過數(shù)值算例表明，該單元在畸變網(wǎng)格下具有較高的精度，優(yōu)于Q4，QM6單元，并且該單元顯示了較強(qiáng)的抗剪切自鎖和體積自鎖時的能力。

[1] BELYTSCHKO，LU Y Y，Gu L.Element free Galerkin methods[J].International Journal for Numerical Methods in Engineering，1994，37(2):229－256.

[2] 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):128－147.

[3] 賈 程，陳國榮，陳卉卉.US－FE－LSPIM 四邊形單元及其在幾何非線性問題中的應(yīng)用[J].計算力學(xué)學(xué)報，2011，28(5):785－791.

[4] Chen XM，Cen S，Long YQ，et al.Membrane elements insensitive to distortion using the quadrilateral area coordinate method[J].Computers and Structures，2004，82(1):35－54.

[5] 鐵摩辛柯，古地爾.彈性理論[M].第3版.徐芝綸，譯.北京:高等教育出版社，1990.