平面應力問題新型廣義有限元法

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

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

關鍵詞： 平面四邊形單元； 等參元； 平面應力； 廣義有限元； 形函數； 分片試驗

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

0引言

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

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

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

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

1.1一般單元

傳統平面矩形單元為4個節點八自由度，見圖1.

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

Fig.1Standard plane quadrilateral element

用4個節點位移表達的位移模式為

u=N1uI+N2uJ+N3uK+N4uL

v=N1vI+N2vJ+N3vK+N4vL

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

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

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

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

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

對于式（1），節點自由度確定的單元節點邊界條件為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個節點的x方向自由度u1為例，分析其他節點自由度固定、僅在節點I施加x方向位移uI時對y方向位移的貢獻.滿足圖1節點位移的條件表達式為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）第一行方程的形函數基為u1（x，y）=uI（1-x）（1-y）/4 （6）將式（6）代入式（5）的第1個方程得2v1xy=0 （7）將式（6）代入式（5）的第2個方程且考慮節點邊界條件得（1-μ）2v1x2+22v1y2+（1+μ）14uI=0

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

v1（-1，1）=0 （8）結合式（7）和式（8），假設第1個自由度uI產生的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完全由節點位移uI引起，設b0=βuI，則 v1（x，y）=β+-2μ+1β+18x2+

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

圖 2平面應力單元因數β與泊松比μ之間的關系

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方向節點位移約束將產生y方向的位移，可以看作x方向位移對y方向位移的貢獻.綜合考慮4個節點的位移附加項可得到位移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）可以證明：在各節點處，N1=N2=N3=N4=1，N5=N6=N7=N8=0，滿足插值函數的Kronecker delta性質；單元中任一點N1+N2+N3+N4=1，N5+N6=0，N7+N8=0，滿足插值函數的歸一性質.

1.2等參單元

傳統等參元中整體坐標與局部坐標的映射關系為x′=N1x′I+N2x′J+N3x′K+N4x′L

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

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）

其中各附加等參形函數為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為面積因數，與單元節點所圍三角形面積相關.圖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.至此，平面應力問題新型廣義等參單元的位移函數已給出，求解過程與傳統的等參元相同.

a）b）圖 3新型4節點廣義等參元

Fig.3New generalized 4nodes isoparametric element

2分片試驗

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

v1=v2=v3=v4=1 （26）計算內部節點5～8的位移.計算結果見表1.

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

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

圖 4等參單元分片試驗

Fig.4Isoparametric element patch test

表 1剛體位移分片測試結果

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

表 2常應變分片測試結果

Tab.2Constant strain patch test results節點坐標理論位移計算結果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單元基本變形模態endprint

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

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

Fig.5Geometrical size of rectangle isoparametric element

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

表 3單元特征值比較

Tab.3Eigenvalue comparison of elements模態傳統等參元廣義等參元新等參元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

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

a）4階b）5階圖 6矩形等參單元的變形模態

Fig.6Deformation modes of rectangle isoparametric element

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

4數值算例

4.1算例1

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

a）b）圖 7端點受剪力作用的結構

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

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

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

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

密度傳統

有限元廣義

有限元新型廣義

有限元理論值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可知：在同一網格密度下，平面應力新型廣義矩形單元比廣義有限元和傳統矩形單元計算精度明顯提高.

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

4.2算例2

MacNeal細長懸臂梁左端所有節點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細長梁在不同網格密度下的端點位移見表5和圖10.結果表明新單元能在一定程度上改善剪切自鎖問題.

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

Tab.5Values of vertical displacement of end point A of MacNeal slender cantilever beam m網格

密度傳統

有限元廣義

有限元新型廣義

有限元理論值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

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

b）不規則網格16×4

c）不規則網格32×8

d）不規則網格64×16

圖 11不規則網格劃分模型

Fig.11Model with Irregular mesh

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

表 6平面應力下不規則網格劃分端點A豎直位移值

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

密度傳統

等參元廣義

等參元新型廣義

等參元理論值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平面應力下不規則網格劃分端點A豎直位移曲線

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

4.4算例4

Cook變截面懸臂梁左端所有節點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梁在不同網格密度下位移值對比見表7和圖14.結果表明新型廣義平面應力等參四邊形單元能提高計算精度.

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

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

密度傳統

等參元廣義

等參元新型廣義

等參元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結束語

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

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

參考文獻：

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

[2]鐘萬勰， 紀崢. 理性有限元[J]. 計算結構力學及其應用， 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]紀崢， 鐘萬勰. 平面理性四節點及五節點四邊形有限元[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]王永富， 鐘萬勰. 平面八節點四邊形理性元[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]紀崢， 劉澤佳， 趙偉. 平面理性八節點曲邊四邊形有限元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]王永富， 鐘萬勰. 空間理性八節點塊體元[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.（下轉第71頁）第24卷 第6期2015年12月計 算 機 輔 助 工 程Computer Aided EngineeringVol.24 No.6Dec. 2015endprint