楊璞++牛紅攀++肖世富

摘要： 結(jié)合廣義有限元和理性有限元的優(yōu)勢，針對平面應力問題提出一種新型廣義四邊形單元.該單元考慮泊松效應，以節(jié)點位移自由度約束彈性力學平面應力方程的半解析解，構(gòu)造單元位移模式的附加項，較準確地反映真實位移場，提高單元的計算精度.推導新型廣義單元及其等參單元的形函數(shù)公式，設(shè)計分片試驗和數(shù)值算例驗證單元的精度.數(shù)值算例結(jié)果表明：在規(guī)則網(wǎng)格和非規(guī)則網(wǎng)格下新單元的計算精度均優(yōu)于傳統(tǒng)有限元和廣義有限元.新單元具有精度高且易于程序?qū)崿F(xiàn)的特點，可推廣應用到實際工程的結(jié)構(gòu)分析中.

關(guān)鍵詞： 平面四邊形單元； 等參元； 平面應力； 廣義有限元； 形函數(shù)； 分片試驗

中圖分類號： O242.21文獻標志碼： A

0引言

有限元法的基本思想是在力學模型上將一個原來連續(xù)的物體離散為有限個單元，對每個單元選擇一種簡單的函數(shù)來表示單元內(nèi)位移的分布規(guī)律，并按能量原理（變分原理）建立單元節(jié)點力和節(jié)點位移之間的關(guān)系，最后將所有單元集成并求解節(jié)點的位移.[1]

傳統(tǒng)線性有限元直接以數(shù)學為基礎(chǔ)采用雙線性多項式構(gòu)造形函數(shù)，其位移模式既沒有考慮彈性力學控制方程，也不能反映不同自由度間位移的相互影響，單元精度較低.為改善計算精度，20世紀90年代鐘萬勰院士團隊研發(fā)理性有限元，核心思想是從彈性力學出發(fā)建立位移模式，其插值函數(shù)的各項均來源于彈性力學精確解，有著明確的物理意義[27]，但理性有限元形函數(shù)繁復，給實際操作帶來很大困難[89].近年來，張洪武教授團隊在理性有限元的基礎(chǔ)上發(fā)展廣義有限元方法，核心思想是改變傳統(tǒng)有限元只利用單元節(jié)點同方向自由度的位移構(gòu)造方式，而以單元所有節(jié)點自由度來構(gòu)造單元位移[1011]，其位移模式沿用傳統(tǒng)有限元的多項式插值函數(shù)，未考慮彈性力學控制方程.

本文針對平面應力問題，提出一種新型廣義四邊形單元，綜合考慮理性有限元和廣義有限元的優(yōu)勢，從彈性力學平面應力方程出發(fā)建立形函數(shù)，并在不改變單元自由度的情況下為形函數(shù)引入位移附加項，提高單元計算精度.該新單元只在傳統(tǒng)有限元形函數(shù)基礎(chǔ)上增加位移附加項，程序操作簡單易于實現(xiàn).

1平面應力新型廣義四邊形單元

1.1一般單元

傳統(tǒng)平面矩形單元為4個節(jié)點八自由度，見圖1.

圖 1標準的平面四邊形單元

Fig.1Standard plane quadrilateral element

用4個節(jié)點位移表達的位移模式為

u=N1uI+N2uJ+N3uK+N4uL

v=N1vI+N2vJ+N3vK+N4vL

（1）其中，插值函數(shù)為雙線性函數(shù)，N1=（1-ξ）（1-η）/4

N2=（1+ξ）（1-η）/4

N3=（1+ξ）（1+η）/4

N4=（1-ξ）（1+η）/4 （2）該形函數(shù)的構(gòu)造局限于數(shù)學層面，位移模式只包含節(jié)點同方向自由度，計算精度較低.

