基于虛功原理的梁柱有限元分析

基于虛功原理的梁柱有限元分析

在建筑結(jié)構(gòu)分析中，梁柱是一個(gè)常見(jiàn)的結(jié)構(gòu)。在用計(jì)算機(jī)分析時(shí),梁柱單元的切線剛度矩陣就很重要了。目前,國(guó)內(nèi)采用較多的是CenapOran的梁柱理論。以虛功原理為理論基礎(chǔ)的有限元是一種有效的力學(xué)分析方法,遺憾的是國(guó)內(nèi)很少見(jiàn)到使用它對(duì)梁柱進(jìn)行非線性分析的文章。如果按照通常有限元分析步驟,應(yīng)用多項(xiàng)式插值函數(shù)推導(dǎo)得到的幾何剛度矩陣計(jì)算時(shí),在荷載較小時(shí)候就發(fā)散了。這也導(dǎo)致了在進(jìn)行梁柱分析中,國(guó)內(nèi)仍然普遍采用以梁柱理論為基礎(chǔ)的梁柱單元,而以虛功原理為基礎(chǔ)的梁柱單元到目前為止仍然采用較少。本文推導(dǎo)出了幾何剛度矩陣?？紤]了因節(jié)點(diǎn)位移導(dǎo)致節(jié)點(diǎn)力的不平衡,推出了節(jié)點(diǎn)力的修正矩陣KE。當(dāng)考慮了節(jié)點(diǎn)力的修正矩陣后,數(shù)值計(jì)算結(jié)果與前人的計(jì)算結(jié)果符合良好。1梁?jiǎn)卧臉?gòu)形標(biāo)記在梁柱的分析中,一般作如下假定:(1)平截面假定,即彎曲變形前為平面的截面,變形后仍然為平面;(2)截面內(nèi)的剪應(yīng)變形不計(jì);(3)等截面桿件。本文中,假定梁為雙對(duì)稱軸的截面,分別以兩主軸和梁的縱軸建立坐標(biāo)系,如圖1所示。根據(jù)假定,距梁左端距離為x截面上的任意一點(diǎn)的位移關(guān)系有uxr=ux-yuy,x-zuz,x,uyr=uy-zθx,uzr=uz+yθx(1)上式中,ux,uy,uz軸線OAOB的位移,θx為截面的轉(zhuǎn)角,uxr,uyr,uzr是截面上任一點(diǎn)的位移;則應(yīng)變?yōu)棣舩x=exx+ηxx,εxy=exy+ηxy,εxz=exz+ηxz(2)其中εij=12[?ui?xj+?uj?xi],ηij=12[?u1?xi?u1?xj+?u2?xi?u2?xj+?u3?xi?u3?xj]exx=ux,x-yuy,xx-zuz,xx(3)ηxx=12(U2x,x+U2y,x+U2z,x+y2U2y,xx+z2U2z,xx)+12[(y2+z2)θ2x,x+2(yuz,x-zuy,x)θx,x-2(yuy,xx+zuz,xx)ux,x+2yzuy,xxuz,xx]exy=12zθx,x,ηxy=12(-ux,xuy,x+yuy,xxuy,x+zuy,xuz,xx+uz,xθx+yθx,xθx)(4)exz=12yθx,x,ηxz=12(-ux,xuz,x+yuy,xxuz,x+zuz,xxuz,x-uy,xθx+zθx,xθx)(5)εij=12[?ui?xj+?uj?xi],ηij=12[?u1?xi?u1?xj+?u2?xi?u2?xj+?u3?xi?u3?xj]exx=ux,x?yuy,xx?zuz,xx(3)ηxx=12(U2x,x+U2y,x+U2z,x+y2U2y,xx+z2U2z,xx)+12[(y2+z2)θ2x,x+2(yuz,x?zuy,x)θx,x?2(yuy,xx+zuz,xx)ux,x+2yzuy,xxuz,xx]exy=12zθx,x,ηxy=12(?ux,xuy,x+yuy,xxuy,x+zuy,xuz,xx+uz,xθx+yθx,xθx)(4)exz=12yθx,x,ηxz=12(?ux,xuz,x+yuy,xxuz,x+zuz,xxuz,x?uy,xθx+zθx,xθx)(5)梁?jiǎn)卧跁r(shí)刻t=0的初始構(gòu)形記為C1,在當(dāng)前時(shí)刻t的構(gòu)形記為C1,把欲求的時(shí)刻t+Δt的構(gòu)形記為C2。在以下的表述中,左側(cè)的上下標(biāo),均表示構(gòu)形;其中左上標(biāo)表示所示的量在那個(gè)構(gòu)形中出現(xiàn),左下標(biāo)表示所示的量在那個(gè)構(gòu)形中測(cè)量,缺省左下標(biāo)表示所示量在C1中測(cè)量得到,缺省左上標(biāo)表示所示的量為C1到C2的增量。在上面的表述中右側(cè)上下標(biāo)表示微分。則以t時(shí)刻的構(gòu)形描述的虛功原理的形式為∫1VCijrs1ersδeij1dV+∫1V01111Sijδηij1dV+01111R=02121R(6)式中Sij為時(shí)刻t的第二類Piola-Kirchhoff應(yīng)力,01111R=∫1V01111SijδeijdV,Cijrs為本構(gòu)關(guān)系;對(duì)于各向同性材料,有Sxx=Eεxx,Sxy=2Gεxy,Sxz=2Gεxz(7)若ux,θx考慮線性插值,uy,θy,uz,θz采用Hermite插值,則ux=Ν1～Τux?θx=Ν1～Τθx,uy=Ν2～Τuy,uz=Ν3～Τuz(8)ux=N1～Tux?θx=N1～Tθx,uy=N2～Tuy,uz=N3～Tuz(8)其中～Τux=[uxAuxB]?～Τuy=[uyAθzAuyBθzB]?～Τθx=[θxAθxB]?～Τuz=[uzAθyAuzBθyB]?Ν1=[1-xlxl]?Ν2=[n1n2n3n4]?Ν3=[n1-n2n3-n4]n1=1-3[xl]2+2[xl]3,n2=lx-2l[xl]2+l[xl]3,n3=3[xl]2-2[xl]3,n4=-l[xl]2+l[xl]3～Tux=[uxAuxB]?～Tuy=[uyAθzAuyBθzB]?～Tθx=[θxAθxB]?～Tuz=[uzAθyAuzBθyB]?N1=[1?xlxl]?N2=[n1n2n3n4]?N3=[n1?n2n3?n4]n1=1?3[xl]2+2[xl]3,n2=lx?2l[xl]2+l[xl]3,n3=3[xl]2?2[xl]3,n4=?l[xl]2+l[xl]3將(8)代入(3)(4)(5)分別得exx、exy、exz,把(3)(4)(5)(7)代入(6)得線性剛度矩陣KL和幾何剛度矩陣KG(見(jiàn)附錄A2);許多文獻(xiàn)忽略了軸向應(yīng)變的非線性項(xiàng),即(3)中的ηxx,得到了幾何剛度矩陣,本文考慮了這一項(xiàng)的影響,因而理論上更為完備,計(jì)算精度更高。2f以2f1f為單位的不平衡力fnlGattass指出了用Updated-Lagrangian求解時(shí),C1中平衡的節(jié)點(diǎn)力,在C2中不再平衡,下面作一引述。用Updated-Lagrangian求解,(6)忽略了應(yīng)變的二次及三次項(xiàng),完整的表述如下021R=∫1VCijrsersδeij1dV+∫1V011Sijδηij1dV+∫1V1Sijδeij1dV+Q+C(9)其中,Q=∫1VCijkl(eklδηij+ηklδeij)1dV,C=∫1VCijklηklδηkl1dV(9)式準(zhǔn)確的表述了力在C2內(nèi)的平衡,當(dāng)C1、C2相隔很近時(shí),忽略Q、C即得(6),當(dāng)位移插值以后,(9)變成了節(jié)點(diǎn)力的平衡,2F=1F+(ΚL+ΚΝL)△+Q+C(10)或者改寫(xiě)為,2F=1F+FL+FΝL)+Q+C(11)(10)式中,KL為線性剛度矩陣,KNL即為幾何剛度矩陣KG,Δ是C1和C2之間的節(jié)點(diǎn)增量位移。注意(10)(11)描述的是節(jié)點(diǎn)在C2內(nèi)的平衡。(11)式忽略Q、C,若僅考慮Δ的二階分析,平衡仍然成立,即2F=1F+FL+FΝL(12)(12)式中,FL在C1中平衡,在C2中仍然平衡;而FNL從C1到C2不再平衡,因而FNL中包含了1F從C1到C2的不平衡力ΔF,因而必須從FNL中扣除ΔF方能是2F-1F的值,即增量荷載。以上的具體分析參見(jiàn)文獻(xiàn)。下面計(jì)算1F從C1到C2產(chǎn)生的不平衡力。如圖2、3所示,圖2為構(gòu)件的C1,圖3為構(gòu)件的C2。在C2中計(jì)算合外力,計(jì)算彎矩時(shí),對(duì)A端的坐標(biāo)軸取矩得,∑X=0,∑Y=0,∑Z=0(13)∑Μx=ΜzA+ΜzBluzB-ΜzA+ΜzBluzA+ΜyA+ΜyBluyA-ΜyA+ΜyBluyB≠0(14)∑Μy=ΜyA+ΜyBluxB-ΜyA+ΜyBluxA+FxBuzB-FxBuzA-ΜzAθxA-ΜzBθxB-12ΜxBθzA+12ΜxBθzB≠0(15)∑Μz=ΜzA+ΜzBluxA-ΜzA+ΜzBluxB-FxBuyB+FxBuyA+ΜyAθxA+ΜyBθxB+12ΜxBθyA-12ΜxBθyB≠0(16)(15)式中,由于θx產(chǎn)生的彎矩增量,參見(jiàn)文。分別在(14)(15)中令uxA,uxB,uyA,uyB,uzA,uzB,θxA,θxB不為零,得KE(見(jiàn)附錄A2)。3節(jié)點(diǎn)合力及幾何數(shù)據(jù)集的計(jì)算本文旨在說(shuō)明幾何剛度矩陣的推導(dǎo),和考慮不平衡力修正剛度矩陣的得來(lái),考慮到荷載增量法簡(jiǎn)單易行,因而本文使用荷載增量法求解。i讀入數(shù)據(jù),并初始化:單元幾何剛度矩陣KL+KNL=0,不平衡力調(diào)整矩陣KE=0,單元節(jié)點(diǎn)力F=0,單元節(jié)點(diǎn)位移δ=0;ii增加荷載增量因子;iii讀出單元的結(jié)點(diǎn)位移δ、KL+KNL,計(jì)算結(jié)點(diǎn)力F=(KL+KNL+KE)δ,求出節(jié)點(diǎn)合力R;iv形成KL;利用F形成KNL和KE,并存儲(chǔ)KL+KNL和KE;v用KE,KNL集成總剛K;vi讀入荷載因子,得到第i步的外荷載F,得增量荷載為ΔF=F-R;vii計(jì)算增量位移δ,求出單元的結(jié)點(diǎn)位移,更新節(jié)點(diǎn)幾何數(shù)據(jù);viii檢查增量位移的范數(shù)達(dá)到收斂標(biāo)準(zhǔn)否,若否,轉(zhuǎn)iii;ix檢查增量因子達(dá)到預(yù)定的值否,若否轉(zhuǎn)ii;x結(jié)果輸出。4梁柱單元的面特性如圖4,所示的剛架Williams做過(guò)理論分析并和試驗(yàn)作了比較,材料和截面特性為:E=7.1×104N/mm2,A=1.180cm2,I=0.0375cm4,采用