無網(wǎng)格局部Petrov-Galerkin法分析板彎曲的剪切自鎖問題

夏 平,胡瑋軍,程玉蘭

(1.湖南工程學(xué)院機械工程學(xué)院,湘潭 411101;2.邵陽學(xué)院機械與能源工程系,邵陽 422000)

無網(wǎng)格局部Petrov-Galerkin法分析板彎曲的剪切自鎖問題

夏 平1,胡瑋軍2,程玉蘭1

(1.湖南工程學(xué)院機械工程學(xué)院,湘潭 411101;2.邵陽學(xué)院機械與能源工程系,邵陽 422000)

用徑向基函數(shù)構(gòu)造無網(wǎng)格局部Petrov-Galerkin方法的形函數(shù),插值函數(shù)具有Kronecker delta函數(shù)性質(zhì),因此可以很方便地施加本質(zhì)邊界條件.分析了板彎曲時剪切自鎖現(xiàn)象產(chǎn)生的原因,利用無網(wǎng)格局部Petrov-Galerkin方法對兩對邊固支另對邊簡支中厚板的彎曲進行了分析和計算.發(fā)現(xiàn)無網(wǎng)格方法相對于有限元法等傳統(tǒng)數(shù)值方法能更好地避免剪切鎖死現(xiàn)象,并提出了避免剪切自鎖現(xiàn)象發(fā)生的有效措施.

無網(wǎng)格局部Petrov-Galerkin方法;徑向基函數(shù);彎曲問題;剪切自鎖

0 引 言

無網(wǎng)格方法是一類在有限元法等傳統(tǒng)數(shù)值方法基礎(chǔ)上,針對其網(wǎng)格單元存在的問題而提出的一類新的數(shù)值方法.無網(wǎng)格方法目前已有數(shù)十種之多[1,2],其中無網(wǎng)格Garlerkin法[3]和無網(wǎng)格Petrov-Garlerkin法[4]是最常用的兩種無網(wǎng)格方法.

無網(wǎng)格方法是基于節(jié)點信息而不是基于單元或網(wǎng)格信息,因此有利于避免剪切自鎖問題.用有限元法求解薄板問題時,常常會在板很薄的情況下發(fā)生剪切鎖死現(xiàn)象,于是采用一系列的方法來避免剪切自鎖現(xiàn)象的發(fā)生,如縮減積分法和假設(shè)剪切應(yīng)變方法等[5].本文通過用無網(wǎng)格局部Petrov-Garlerkin法對中厚板的彎曲進行計算來對剪切自鎖問題進行分析.

1 中厚板的基本方程

在圖1所示笛卡爾坐標(biāo)系中,取 xy為板的中面,取z軸與xy面垂直向下,設(shè)板厚度為h,x和y方向的長度分別為a,b.板內(nèi)任一點在變形后沿 x,y,z軸方向的位移分量記為u,v和w,用三個廣義位移 θx(x,y),θy(x,y),w(x,y)來表示三個位移分量,存在關(guān)系

其中θx、θy分別是垂直中面的直法線在xz和yz平面內(nèi)的轉(zhuǎn)角,w(x,y)為中面撓度.軸向位移和角位移的正方向如圖1所示.

圖1 矩形板

對于各向同性板,廣義內(nèi)力與廣義應(yīng)變之間的關(guān)系形式上可以表示為

或

式中,σ為廣義內(nèi)力向量

ε為廣義應(yīng)變向量

其中,εf為彎曲應(yīng)變向量,εs為剪切應(yīng)變向量

在式(3)中,D為板材料的彈性矩陣

其中D為板的彎曲剛度,C為板的剪切剛度.

式中,E、G和μ分別為板材料的彈性模量、剪切模量和泊松比.

2 剪切自鎖產(chǎn)生的原因