由彈性力學理論的泊松效應可知，一個方向的變形將引起另一個方向的位移.設(shè)計開發(fā)新型平面四邊形單元，考慮彈性體的泊松效應，由單元的所有自由度構(gòu)造單元位移模式，即單元內(nèi)任意一點的位移分量不僅與4個節(jié)點位移分量方向的節(jié)點自由度相關(guān)，而且與4個節(jié)點處位移分量方向正交的節(jié)點自由度相關(guān).該開發(fā)思想與廣義有限元[1011]一致，區(qū)別是在建立形函數(shù)時進一步借鑒理性有限元思想，在傳統(tǒng)有限元位移模式的基礎(chǔ)上，根據(jù)節(jié)點位移自由度約束彈性力學平面應力方程，確定單元位移模式的附加項.新平面四邊形單元的位移模式既考慮節(jié)點上所有的自由度，又考慮彈性力學控制方程，能較準確反映真實位移場，因此可能有更高的計算精度.以下針對平面應力情況進行形函數(shù)構(gòu)造.

對于式（1），節(jié)點自由度確定的單元節(jié)點邊界條件為u（-1，-1）=uI，u（1，-1）=uJ，

u（1，1）=uK，u（-1，1）=uL

v（-1，-1）=vI，v（1，-1）=vJ，

v（1，1）=vK，u（-1，1）=vL （3）為考察一個自由度對另一自由度位移的影響，以第1個節(jié)點的x方向自由度u1為例，分析其他節(jié)點自由度固定、僅在節(jié)點I施加x方向位移uI時對y方向位移的貢獻.滿足圖1節(jié)點位移的條件表達式為u1（-1，-1）=uI

u1（1，-1）=u1（1，1）=u1（-1，1）=0

v1（-1，-1）=v1（1，-1）=v1（1，1）=

v1（-1，1）=0 （4）線彈性平面應力問題的齊次位移平衡微分方程為

22ux2+（1-μ）2uy2+（1+μ）2vxy=0

（1-μ）2vx2+22vy2+（1+μ）2uxy=0

（5）取滿足式（4）第一行方程的形函數(shù)基為u1（x，y）=uI（1-x）（1-y）/4 （6）將式（6）代入式（5）的第1個方程得2v1xy=0 （7）將式（6）代入式（5）的第2個方程且考慮節(jié)點邊界條件得（1-μ）2v1x2+22v1y2+（1+μ）14uI=0

v1（-1，-1）=v1（1，-1）=v1（1，1）=

v1（-1，1）=0 （8）結(jié)合式（7）和式（8），假設(shè)第1個自由度uI產(chǎn)生的y方向的位移為v1（x，y）=b0+b1x+b2y+b3x2+b4y2 （9）代入式（8）得b1=b2=0

b3=-2b0/（μ+1）+uI/8

b4=-（μ-1）b0/（μ+1）-uI/8 （10）由于位移v1完全由節(jié)點位移uI引起，設(shè)b0=βuI，則 v1（x，y）=β+-2μ+1β+18x2+

-μ-1μ+1β-18y2uI （11）代入坐標0，0得v1（0，0）=βuI （12）式（12）表明：位移v1中的常數(shù)項是由于泊松效應，由單元節(jié)點自由度位移分量引起的單元中心點該位移分量正交方向上的位移.通過ANSYS有限元計算并采用響應面方法，擬合得到的β與泊松比μ的關(guān)系式為β=0.152 1+0.023 78μ （13）式（13）的數(shù)據(jù)擬合圖見圖2.endprint

圖 2平面應力單元因數(shù)β與泊松比μ之間的關(guān)系

Fig.2Relationship between plane stress element factor β and Poisson ratio μ

將式（13）代入式（11）得v1（x，y）=（c1（1-x2）+c2（1-y2））uI （14）其中c1=2μ+1（0.152 1+0.023 78μ）-18

c2=μ-1μ+1（0.152 1+0.023 78μ）+18 （15）同理，其他自由度的位移可表示為v2（x，y）=-（c1（1-x2）+c2（1-y2））uJ

