賈 程

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

0 引言

傳統(tǒng)的等參有限單元模型滿足變分一致性要求，能夠較好的分析許多實際的問題，但是存在剛度矩陣過剛、對網(wǎng)格畸變敏感等問題。剛度矩陣過剛帶來了自鎖問題和計算精度下降問題，對于線性三角形單元尤其明顯。對網(wǎng)格畸變敏感要求劃分問題域的網(wǎng)格質(zhì)量較高，如四邊形單元要求其夾角盡量保持90°，這給前處理帶來巨大的工作量。為了改善上述缺陷，學(xué)者們提出了許多方法:hybrid-Trefftz單元法[1]、無網(wǎng)格法[2]等，然而，這些方法還有局限性或不夠完善的地方。

在實際應(yīng)用時，工程師比較偏愛三角形單元。因為它們能夠被自動的生成。本文使用非對稱有限單元公式的概念，將線性三角形單元形函數(shù)和有限元無網(wǎng)格耦合三角形單元(FE-LSPIM TRI3單元)形函數(shù)分別作檢驗函數(shù)和試函數(shù)，構(gòu)造出基于彼得羅夫伽遼金公式的US-FE-LSPIM TRI3三角形單元。

1 US-FE-LSPIM TRI3三角形單元公式

1.1 FE-LSPIM TRI3三角形單元的形函數(shù)

圖1 支持域的定義

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

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

其中，φi=[φi1φi2φi3…φiN]，i=1，2，3，ui=[u1u2u3L uN]T。

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

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

將方程(2)代入式(1)得:

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

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

令該單元的第一個節(jié)點記為節(jié)點i。由于傳統(tǒng)的最小二乘近似a=(PTP)－1PTui使得節(jié)點i的位移近似值不等于該點的位移值，即ui≠p(xi，yi)a，導(dǎo)致位移條件施加比較困難。故為了使節(jié)點i的位移近似值等于該點的位移值，利用方程(5)中的第一個方程解出a1，再從方程(5)其余的方程中消去a1，并由最小二乘法得:

其中，φi為局部節(jié)點位移函數(shù)的形函數(shù);1為所有元素均為1的(N－1)行列向量;和的定義和文獻(xiàn)[3]中相似，只是少了x2y和xy2項。利用式(4)進(jìn)而可以求出單元的形函數(shù)矩陣ψ。

根據(jù)式(6)可以得出，在(xi，yi)點存在 ui(xi，yi)=ui，再根據(jù)式(4)可得FE-LSPIM TRI3三角形單元的形函數(shù)ψ具有Kronecker delta性質(zhì)，能夠像普通的有限元一樣直接施加位移邊界條件。

1.2 US-FE-LSPIM TRI3三角形單元公式

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

其中，δu為虛位移域;δε為響應(yīng)的虛應(yīng)變域;σ為真實應(yīng)力。

由于ψ為線性三角形單元形函數(shù)和點插值無網(wǎng)格法形函數(shù)的復(fù)合函數(shù)，具有高階的完備性，ψ插值的應(yīng)力函數(shù)σ^比線性三角形單元的應(yīng)力函數(shù)更加適合作試函數(shù)。對于檢驗函數(shù)，只要其滿足邊界條件和域內(nèi)連續(xù)條件，故選擇線性三角形單元形函數(shù)插值的虛位移域δ做檢驗函數(shù)。

將它們代入式(7)，并整理得:

2 檢驗算例

2.1 受彎矩作用的懸臂梁

如圖2a)所示懸臂梁，在端部受常彎矩M=20C2作用，彎矩形式分布在右端，尺寸L=100，C=10，彈性模量 E=1.0 ×107，泊松比 μ =0.3，ρ=8 000 kg/m3，考慮平面應(yīng)力狀態(tài)。

圖2 端部受彎矩的懸臂梁

采用2×6×2，2×10×2，4×16×2，6×20×2，8×28×2五種形式的三角形網(wǎng)格劃分問題域，為了與四邊形等參元比較，本問題也采用相應(yīng)的2×6，2×10，4×16，6×20，8×28五種四邊形網(wǎng)格劃分。圖2b)為2×6×2網(wǎng)格劃分示例。A點的位移計算結(jié)果列于表1。

表1 A點的位移

從表1中可以看出，US-FE-LSPIM TRI3單元能夠很好的重構(gòu)常彎矩下的解，而傳統(tǒng)的線性三節(jié)點三角形單元(T3)和四節(jié)點等參四邊形單元(Q4)的結(jié)果不夠精確。

2.2 Cook懸臂梁

本例為變截面的懸臂梁，在端部受一個分布剪力F=1/16作用，如圖3a)所示。彈性模量E=1.0，泊松比μ=1/3。本問題采用2×2×2，4×4×2，8×8×2，16×16×2四種三角形網(wǎng)格劃分形式。圖3b)所示為4×4×2的形式。

圖3 Cook懸臂梁

不同網(wǎng)格下A點的最小主應(yīng)力列于表2。相比于T3和Q4單元，US-FE-LSPIM TRI3單元顯示了最佳的結(jié)果。

表2 A點的最小主應(yīng)力

3 結(jié)語

使用非對稱有限單元法的概念，構(gòu)造出基于彼得羅夫伽遼金公式的US-FE-LSPIM TRI3三角形單元。由于線性三角形單元形函數(shù)和FE-LSPIM TRI3單元形函數(shù)插值的位移函數(shù)能夠滿足域內(nèi)位移連續(xù)性和完備性的要求，選擇它們分別作為檢驗函數(shù)和試函數(shù)。數(shù)值算例顯示，該US-FE-LSPIM TRI3三角形單元精度較高，優(yōu)于傳統(tǒng)的線性三角形單元和四邊形等參元。

[1] FREITAS J A，BUSSAMRA F L S.Three-dimensional hybrid-Trefftz stress elements[J].International Journal for Numerical Methods in Engineering，2000，47(5):927-950.

[2] 張 雄，胡 煒，潘小飛，等.加權(quán)最小二乘無網(wǎng)格法[J].力學(xué)學(xué)報，2003，35(4):425-431.

[3] 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.

[4] Timoshenko S P，Goodier J N.Theory of Elasticity[M].3rd edn.New York:McGraw-Hill，1934.

[5] 庫 克，馬爾庫斯，普利沙.有限元分析的概念與應(yīng)用[M].關(guān)正西，強(qiáng)洪夫，譯.第4版.西安:西安交通大學(xué)出版社，2007.