中厚板的應(yīng)變能包括彎曲應(yīng)變能和剪切應(yīng)變能兩部分.由中厚板的虛功原理,即彎曲虛應(yīng)變能和剪切虛應(yīng)變能之和等于外載荷的虛功可得

當(dāng)板變得比較薄的時候,方程(8)左邊的第二項即與剪切應(yīng)變相關(guān)的剪切應(yīng)變能項可能變得比與第一項彎曲應(yīng)變相關(guān)的彎曲應(yīng)變能大.在數(shù)值模擬中,都使用同階多項式函數(shù)近似位移變量w、θx和θy,因此當(dāng)板變得很薄時,第二項會變得比第一項大很多,這樣就會導(dǎo)致出現(xiàn)錯誤的結(jié)果,這種現(xiàn)象稱為剪切自鎖.有限元法中一般使用線性函數(shù)來插值位移變量 w 、θx和θy,由于 ?w/?x與 ?w/?y和 θy與 的階次不協(xié)調(diào),剪切應(yīng)變能比彎曲應(yīng)變能大很多,這時就很容易導(dǎo)致剪切自鎖現(xiàn)象的發(fā)生.

3 避免剪切自鎖現(xiàn)象的有效措施

3.1 采用高次多項式基函數(shù)

在無網(wǎng)格局部Petrov-Garlerkin法中,采用徑向基函數(shù)耦合多項式基函數(shù)來構(gòu)造形函數(shù),通過提高多項式基函數(shù)的階數(shù),可以很容易構(gòu)造任意階次的形函數(shù),因此能有效的避免剪切自鎖現(xiàn)象的發(fā)生.

為檢驗高次多項式基函數(shù)能夠有效的避免剪切自鎖現(xiàn)象,在此考慮對邊固支另對邊簡支中厚方板,取邊長為a=1 m,彈性模量E=200 GPa,泊松比μ=0.3,受均布面載荷ε0=1MPa.利用MQ徑向基函數(shù)附加二次多項式進行插值近似計算.板上布置225(15×15)個規(guī)則離散節(jié)點,具體計算方法和步驟參考文獻[6,7].假設(shè)板厚跨比為0.001,其中點撓度的級數(shù)解[8]為wmax=0.00192q0a4/Q.三個位移變量w、θx和θy使用相同的插值函數(shù)進行插值,插值函數(shù)使用不同階次的多項式基函數(shù)進行構(gòu)造.圖2給出了沿x軸方向(兩固支邊中點連線上)在不同階數(shù)多項式基函數(shù)的情況下?lián)隙鹊淖兓€.從圖中可以看出使用一階多項式基函數(shù)時(m=3),即線性多項式時,發(fā)生明顯的剪切自鎖現(xiàn)象.當(dāng)使用二次多項式基函數(shù)時(m=6),基本上消除了剪切自鎖現(xiàn)象.但是,當(dāng)多項式基函數(shù)的項數(shù)取值繼續(xù)增大時,則需要更多的節(jié)點來插值位移變量,因此影響域的大小要明顯增大,這樣會增大計算量.

圖2 不同階數(shù)多項式基函數(shù)下?lián)隙鹊淖兓€

3.2 位移變量分別插值

另一種避免剪切自鎖現(xiàn)象的有效方法是對三個位移變量采用不同階次的形函數(shù)進行插值,即轉(zhuǎn)角采用比撓度低一階的插值函數(shù).由于w的連續(xù)性要求比θx和θy高一個階次,因此可以采用如下插值方式

該方法能夠有效地消除剪切自鎖現(xiàn)象的發(fā)生.對于上面的算例,表1列出了在不同階數(shù)多項式基函數(shù)和不同厚跨比時,位于板中點的最大撓度值.

表1 分別插值時板中點的最大撓度值w(×10-3D/q0a4)

從表1中可以看出對撓度變量和兩個轉(zhuǎn)角變量采用不同階次的形函數(shù)進行分別插值時,即使使用線性多項式基函數(shù)(m=3)也能夠有效的避免剪切自鎖現(xiàn)象的發(fā)生.