v3（x，y）=（c1（1-x2）+c2（1-y2））uK

v4（x，y）=-（c1（1-x2）+c2（1-y2））uL （16）施加x方向節(jié)點位移約束將產(chǎn)生y方向的位移，可以看作x方向位移對y方向位移的貢獻.綜合考慮4個節(jié)點的位移附加項可得到位移v的表達式為v=N1vI+N2vJ+N3vK+N4vL+

N5（uI+uK）+N6（uJ+uL） （17）其中N5=c1（1-x2）+c2（1-y2）

N6=-（c1（1-x2）+c2（1-y2）） （18）同理可得到水平位移u的表達式為u=N1uI+N2uJ+N3uK+N4uL+

N7（vI+vK）+N8（vJ+vL） （19）其中：N7=c2（1-x2）+c1（1-y2）

N8=-（c2（1-x2）+c1（1-y2）） （20）寫成矩陣表達式為u

v=N1N7N2N8N3N7N4N8

N5N1N6N2N5N3N6N4·

[uIvIuJvJuKvKuLvL]T（21）可以證明：在各節(jié)點處，N1=N2=N3=N4=1，N5=N6=N7=N8=0，滿足插值函數(shù)的Kronecker delta性質(zhì)；單元中任一點N1+N2+N3+N4=1，N5+N6=0，N7+N8=0，滿足插值函數(shù)的歸一性質(zhì).

1.2等參單元

傳統(tǒng)等參元中整體坐標與局部坐標的映射關(guān)系為x′=N1x′I+N2x′J+N3x′K+N4x′L

y′=N1y′I+N2y′J+N3y′K+N4y′L （22）式中：N1，N2，N3和N4為式（2）在局部坐標系下的表達式.根據(jù)第1.1節(jié)中的思路，平面應力問題新型廣義等參元位移模式為

u′=N1u′I+N2u′J+N3u′K+N4uL+N9v′I+

N10v′J+N11v′K+N12v′L

v′=N1v′I+N2v′J+N3v′K+N4v′L+N5u′I+

N6u′J+N7u′K+N8u′L （23）

其中各附加等參形函數(shù)為N5=β1（c1（1-ξ2）+c2（1-η2））

N6=-β2（c1（1-ξ2）+c2（1-η2））

N7=β3（c1（1-ξ2）+c2（1-η2））

N8=-β4（c1（1-ξ2）+c2（1-η2））

N9=β1（c2（1-ξ2）+c1（1-η2））

N10=-β2（c2（1-ξ2）+c1（1-η2））

N11=β3（c2（1-ξ2）+c1（1-η2））

N12=-β4（c2（1-ξ2）+c1（1-η2）） （24）式中：βi為面積因數(shù)，與單元節(jié)點所圍三角形面積相關(guān).圖3b中，定義SΔJKL，SΔKLI，SΔLIJ和SΔIJK分別為三角形JKL，KLI，LJK和IJK的面積，S為四邊形總面積，則βi的具體表達式為β1=2SΔJKL/S

β2=2SΔKLI/S

β3=2SΔLIJ/S

β4=2SΔIJK/S （25）在標準單元中，β1=β2=β3=β4=1.至此，平面應力問題新型廣義等參單元的位移函數(shù)已給出，求解過程與傳統(tǒng)的等參元相同.

a）b）圖 3新型4節(jié)點廣義等參元

Fig.3New generalized 4nodes isoparametric element

2分片試驗

對平面應力問題新型廣義等參四邊形單元的收斂性進行分片測試.分片單元見圖4.剛體位移測試：假定剛體位移模式為u=1，v=1，即在所有外節(jié)點處u1=u2=u3=u4=1

v1=v2=v3=v4=1 （26）計算內(nèi)部節(jié)點5～8的位移.計算結(jié)果見表1.

常應變測試：對邊界上的節(jié)點給以常應變位移場 u=1+3x+y

