索結(jié)構(gòu)空間懸鏈線半解析非線性有限元分析

索結(jié)構(gòu)空間懸鏈線半解析非線性有限元分析

懸索結(jié)構(gòu)主要表現(xiàn)為高度柔性和強(qiáng)幾何非線性。在索的有限單元法分析中最開始主要采用直線型桿法,忽略索的抗扭能力和抗彎能力,用不傳遞彎矩和扭距的許多相連直桿來模擬索。當(dāng)單元足夠多時(shí),所模擬的索單元就趨近于索的真實(shí)力學(xué)性能了。由于該法不能考慮垂度引起的非線性,因此,Ernst提出用修正材料的彈性模量考慮索垂度引起的非線性。為了追求高精度,一般還將索模擬成多節(jié)點(diǎn)的曲線索單元,如三節(jié)點(diǎn)二次多項(xiàng)式曲線單元,四節(jié)點(diǎn)等參索單元,五節(jié)點(diǎn)曲線索單元等。由于多節(jié)點(diǎn)索單元的前后處理不便,一些研究者探索將索模擬成二節(jié)點(diǎn)的曲線索單元。其曲線線型一般有拋物線和懸鏈線2種。拋物線法由于其假定的線形與實(shí)際索的線形略有不同,只是索單元的一種近似,適合于中小垂度的索分析。Z.H.Zhu基于空間曲梁理論,實(shí)現(xiàn)了松弛應(yīng)力和大垂度的索受力分析?；谟邢拊碚摰乃鲉卧獛缀醵际窃趹宜鹘Y(jié)構(gòu)的初始態(tài)已經(jīng)確定的基礎(chǔ)上建立的,在已確定的初始幾何上,基于總體Lagrangian描述建立位移函數(shù)和有限元方程,這樣建立起來的有限元方程用于懸索結(jié)構(gòu)荷載態(tài)的分析比較方便,但不能直接或不便用于懸索結(jié)構(gòu)初始幾何及內(nèi)力的確定。事實(shí)上,懸索結(jié)構(gòu)設(shè)計(jì)的關(guān)鍵就在于其初始幾何及內(nèi)力的確定。而上述有限元方法依賴的索初始形態(tài)通過重力作用下的平面懸鏈線或拋物線形確定。目前用懸鏈線方法進(jìn)行分析時(shí),主要采用重力作用下平面懸鏈線理論結(jié)果。對(duì)空間索單元的平面外剛度系數(shù),目前文獻(xiàn)中存在多種取值方法,不恰當(dāng)?shù)娜≈悼赡苁褂?jì)算收斂困難或收斂緩慢。上述推導(dǎo)與分析思路是在平面懸鏈線單元的基礎(chǔ)上增加剛度平面外剛度系數(shù)而得到。文中不依賴平面懸鏈線結(jié)果,直接推導(dǎo)了空間均布荷載作用下懸鏈線單元的解析解,并將其應(yīng)用于空間懸鏈線單元的有限單元分析。該方法不需要預(yù)知索的初始形態(tài),可以一次性精確求解在任意方向空間均布和集中荷載下索的位形及內(nèi)力。1在空間分布負(fù)載的功能下，尋求分析的解決方案1.1基本假設(shè)1)索在彈性階段工作,符合胡克定律;2)大位移小應(yīng)變變形;3)忽略索的抗彎能力。1.2未伸長(zhǎng)索上的荷載q1)索段的平衡微分方程在三維笛卡爾坐標(biāo)中,任意坐標(biāo)及方向的曲線索段如圖1所示,在3個(gè)坐標(biāo)軸方向ei(i=1,2,3)都承受著分布荷載q。三維索段變形后微段dS*的平衡方程為dΤdS*+q*=0,(1)dTdS?+q?=0,(1)式中,S*為A點(diǎn)到B點(diǎn)的有應(yīng)力長(zhǎng)度,T為張力矢量,用笛卡爾坐標(biāo)系的分量表達(dá)式為T=Tiei。矢量q*為伸長(zhǎng)后單位索長(zhǎng)上的分布荷載。為了使用方便,一般采用A點(diǎn)到B點(diǎn)的無應(yīng)力長(zhǎng)度S和未伸長(zhǎng)索上的荷載q,有q*S*=qS,(2)q?S?=qS,(2)假設(shè)qi是沿索長(zhǎng)均勻分布,則有:Τi=Τ0i-qiS(i=1,2,3),(3)Ti=T0i?qiS(i=1,2,3),(3)式中,T0i為初始點(diǎn)的張力矢量。2)物理方程對(duì)于符合胡克定律的材料,有應(yīng)變?chǔ)?T/EA,這里假設(shè)γ為小應(yīng)變,因此,有dS*dS=1+γ=1+ΤAE,(4)式中,T為索的張力合力。3)幾何協(xié)調(diào)條件張力矢量是沿著單位切線矢量方向的,能夠得到三維笛卡爾坐標(biāo)系下坐標(biāo)dxi和dS的關(guān)系,dxidS=dxidS*dS*dS=ΤiΤdS*dS,(5)聯(lián)合式(3)、(4)、(5),可得:xi=x0i+SEA(Τ0i-qiS2)+1Q3[Τ0m(Q2δim-qiqm)?lnQΤ-qjΤjQΤ0-qjΤ0j-qiQ(Τ-Τ0)](i,m,j=1,2,3)?(6)式中,下標(biāo)“0”表示索起始位置;系數(shù)Q是qi的量值,即:i=1,2,3;δim是克羅內(nèi)克符號(hào),即:δim={10(i=m)(i≠m)(i,m=1,2,3)。索伸長(zhǎng)后的長(zhǎng)度為S*=S(1+Τ2EA)+12EAQ3[QΤ0iqi(Τ0-Τ)+(Q2Τ20-Τ0iqiΤ0jqj)lnQΤ-Τ0kqk+Q2SQΤ0-Τ0mqm]。(7)假如兩端的力都已指定,則這些力必須滿足整體平衡方程,即式(3)。對(duì)于沿坐標(biāo)系方向的坐標(biāo)xi,可由在式(6)中令起始坐標(biāo)x0i=0得到。2空間負(fù)荷功能的支撐單元2.1有限元分析的方法對(duì)于任意的空間荷載q1、q2及q3代入式(6)可以推導(dǎo)出3個(gè)方向的索變形后長(zhǎng)度表達(dá)式。T0j(j=1,2,3)為起始點(diǎn)力,表示為-F1、-F2、-F3;Tj(j=1,2,3)為終點(diǎn)力,表示為F4、F5、F6。將式(6)展開,有Lx=-F1L0-12q21L0EA+q1(Τ1-Τ2)Q2-1Q3[-F1(q22+q23)+F3q2q3+F2q1q2]?lnuvLy=-F2L0-12q22L0EA+q2(Τ1-Τ2)Q2-1Q3[F2(q21+q23)-F3q1q3-F1q1q2]?lnuvLz=-F3L0-12q23L0EA+q3(Τ1-Τ2)Q2-1Q3[F3(q21+q22)+F1q1q3+F2q2q3]?lnuv}?(8)式中:Lx、Ly以及Lz分別表示3個(gè)方向索投影長(zhǎng)度;L0為無應(yīng)力長(zhǎng)度;F1、F2、F3分別為索兩端節(jié)點(diǎn)力(見圖2);其余參數(shù)意義見式(9)。u=QΤ2-q1F4-q2F5-q3F6v=QΤ1+q1F1+q2F2+q3F3Τ1=√F21+F22+F23Τ2=√F23+F24+F25Q=√q21+q22+q23}?(9)對(duì)于受任意均布荷載下的單根索,利用式(8)參考平面懸鏈線的迭代分析即可求得節(jié)點(diǎn)力。對(duì)于受有集中荷載作用的索,必須分成多個(gè)單元,針對(duì)每個(gè)單元,由式(8)分別對(duì)F1、F2、F3求偏導(dǎo)數(shù),可得柔度系數(shù)f11=?Lx?F1=-L0EA+F1q1Q2Τ1+F4q1Q2Τ2-q22+q23Q3lnuv+(Τ2q1-QF4Τ2u-Τ1q1+QF1Τ1v)w1Q3?f12=?Lx?F2=F2q1Q2Τ1+F5q1Q2Τ2+q1q2Q3lnuv+((Τ2q2-QF5)Τ2u-(Τ1q2+QF2)Τ1v)w1Q3?f13=?Lx?F3=F3q1Q2Τ1+F6q1Q2Τ2+q1q2Q3lnuv+((Τ2q3-QF6)Τ2u-(Τ1q3+QF3)Τ1v)w1Q3?f22=?Ly?F2=-L0EA+F2q2Q2Τ1+F5q2Q2Τ2-q21+q23Q3lnuv+((Τ2q2-QF5)Τ2u-(Τ1q2+QF2)Τ1v)w2Q3?f23=?Ly?F3=F3q2Q2Τ1+F6q2Q2Τ2+q2q3Q3lnuv+((Τ2q3-QF6)Τ2u-(Τ1q3+QF3)Τ1v)w2Q3?f33=?Lz?F3=-L0EA+F3q3Q2Τ1+F6q3Q2Τ2-q21+q22Q3lnuv+((Τ2q3-QF6)Τ2u-(Τ1q3+QF3)Τ1v)w3Q3?f21=?Ly?F1=f12;f31=?Lz?F1=f13?f32=?Lz?F2=f23?(10)式中:w1=F2q1q2-F1q22+F3q1q3-F1q23w2=F1q1q2-F2q21+F3q2q3-F2q23w3=F1q1q3-F3q21+F2q2q3-F3q22}?(11)則剛度矩陣為Κ=[f11f12f13f21f22f23f31f32f33]-1?(12)根據(jù)位移互等定理,任意空間均布荷載作用下空間索單元的切線剛度矩陣可表示為Κe=[Κ-Κ-ΚΚ]。(13)與平面索單元一樣,空間荷載作用下索單元的切線剛度矩陣也是對(duì)稱的。對(duì)于該剛度矩陣的求解,通常是預(yù)先假定一端節(jié)點(diǎn)力F1、F2、F3,根據(jù)式(3)可以確定F4、F5、F6,根據(jù)式(8)求得位移的增量Δd,進(jìn)而確定節(jié)點(diǎn)力的增量,反復(fù)迭代直至位移或者力的增量歐幾里德范數(shù)足夠小為止,通常情況下,只需要迭代幾步就可快速收斂。F1、F2、F3的初值可按下式F1=-q2Lx/2λ-q3Lx/2λF2=-q2(S-Lycoth(λ))F3=-q3(S-Lzcoth(λ))}?(14)式中:λ=√3(S2-L2y-L2z)/L2z-1;當(dāng)Lx等于0時(shí),λ=1e6;當(dāng)初始無應(yīng)力長(zhǎng)度S小于弦長(zhǎng)時(shí),λ=0.2。2.2節(jié)點(diǎn)荷載的疊加如果整根索受到空間均布荷載的作用,利用一個(gè)單元即可得到精確解,考慮到索的中間往往存在集中質(zhì)量和集中力,如減振器等,事實(shí)上整根索已經(jīng)不是懸鏈線線形,而是分段懸鏈線線形,因此,需要將索在有集中力處分段,將分段處的集中力和等效節(jié)點(diǎn)荷載疊加成總節(jié)點(diǎn)荷載。每個(gè)索單元兩端節(jié)點(diǎn)處的等效節(jié)點(diǎn)力可以通過對(duì)式(8)的迭代求得節(jié)點(diǎn)力的反號(hào)而得。ˉF={-F1,-F2,-F3,-F4,-F5,-F6}?(15)總等效節(jié)點(diǎn)荷載為等效節(jié)點(diǎn)力與節(jié)點(diǎn)集中荷載之和。如果需要考慮溫度荷載,則可以通過對(duì)(7)式加上溫度變形實(shí)現(xiàn)求解。2.3單元的節(jié)點(diǎn)位移增量計(jì)算由于上述推導(dǎo)均在全局笛卡爾坐標(biāo)系下完成,因此不需要進(jìn)行坐標(biāo)系的轉(zhuǎn)換,通過對(duì)單元的剛度矩陣以及等效結(jié)點(diǎn)荷載向量集總可直接得到全局坐標(biāo)系下的總剛和節(jié)點(diǎn)力。由虛功原理可以得到非線性有限元平衡迭代方程[Κ]t+Δt{Δd}=t+Δt{F}?(16)式中:[Κ]為全局坐標(biāo)系下t+Δt時(shí)刻的切線剛度矩陣,{Δd}為t+Δt時(shí)段索單元的節(jié)點(diǎn)位移增量,t+Δt{F}為t+Δt步迭代過程中的等效節(jié)點(diǎn)荷載向量,通過迭代由(14)求得。索的位形可以通過對(duì)節(jié)點(diǎn)坐標(biāo)的更新而調(diào)整。上式中剛度矩陣[Κ]和索節(jié)點(diǎn)位形有關(guān),等效節(jié)點(diǎn)荷載t+Δt{F}和索節(jié)點(diǎn)位形有關(guān)。因此,上述該方程采用雙重平衡迭代法進(jìn)行求解。1)對(duì)于任意給定空間索的位形,迭代求解每一個(gè)單元的剛度矩陣和等效節(jié)點(diǎn)力;2)組裝成全局坐標(biāo)系下總剛度切線矩陣和等效節(jié)點(diǎn)力,求得節(jié)點(diǎn)位移以及等效節(jié)點(diǎn)荷載向量;3)求解式(16),求得t+Δt時(shí)段索的節(jié)點(diǎn)位移增量,更新索的坐標(biāo)。若位移增量的歐幾里德范數(shù)足夠小,則收斂,可得到節(jié)點(diǎn)位移和節(jié)點(diǎn)力。否則,重復(fù)第(1)步迭代過程,直至收斂。上述平衡方程求解步驟在每一次迭代過程中的等效節(jié)點(diǎn)荷載向量和切線剛度矩陣都是基于彈性懸鏈線精確表達(dá)式一次求解,沒有采用UL列式的增量求解遞推關(guān)系,避免了對(duì)等效節(jié)點(diǎn)荷載的增量求解,減少了誤差累積。計(jì)算表明,單元數(shù)較少的情況下能夠很快收斂至平衡狀態(tài);若單元數(shù)較多,為避免計(jì)算發(fā)散,可對(duì)一次性求解的等效節(jié)點(diǎn)荷載向量乘以一個(gè)較小的系數(shù),這樣可確保收斂,但收斂速度會(huì)減慢。3計(jì)算與分析3.1單元?jiǎng)澐值姆椒ǜ鶕?jù)上述理論,編制程序,首先驗(yàn)證一個(gè)經(jīng)典的平面索的算例。如圖3所示的各種索的形狀。索的無應(yīng)力長(zhǎng)度L0=100m,溫度變化T=100℃,材料幾何特性EA=3×107N,自重W=1N/m,線膨脹系數(shù)α1=0.65×10-5,輸入每個(gè)J點(diǎn)的位置,用程序計(jì)算劃分一個(gè)單元即可得到每個(gè)位置J端點(diǎn)處索的張力,計(jì)算的結(jié)果與文獻(xiàn)中的完全一樣。將圖中索劃分多個(gè)單元,計(jì)算位形和節(jié)點(diǎn)力和一個(gè)單元完全一致,計(jì)算精度與單元?jiǎng)澐謹(jǐn)?shù)無關(guān),可見求解是完全精確的。3.2基于空間懸鏈線單元非線性有限元分析目前索受空間風(fēng)荷載和自重荷載的算例較少,為了說明算法的有效性,采用商用程序ANSYS進(jìn)行驗(yàn)證。計(jì)算模型同樣為一根空間索,材料特性和上例一樣,索的初始無應(yīng)力長(zhǎng)度為100N,兩端節(jié)點(diǎn)坐標(biāo)分別為I(0,0,0)和J(80,0,0)。ANSYS采用類似于文獻(xiàn)的兩節(jié)點(diǎn)直線索單元,并不能直接得到懸鏈線的初始內(nèi)力和位形。因此,ANSYS計(jì)算結(jié)果是在確定了懸鏈線的位形之后施加風(fēng)荷載進(jìn)行求解的。由于ANSYS采用了直線兩節(jié)點(diǎn)單元逼近曲線索,因此劃分了多達(dá)500個(gè)單元。文中算法將單元分成4個(gè)單元。圖4顯示了索中點(diǎn)坐標(biāo)的收斂情況。索的自重荷載和風(fēng)荷載均為1N/m。圖中顯示索中點(diǎn)的初始坐標(biāo)為(5,-40,0),在荷載作用下能快速收斂到最終位形(40.0,-18.7,-18.7)。大量分析表明,利用基于空間懸鏈線單元,當(dāng)索的單元原長(zhǎng)和荷載確定了以后,索單元的中間節(jié)點(diǎn)可以在任意位置,均可快速迭代至最終平衡位置。圖5表示不同風(fēng)速下索中點(diǎn)的位形變化。其中,自重荷載為1N/m,風(fēng)速荷載為0.25～2.0N/m;可見,在風(fēng)荷載下,索的位形具有明顯的非線性特征。圖6表示不同風(fēng)偏角下索中點(diǎn)的位形變化。自重荷載為1N/m,風(fēng)速荷載為1N/m,風(fēng)偏角為0°～90°?？梢钥闯霾煌娘L(fēng)偏角導(dǎo)致索中點(diǎn)位形呈非線性變化。圖7表示不同風(fēng)偏角下索兩端節(jié)點(diǎn)力的變化。隨風(fēng)偏角增大,索一端節(jié)點(diǎn)張力增大,另一端減小,呈非線性變化。圖8為索中間施加不同集中力下位形變化。自重荷載為1N/m,風(fēng)速荷載為1N/m,集中力為1～10N?？梢钥闯?集中力對(duì)索的位形也有較大影響。從上述計(jì)算結(jié)果可見,基于空間懸鏈線單元非線性有限元分析可以一次性實(shí)現(xiàn)索的張力和位形的迭代,可以求解任意方向空間荷載下索的位形和張力。由于本算例為一個(gè)小垂度索,因此,計(jì)算結(jié)果和ANSYS計(jì)算結(jié)果比較接近。但是基于空間懸鏈線單元的有限元分析需要更少的單元,不需要初始位形,不局限中小垂度范圍,效率更高,結(jié)果