剪切自鎖現(xiàn)象能得到有效地避免,就能提高求解板殼問題的效率并增強求解方法的通用性.

4 結(jié) 論

用無網(wǎng)格局部Petrov-Garlerkin法分析中厚板彎曲時的剪切自鎖問題可以得到如下結(jié)論:

(1)在板的厚跨比大于0.01的情況下,即使使用線性多項式基函數(shù)(m=3)也不會出現(xiàn)剪切自鎖現(xiàn)象.但通常選擇二階多項式基函數(shù)(m=6)可以獲得更加精確的結(jié)果.

(2)使用高階多項式基函數(shù),如二階多項式基函數(shù)(m=6)或二階以上的多項式基函數(shù)時,即使板的厚跨比變得很小,也可以基本上消除剪切自鎖.

(3)如果撓度和轉(zhuǎn)角分別使用不同階次的形函數(shù)進行插值,即轉(zhuǎn)角采用比撓度低一階的插值函數(shù),則完全能夠避免剪切自鎖現(xiàn)象的發(fā)生.

[1]曹國金,姜泓道.無單元法研究和應(yīng)用現(xiàn)狀及動態(tài)[J].力學(xué)進展,2002,32(4):526-534.

[2]張 雄,宋康祖,陸明萬.無網(wǎng)格法研究進展及其應(yīng)用[J].計算力學(xué)學(xué)報,2003,20(6):730-742.

[3]Lu Y Y.,Belytschko T,Gu L.A New Implementation of the Element-free Galerkin Method[J].Comput.Methods in Appl.Mech.and Engrg.,1994,113:397-414.

[4]龍述堯.彈性力學(xué)問題的局部Petrov-Galerkin方法[J].力學(xué)學(xué)報,2001,33(4):508-518.

[5]王勖成,邵 敏.有限單元法基本原理和數(shù)值方法(第2版)[M].清華大學(xué)出版社,1997.

[6]夏 平,龍述堯,崔洪雪.用無網(wǎng)格局部徑向點插值法分析中厚板的彎曲問題[J].應(yīng)用力學(xué)學(xué)報,2009,26(2):383-388.

[7]Xia P,Long SY,Cui HX,et al.The Static and Free Vibration Analysis for the Nonhomogeneous Moderately Thick Plate by the Meshless Local Radial Point Interpolation Method[J],Engineering Analysis with Boundary Elements,2009,33(6):770-777.

[8]S.鐵摩辛柯,S.沃諾斯基.板殼理論[M].北京:科學(xué)出版社,1977.

Shear Locking Analysis of Plate Bending Problem by Using Meshless Local Petrov-Galerkin Method

XIA Ping1,HU Wei-jun2,CHENG Yu-lan1
(1.College of Mechanical Engineering,Hunan Institute of Engineering,Xiangtan 411101,China;2.Department of Mechanical and Energy Engineering,Shaoyang Institute,Shaoyang 422000,China)

The shape function of the meshless local Petrov-Galerkin method is constructed by using the radial basis functions and possesses Kronecker Delta function properties.Therefore,the essential boundary conditions can be easily imposed.Causation of shear locking occurring in plate bending is analyzed.Bending problems for plate with two sides simply supported,the other two sides clamped boundary conditions,are analyzed by the meshless local Petrov-Galerkin method in this paper.It is found that the shear locking is easier to avoid in the meshless method than in the finite element method,and the measures of avoiding the shear locking are presented.

meshless local Petrov-Galerkin method;radial basis function;bending problem;shear locking

TU457

A

1671-119X(2011)01-0025-03

2010-09-07

湖南省自然科學(xué)基金資助項目(10JJ3036);湖南省教育廳科研資助項目(09C880)

夏 平(1968-),男,工學(xué)博士,副教授,研究方向:計算固體力學(xué)及其應(yīng)用.