v=2+x+3y （27）采用1個Gauss積分點方案，內(nèi)部節(jié)點位移結(jié)果見表2.結(jié)果表明：新單元能通過剛體位移測試，在1個積分點條件下通過常應變測試.

圖 4等參單元分片試驗

Fig.4Isoparametric element patch test

表 1剛體位移分片測試結(jié)果

Tab.1Rigid body displacement patch test results節(jié)點坐標理論位移計算結(jié)果xyuvuv50.20.2111.001.0060.80.4111.001.0070.80.6111.001.0080.20.8111.001.00

表 2常應變分片測試結(jié)果

Tab.2Constant strain patch test results節(jié)點坐標理論位移計算結(jié)果xyuvuv50.20.21.82.81.802.8060.80.43.84.03.804.0070.80.64.04.64.004.6080.20.82.44.62.404.60

3單元基本變形模態(tài)endprint

為考察平面應力新型廣義等參單元的基本變形模態(tài)，采用圖5的矩形進行等參單元特征值分析.

圖 5矩形等參單元幾何尺寸

Fig.5Geometrical size of rectangle isoparametric element

采用4點Gauss積分方案計算單元剛度矩陣特征值，并與廣義等參元和傳統(tǒng)等參元的特征值進行比較.材料參數(shù)E=1，μ=0.3，結(jié)果見表3.

表 3單元特征值比較

Tab.3Eigenvalue comparison of elements模態(tài)傳統(tǒng)等參元廣義等參元新等參元1階0002階0003階0004階 0.44 0.35 0.255階 0.63 0.51 0.526階 0.63 0.63 0.637階 0.83 0.83 0.838階 1.75 1.75 1.75

根據(jù)特征值檢驗法[12]，不同的合法有效單元的8個特征值應含有3個零特征值（分別對應3個剛體運動）和3個彼此相等的非零特征值（分別對應3個常應變模態(tài)），另外2個不同的特征值分別關(guān)于x軸和y軸彎曲.表3中3種單元的1～3階特征值均為0，6～8階特征值分別相等，僅第4和5階特征值出現(xiàn)不同，說明廣義等參單元和新等參單元均合法有效.2種單元較傳統(tǒng)單元特征值均有所下降，表明2種單元均有更好的柔性，且新等參元的4階特征值較廣義等參元下降更快，表明新單元在關(guān)于x軸彎曲上比廣義等參元具有更好的柔性.第4和5階單元的變形模態(tài)見圖6.

a）4階b）5階圖 6矩形等參單元的變形模態(tài)

Fig.6Deformation modes of rectangle isoparametric element

采用9點Gauss積分方案，新等參元所得各階特征值與4點積分所得值相同，表明新等參元在4積分點下已經(jīng)具有較高的計算精度.

4數(shù)值算例

4.1算例1

算例模型見圖7a，長L=10 cm，高h=2 cm，材料常數(shù)E=2.1×105 MPa，泊松比μ=0.3.結(jié)構(gòu)左端所有節(jié)點2個方向自由度均固定，右端點A集中加載P=1 000 N的剪力作用.對新型廣義有限元、廣義有限元和傳統(tǒng)有限元3種單元用不同密度的網(wǎng)格進行計算結(jié)果比較.在剪力作用下該模型右端點A的平面應力豎直位移表達式[1]為v=P6EI3Lx2-x3+Ph28IGx （26）可得A點的豎直位移理論值為v=2.474×10-6 m.

a）b）圖 7端點受剪力作用的結(jié)構(gòu)

Fig.7Structure of which end point is subjected to shear force

在平面應力情況下不同網(wǎng)格密度計算所得A點的豎直位移對比見表4，隨網(wǎng)格加密的收斂曲線見圖8，其中傳統(tǒng)有限元值經(jīng)ANSYS算出.

表 4平面應力下端點A豎直位移值

Tab.4Values of vertical displacement of end point A in plane stress condition 10-6 m網(wǎng)格

密度傳統(tǒng)

有限元廣義

有限元新型廣義

有限元理論值10×2-2.169-2.340-2.424-2.47420×4-2.369-2.418-2.440-2.47440×8-2.429-2.442-2.448-2.47480×16-2.447-2.450-2.452-2.474

圖 8平面應力下端點A豎直位移曲線

Fig.8Vertical displacement curves of point A in plane stress condition

由圖8可知：在同一網(wǎng)格密度下，平面應力新型廣義矩形單元比廣義有限元和傳統(tǒng)矩形單元計算精度明顯提高.

方向不變性測試：該模型旋轉(zhuǎn)90°后（見圖7b），所得結(jié)果與原結(jié)果相同，說明新單元具有方向不變性.

4.2算例2

MacNeal細長懸臂梁左端所有節(jié)點2個方向自由度均固定，右端點A加載向下集中載荷P=1 N，彈性模量E=1×106 Pa，泊松比μ=0.3.模型幾何尺寸見圖9.

圖 9MacNeal細長懸臂梁受集中載荷模型，m

Fig.9Model of MacNeal slender cantilever beam under

concentrated load， m

MacNeal細長梁在不同網(wǎng)格密度下的端點位移見表5和圖10.結(jié)果表明新單元能在一定程度上改善剪切自鎖問題.

表 5MacNeal細長懸臂端點A豎直位移值

Tab.5Values of vertical displacement of end point A of MacNeal slender cantilever beam m網(wǎng)格

密度傳統(tǒng)

有限元廣義

有限元新型廣義

有限元理論值6×1-0.010-0.009-0.012-0.10812×2-0.031-0.028-0.036-0.10824×4-0.067-0.063-0.072-0.10848×8-0.094-0.092-0.096-0.10896×16-0.104-0.103-0.105-0.108

圖 10MacNeal細長懸臂端點A豎直位移曲線

Fig.10Vertical displacement curves of end point A of MacNeal slender cantilever beamendprint

4.3算例3

為測試新等參單元對不規(guī)則四邊形網(wǎng)格的計算精度，本算例仍使用算例1的模型，采用平面應力單元對不規(guī)則網(wǎng)格劃分模型進行計算.網(wǎng)格劃分情況見圖11.a）不規(guī)則網(wǎng)格8×2

b）不規(guī)則網(wǎng)格16×4

c）不規(guī)則網(wǎng)格32×8

d）不規(guī)則網(wǎng)格64×16

圖 11不規(guī)則網(wǎng)格劃分模型

Fig.11Model with Irregular mesh

平面應力下不同網(wǎng)格密度模型計算所得端點豎直位移值對比見表6和圖12.計算結(jié)果表明新單元計算精度高于傳統(tǒng)有限元和廣義有限元.

表 6平面應力下不規(guī)則網(wǎng)格劃分端點A豎直位移值

Tab.6Vertical displacement values of end point A with irregular mesh in plane stress condition 10-6 m網(wǎng)格

密度傳統(tǒng)

等參元廣義

等參元新型廣義

等參元理論值8×2-1.374-1.753-1.982-2.47416×4-2.002-2.219-2.323-2.47432×8-2.319-2.390-2.419-2.47464×16-2.418-2.437-2.444-2.474

圖 12平面應力下不規(guī)則網(wǎng)格劃分端點A豎直位移曲線

Fig.12Vertical displacement curves of end point A with irregular mesh in plane stress condition

4.4算例4

Cook變截面懸臂梁左端所有節(jié)點2個方向自由度均固定，右端中點A加載向上的集中載荷P=1 000 N，材料彈性模量E=2.1×105 MPa，μ=0.3，模型尺寸見圖13.

圖 13Cook變截面懸臂梁端點受集中載荷模型，cm

Fig.13Model of Cook variable cross section cantilever beam of which end point is under concentrated load， cm

Cook梁在不同網(wǎng)格密度下位移值對比見表7和圖14.結(jié)果表明新型廣義平面應力等參四邊形單元能提高計算精度.

表 7Cook變截面懸臂梁受集中載荷端點A豎直位移值

Tab.7Vertical displacement values of Cook variable cross section cantilever beam end point A under concentrated load 10-7 m網(wǎng)格

密度傳統(tǒng)

等參元廣義

等參元新型廣義

等參元2×20.6530.7000.7694×40.8890.9150.9448×81.0001.0101.01616×161.0481.0511.051圖 14Cook變截面懸臂梁受集中載荷端點A位移曲線

Fig.14Vertical displacement curves of Cook variable cross section cantilever beam end point A under concentrated load

5結(jié)束語

從彈性力學平面應力方程出發(fā)推導形函數(shù)，同時在插值函數(shù)中引入位移附加項，建立平面應力問題一種新的四邊形單元.進一步在位移附加項中引入面積因數(shù)，建立平面應力新型廣義等參四邊形單元.新等參單元能收斂于真實解，數(shù)值算例表明其計算精度高于傳統(tǒng)四邊形單元和廣義四邊形單元，且收斂速度更快.

本文方法只是在傳統(tǒng)的形函數(shù)上增加位移附加項，沒有增加單元自由度，程序修改簡單，實現(xiàn)方便，易于推廣.今后將在三維新型廣義單元等方面進一步開展研究工作.

參考文獻：

[1]王勖成. 有限元法[M]. 北京： 清華大學出版社， 2003.

[2]鐘萬勰， 紀崢. 理性有限元[J]. 計算結(jié)構(gòu)力學及其應用， 1996， 13（1）： 18.

ZHONG Wanxie， JI Zheng. Rational finite element RCQ8[J]. Comput Struct Mech & Appl， 1996， 13： 18.

[3]鐘萬勰， 紀崢. 平面理性元的收斂性證明[J]. 力學學報， 1997， 29（6）： 676685.

ZHONG Wanxie， JI Zheng. Convergence proof of plane rational finite element[J]. Acta Mech Sinica， 1997， 13（6）： 676685.

[4]紀崢， 鐘萬勰. 平面理性四節(jié)點及五節(jié)點四邊形有限元[J]. 計算力學學報， 1997， 14（1）： 1927.

JI Zheng， ZHONG Wanxie. Rational plane quadrilateral four and five nodes finite elements [J]. Chinese J Comput Mech， 1997， 14（1）： 1927.

[5]王永富， 鐘萬勰. 平面八節(jié)點四邊形理性元[J]. 固體力學學報， 2002， 23（1）： 109113.

WANG Yongfu， ZHONG Wanxie. Plane 8 nodes rational finite element[J]. Acta Mech Solida Sinica， 2002， 23（1）： 109113.

[6]紀崢， 劉澤佳， 趙偉. 平面理性八節(jié)點曲邊四邊形有限元RCQ8[J]. 大連理工大學學報， 2000， 40（1）： 4044.

JI Zheng， LIU Zejia， ZHAO Wei. 8Node curve quadrilateral rational finite element [J]. J Dalian Univ Technol， 2000， 40（1）： 4044.

[7]王永富， 鐘萬勰. 空間理性八節(jié)點塊體元[J]. 應用力學學報， 2003， 20（3）： 131135.

WANF Yongfu， ZHONG Wanxie. A rational hexahedron 8node finite element[J]. Chin J Appl Mech， 2003， 20（3）： 131135.

[8]ZHANG H W， ZHENG Y G， WU J K， et al. Generalized fournode plane rectangular and quadrilateral elements and their applications in the multiscale analysis of heterogeneous structures[J]. Int J Multiscale Comput Eng， 2013， 11（1）： 7191.（下轉(zhuǎn)第71頁）第24卷 第6期2015年12月計 算 機 輔 助 工 程Computer Aided EngineeringVol.24 No.6Dec. 2015endprint