梁浩博,劉小會(huì),2,楊曙光,閔光云,伍 川
(1.重慶交通大學(xué) 土木工程學(xué)院,重慶400074;2.重慶交通大學(xué) 省部共建山區(qū)橋梁及隧道工程國家重點(diǎn)實(shí)驗(yàn)室,重慶400074;3.國網(wǎng)河南省電力公司 電力科學(xué)研究院,鄭州450052)
基于經(jīng)濟(jì)發(fā)展的需求,國內(nèi)外很多工程領(lǐng)域都廣泛使用大跨度的索結(jié)構(gòu),例如橋梁工程領(lǐng)域的斜拉橋和懸索橋、電氣工程領(lǐng)域的高壓輸電線路等。然而索結(jié)構(gòu)因其自身重量輕、柔度大以及低阻尼的特點(diǎn),很容易在外力荷載作用下出現(xiàn)大幅度的振動(dòng),呈現(xiàn)出非常明顯的非線性現(xiàn)象。索的這種大幅振動(dòng)極易導(dǎo)致整個(gè)工程結(jié)構(gòu)出現(xiàn)損壞。因此,關(guān)于索結(jié)構(gòu)的非線性動(dòng)力學(xué)研究是個(gè)很重要的課題。
關(guān)于索的非線性動(dòng)力學(xué)研究早在18 世紀(jì)就開展了。隨后,Routh[1]推導(dǎo)出懸索的運(yùn)動(dòng)方程并得到懸索固有頻率。Irvine 等[2]在1974年推導(dǎo)出水平索固有振動(dòng)方程式,并對(duì)拉索的模態(tài)、頻率等參數(shù)進(jìn)行系統(tǒng)的分析研究。文獻(xiàn)[3]基于哈密頓變分原理并采用Galerkin 離散法推導(dǎo)出單檔覆冰懸索運(yùn)動(dòng)方程。文獻(xiàn)[4]引入無量綱參數(shù)推導(dǎo)出受溫度變化影響的懸索面內(nèi)非線性偏微分方程,利用Galerkin 法得到離散后的微分方程。文獻(xiàn)[5]同樣采用Galerkin法,建立一個(gè)兩自由度懸索離散模型,對(duì)懸索的1:2內(nèi)共振進(jìn)行了研究。文獻(xiàn)[6]采用Galerkin方法得到懸索的時(shí)滯微分方程,以對(duì)其主共振進(jìn)行分析。
由于大跨度輸電線路實(shí)質(zhì)上也屬于索結(jié)構(gòu),因此關(guān)于輸電線路的學(xué)術(shù)研究同樣可以為拉索研究帶來指導(dǎo)性意見。Rega等[7]引入無量綱參數(shù)并基于哈密頓變分原理推導(dǎo)出輸電線面內(nèi)非線性偏微分方程,接著應(yīng)用Galerkin 法和多尺度法求解該方程。汪峰等[8]基于Galerkin 模態(tài)截?cái)喾ê投喑叨确ㄍ茖?dǎo)出斜拉索面內(nèi)振動(dòng)方程并得到幅頻響應(yīng)方程,分析了梯度溫度場(chǎng)等因素對(duì)振動(dòng)的影響。文獻(xiàn)[9]基于懸鏈法和熱應(yīng)力理論推導(dǎo)出考慮溫度效應(yīng)的覆冰導(dǎo)線偏微分方程,接著采用Galerkin 法和多尺度法將其離散和求解。受限于非線性理論的發(fā)展,上述眾多學(xué)者關(guān)于理論方程的求解都是基于模態(tài)疊加法對(duì)系統(tǒng)進(jìn)行了線性化處理,最終得到結(jié)構(gòu)的近似分析解。該處理方法雖然在一定程度上簡(jiǎn)化了計(jì)算,但對(duì)于這些具有幾何非線性的結(jié)構(gòu)是否可以直接使用,會(huì)不會(huì)對(duì)結(jié)果造成誤差等等,這些問題仍有待商榷。
國內(nèi)學(xué)者趙躍宇等[10]以懸索為例,并引入新的無量綱參數(shù)得到無量綱懸索面內(nèi)運(yùn)動(dòng)方程,對(duì)結(jié)構(gòu)非線性動(dòng)力學(xué)中常用的Galerkin法和直接法進(jìn)行了對(duì)比,發(fā)現(xiàn)離散法并不適用于所有結(jié)構(gòu)非線性動(dòng)力學(xué)的求解。文獻(xiàn)[11]中,作者采用小參數(shù)法得到久期項(xiàng)導(dǎo)致的懸索振動(dòng)形態(tài),并與Galerkin 離散法得到的相應(yīng)振動(dòng)形態(tài)進(jìn)行對(duì)比,最終驗(yàn)證了小參數(shù)法的可行性。然而上述文章均是以兩種算法的對(duì)比為主要目的,針對(duì)Galerkin 離散法本身的適用性并未作出詳細(xì)區(qū)分。
另外,Rega 等[7]和趙躍宇等[10]分別引入了不同的無量綱參數(shù)將運(yùn)動(dòng)方程無量綱化,但卻并未有學(xué)者對(duì)這兩種無量綱參數(shù)對(duì)系統(tǒng)非線性的影響以及兩種無量綱方法是否適用于所有結(jié)構(gòu)非線性系統(tǒng)做詳細(xì)的研究。
因此,本文基于簡(jiǎn)單的單跨拉索力學(xué)模型以及振動(dòng)理論,針對(duì)Galerkin 離散法和兩種無量綱方法的適用性進(jìn)行詳細(xì)討論。
設(shè)單跨拉索兩端處于同一水平高度,在自重作用下位于如圖1所示的XOY平面內(nèi),X軸沿拉索弦向方向,Y軸沿拉索豎向垂度方向。拉索跨徑為l,其變形可以分為兩部分:
圖1 單跨水平拉索力學(xué)模型
(1) 拉索在自重作用下處于靜平衡狀態(tài),垂度為d;
(2)拉索在面內(nèi)外部荷載作用下產(chǎn)生處于XOY平面內(nèi)的位移。
考慮到拉索在受到幾何非線性和三向耦合時(shí)所引發(fā)的動(dòng)力學(xué)問題很復(fù)雜,本文將模型簡(jiǎn)化,做出如下假定:
(2)本構(gòu)關(guān)系符合胡克定律;
(3) 不考慮彎曲、扭轉(zhuǎn)以及剪切,忽略其軸向運(yùn)動(dòng);
(4)不考慮張力沿拉索軸向的變化。
圖1所示理論模型中,考慮沿拉索Y向施加均布簡(jiǎn)諧激勵(lì)荷載Pcos(Ωt),P為激勵(lì)幅值,Ω為簡(jiǎn)諧激勵(lì)的頻率。取圖1中拉索處于靜平衡狀態(tài)下的微段ds進(jìn)行研究,可以得到單跨拉索的面內(nèi)運(yùn)動(dòng)方程:
式中:H為拉索在自重作用下達(dá)到靜平衡狀態(tài)時(shí)的水平張力,v為沿拉索面內(nèi)Y向的位移,E為拉索彈性模量,A為拉索橫截面積,m為拉索單位長(zhǎng)度質(zhì)量,μ為黏性阻尼系數(shù)。
引入無量綱參數(shù):
將上述參數(shù)代入式(1),得到拉索面內(nèi)無量綱運(yùn)動(dòng)方程:
為了方便后文書寫,式(3)中各無量綱參數(shù)的上劃線“-”均已省略(后文同)。式(3)中“v′”表示對(duì)無量綱時(shí)間t求1 階導(dǎo),“v″”表示對(duì)無量綱位移x求2階導(dǎo)。
由于v是關(guān)于時(shí)間t和位移x的參數(shù),根據(jù)Galerkin理論,可將其分離變量為:
其中:φ(x)為模態(tài)函數(shù),q(t)為時(shí)間函數(shù)。
將式(4)代入式(3)中,并在方程兩端乘以模態(tài)函數(shù)φ(x),接著在x∈[0,1]內(nèi)進(jìn)行積分,可得到:
式中:
其中:λ2為Irvine 參數(shù),見式(7)。該參數(shù)可以反映拉索的垂度效應(yīng)。λ的值越大,拉索的垂度越大。ωn即為無量綱化后的頻率,ω表示面內(nèi)自振等效圓頻率,f為拉索自振頻率。式(5)中q上方的點(diǎn)“˙”和“˙˙”,分別表示對(duì)無量綱時(shí)間t求1階導(dǎo)和2階導(dǎo)。
g是重力加速度。
引入另一種無量綱參數(shù):
式中:δ為垂跨比,張力H可表示為:
取其無量綱靜態(tài)垂度曲線為:
將上述無量綱參數(shù)代入式(1),可得拉索面內(nèi)無量綱運(yùn)動(dòng)方程:
式中:
同樣,該式通過模態(tài)疊加法可得到:
式中:
由于兩種無量綱方法主要區(qū)別在于豎向位移v所除是d還是l,故下文均以無量綱d″和無量綱l″來分別表示兩種不同的無量綱方法。
研究拉索在外激勵(lì)作用下產(chǎn)生受迫振動(dòng)的現(xiàn)象,首先要研究其頻率和對(duì)應(yīng)模態(tài)。其中模態(tài)分為兩種形式,對(duì)稱模態(tài)與反對(duì)稱模態(tài)。本文主要研究對(duì)稱模態(tài)。
取方程(3)的線性部分(忽略阻尼項(xiàng))及邊界條件:
令v=φ(x)eiωt并將其代入式(15)中,結(jié)合邊界條件可得拉索面內(nèi)模態(tài)函數(shù):
取正交模態(tài):
ωn可由超越方程式(19)求得:
觀察式(1)的拉索面內(nèi)運(yùn)動(dòng)方程發(fā)現(xiàn),方程的形式較為復(fù)雜,且變量較多,不利于進(jìn)行數(shù)值求解。而引入無量綱參數(shù)就可以減少變量數(shù)目,便于后續(xù)求解。同時(shí),由于對(duì)參數(shù)進(jìn)行了無量綱化,更方便做不同工況下的參數(shù)對(duì)比。
但是可以注意到兩種方法關(guān)于面內(nèi)豎向位移v的處理方法是不同的,這是否會(huì)導(dǎo)致兩種方法所得結(jié)果出現(xiàn)一定的差異,本文考慮將這兩種方法做對(duì)比(拉索物理參數(shù)均取自文獻(xiàn)[12],詳見表1)。
表1 拉索物理參數(shù)
文中系統(tǒng)的黏性阻尼系數(shù)μ均采用瑞利阻尼,具體計(jì)算方法如下:
取α為0.02,β為0,將其代入公式中,即可得到阻尼比ξ。而阻尼μ則可由公式μ=2ξωm得到。
由式(7)可知:受到張力和垂度的影響,不同參數(shù)對(duì)應(yīng)的拉索非線性強(qiáng)度各不相同,故本文選取λ為等于1π、3π、7π、9π等4種工況進(jìn)行研究。4種工況下的具體參數(shù)見表2。
表2 不同Irvine參數(shù)對(duì)應(yīng)的物理參數(shù)
將表1 和表2 具體參數(shù)導(dǎo)入MATLAB 程序中,即可運(yùn)行得到系統(tǒng)時(shí)間歷程曲線,如圖2所示。
2018年1—8月,全國國有及國有控股企業(yè)(本月報(bào)所稱全國國有及國有控股企業(yè),包括中央管理企業(yè)、中央部門和單位所屬企業(yè)以及36個(gè)省、自治區(qū)、直轄市、計(jì)劃單列市的地方國有及國有控股企業(yè),不含國有金融類企業(yè)。以下簡(jiǎn)稱國有企業(yè))經(jīng)濟(jì)運(yùn)行態(tài)勢(shì)良好。償債能力和盈利能力比上年同期均有所提升,利潤增幅高于收入10.4個(gè)百分點(diǎn),鋼鐵等行業(yè)利潤增幅較高。
通過觀察圖2可知:無量綱d方法相較于無量綱l方法所得幅值略大,尤其在λ=1π 時(shí),差距較為明顯,但隨著λ值的增大,這種差距逐漸縮小,直至忽略不計(jì);同時(shí)由表3可知:兩種無量綱方法均會(huì)使系統(tǒng)非線性系數(shù)增大,尤其以無量綱l方法增加最大。
圖2 兩種無量綱方法時(shí)間歷程曲線
表3 不同Irvine參數(shù)對(duì)應(yīng)的非線性項(xiàng)
由于拉索運(yùn)動(dòng)方程比較復(fù)雜,一般難以求得其精確解,且拉索振動(dòng)屬于弱非線性問題,因此可采取近似的解析方法求解。傳統(tǒng)的解析方法主要包括:漸進(jìn)法、多尺度法、平均法等,而多尺度法是其中應(yīng)用得最廣泛的方法。本文亦采用多尺度法來分析系統(tǒng)響應(yīng)。
為了滿足多尺度法的求解條件,首先對(duì)式(1)直接采用模態(tài)疊加法進(jìn)行離散。接著,為了使非線性項(xiàng),外部激勵(lì)和阻尼出現(xiàn)在同階中,引入無量綱小量ε,經(jīng)過化簡(jiǎn)后,式(1)即可轉(zhuǎn)化為:
式中:
上式(21)中的等效圓頻率、平方和立方非線性項(xiàng)系數(shù)以及外激勵(lì)中所包含的I的取值可表示為:
為了方便書寫,后文已將式(20)和式(21)中的“*”忽略。
引入以下參數(shù):
其中:σ為調(diào)諧參數(shù)。
將式(23)代入式(20)中,并比較ε同次冪項(xiàng)的系數(shù),可得:
式(24)的解可用指數(shù)形式表示為:
此時(shí),A可寫成極坐標(biāo)形式:
式中:a表示振幅,β表示相位。
將其代入式(25)中,即可得到:
消除式(28)中的久期項(xiàng),可得:
將上式進(jìn)行化簡(jiǎn)可得:
并將A代入式(30)可得到:
分離實(shí)部和虛部可得:
引入:
方程即可轉(zhuǎn)化為:
令D1a=0,D1γ=0可得出以下方程:
上述方程兩端各自平方并相加即可得幅頻響應(yīng)方程:
由式(39)并結(jié)合表1中的拉索物理參數(shù),即可得出拉索的幅頻響應(yīng)曲線,見圖3。
圖3 幅頻響應(yīng)曲線
如圖3 所示:橫坐標(biāo)σ表示調(diào)諧參數(shù)(無量綱),縱坐標(biāo)a表示幅值(單位為m)。實(shí)線表示穩(wěn)定值,虛線表示不穩(wěn)定值。
結(jié)合曲線和式(23)所引入的小量可知:當(dāng)調(diào)諧參數(shù)σ=0時(shí),外激勵(lì)頻率與系統(tǒng)固有圓頻率相等,此時(shí)拉索產(chǎn)生主共振。同時(shí),隨著λ增大,系統(tǒng)非線性不斷增強(qiáng),曲線出現(xiàn)跳躍現(xiàn)象,幅值也出現(xiàn)明顯的衰減。
當(dāng)σ=0時(shí),將參數(shù)代入式(39)中可得到系統(tǒng)穩(wěn)定后的理論幅值為:(1π)0.589 27 m;(3π)0.210 062 m;(7π)0.062 395 m;(9π)0.051 052 m。
將該理論幅值與1.3 節(jié)中兩種無量綱方法穩(wěn)定后的幅值對(duì)比后發(fā)現(xiàn):無量綱d方法與理論解的誤差依次為:22 %(1π),3.797 %(3π)、2.124 %(7π)、0.878%(9π);無量綱l方法與理論解的誤差依次為:4.323 %(1π)、3.688 %(3π)、2.624 %(7π)、3.005 %(9π),總體來看無量綱l方法所得結(jié)果更接近理論解。
由于拉索在受到覆冰、外部激勵(lì)等因素的影響時(shí),易產(chǎn)生較大的振動(dòng)幅度。因此不易通過實(shí)驗(yàn)來獲取數(shù)據(jù),以驗(yàn)證理論,而這時(shí)采用有限元方法數(shù)值模擬拉索的振動(dòng)便成為了研究拉索受迫振動(dòng)的重要手段[13-14]。
然而在實(shí)際工程中,拉索振動(dòng)時(shí)具有明顯的幾何非線性效應(yīng),通過有限元模擬所得的結(jié)果是否會(huì)考慮到這種非線性效應(yīng)的影響?是否會(huì)產(chǎn)生誤差?因此,本文通過建立有限元模型,來模擬拉索的受迫振動(dòng),并將結(jié)果與數(shù)值解作對(duì)比。
采用桿單元模擬單跨拉索,網(wǎng)格數(shù)量為100。由于本文僅討論研究拉索在垂直方向上的運(yùn)動(dòng),故在拉索兩端施加邊界條件,將其完全固定。
采用參數(shù)化建模后,第一步:在拉索上施加重力荷載,以對(duì)其進(jìn)行非線性靜力變形分析。第二步:頻率分析,提取拉索的前30 階固有頻率及模態(tài),并與基于理論所得的固有頻率作對(duì)比,見表4。第三步,沿拉索豎向施加均勻外力激勵(lì),以使其產(chǎn)生面內(nèi)的豎向振動(dòng),力的大小設(shè)置為0.7 N(相當(dāng)于施加大小為5.6 m/s風(fēng)速的風(fēng)荷載)。阻尼項(xiàng)選取瑞利阻尼。
表4 不同垂度下拉索第一階面內(nèi)對(duì)稱模態(tài)自振頻率對(duì)比
重點(diǎn)考慮最危險(xiǎn)的情形,即拉索產(chǎn)生主共振。將表1和表2的拉索物理參數(shù)分別輸入ABAQUS軟件和MATLAB程序中,即可運(yùn)行得到系統(tǒng)時(shí)間歷程曲線。由第2節(jié)在本文方法基礎(chǔ)上分別使用兩種無量綱參數(shù)所得數(shù)值解與理論解相比的誤差可知:無量綱l方法所得結(jié)果與理論解吻合良好。限于篇幅,本文僅給出使用無量綱l方法時(shí)本文方法所得數(shù)值解與有限元解的對(duì)比結(jié)果,如圖4至圖7所示。
圖4 (λ=1π)拉索面內(nèi)對(duì)稱模態(tài)共振時(shí)間歷程曲線
圖5 (λ=3π)拉索面內(nèi)對(duì)稱模態(tài)共振時(shí)間歷程曲線
圖6 (λ=7π)拉索面內(nèi)對(duì)稱模態(tài)共振時(shí)間歷程曲線
圖7 (λ=9π)拉索面內(nèi)對(duì)稱模態(tài)共振時(shí)間歷程曲線
如圖4、圖5所示:當(dāng)λ取1π、3π時(shí),有限元解出現(xiàn)往幅值為0的水平線上方偏移的現(xiàn)象,而本文方法的計(jì)算結(jié)果則出現(xiàn)向下偏離的現(xiàn)象(偏移量見表5)。造成該現(xiàn)象的主要原因是由于系統(tǒng)中平方非線性項(xiàng)的存在,如:式(5)中的c2、式(13)中的γn以及式(20)中的εc2項(xiàng),具體數(shù)值見表3;該項(xiàng)越大,偏上或偏下的程度也就越大。同時(shí),由于λ取1π、3π 時(shí)系統(tǒng)的振動(dòng)幅值較大,其偏移也就相對(duì)更加明顯。當(dāng)λ=7π和9π時(shí),由于本身幅值很小,加之平方非線性項(xiàng)在7π 和9π 時(shí)逐步減小,故其偏移程度也相對(duì)很小。
表5 共振時(shí)拉索的時(shí)間歷程曲線偏移量
如圖4至圖7所示并結(jié)合表6、表7可以發(fā)現(xiàn):本文方法所得數(shù)值解與上節(jié)幅頻響應(yīng)函數(shù)計(jì)算得出的理論解吻合良好,而有限元解與本文方法所得數(shù)值解相比則出現(xiàn)了不同程度的幅值差距。這種差距呈現(xiàn)出隨著λ增大而增大的趨勢(shì)。
表6 無量綱l方法所得數(shù)值解與理論解對(duì)比
表7 有限元解與理論解對(duì)比
如圖6、圖7 所示:本文方法所得數(shù)值解較之有限元解更大,其原因在于當(dāng)使用有限元軟件進(jìn)行模擬時(shí),能夠真實(shí)的體現(xiàn)非線性動(dòng)力響應(yīng),而本文方法得到的數(shù)值解是采用了Galerkin離散法近似化的運(yùn)動(dòng)方程,這種近似化的方法弱化了非線性。綜合上述系統(tǒng)非線性越強(qiáng),幅值就越小的結(jié)論,最終出現(xiàn)了圖6、圖7所示結(jié)果。
就實(shí)際情況而言,拉索因發(fā)生共振而造成破壞的情況極少,大多是在覆冰和風(fēng)載兩種主要影響因素下所造成的非共振區(qū)的大幅度振動(dòng)。因此本節(jié)通過將施加在拉索上的外激勵(lì)頻率(ω=4 rad/s)偏離其固有圓頻率并采取與上一節(jié)相同的模式來進(jìn)行對(duì)比。其中,λ=1π、3π、7π、9π 時(shí)頻率偏離大小分別為:0.501 9 rad/s、-0.591 6 rad/s、0.137 7 rad/s、0.433 0 rad/s(外激勵(lì)頻率減去固有圓頻率)。
如圖9 至圖12 所示:其部分現(xiàn)象與共振區(qū)現(xiàn)象一致。而圖9所示結(jié)果出現(xiàn)了有限元解與本文方法所得數(shù)值解差距較大的現(xiàn)象。初步結(jié)合表2 判斷,造成該現(xiàn)象的原因是所施加的外荷載激勵(lì)頻率小于λ=3π時(shí)的固有圓頻率(ω=4.59 rad/s)。為驗(yàn)證該判斷,本文取ω=4.8 rad/s 時(shí),采用無量綱l方法進(jìn)行再次對(duì)比,結(jié)果如圖12所示,表明判斷正確。
圖8 (λ=1π)拉索面內(nèi)對(duì)稱模態(tài)非共振時(shí)間歷程曲線
圖9 (λ=3π)拉索面內(nèi)對(duì)稱模態(tài)非共振時(shí)間歷程曲線
圖10 (λ=7π)拉索面內(nèi)對(duì)稱模態(tài)非共振時(shí)間歷程曲線
圖11 (λ=9π)拉索面內(nèi)對(duì)稱模態(tài)非共振時(shí)間歷程曲線
圖12 (λ=3π)拉索面內(nèi)對(duì)稱模態(tài)非共振(ω=4.8 rad/s)時(shí)間歷程曲線
由圖(12)可以發(fā)現(xiàn):當(dāng)所施加外部激勵(lì)頻率為4 rad/s 時(shí),圖中所呈現(xiàn)出的現(xiàn)象并無明顯規(guī)律。分析可知:在非共振區(qū),本文方法所得數(shù)值解與有限元解的差距,主要取決于外激勵(lì)頻率的大小與其圓頻率的差距。外激勵(lì)頻率與圓頻率的差距越大,本文方法所得數(shù)值解與有限元解的差距也就越大。同時(shí),隨著λ值的變大,其時(shí)間歷程曲線就相對(duì)不規(guī)則,甚至出現(xiàn)相位差,如圖11所示。
本文通過應(yīng)用Galerkin離散法得出拉索無量綱面內(nèi)振動(dòng)方程,并運(yùn)用多尺度法得到系統(tǒng)的幅頻響應(yīng)曲線,最后分別通過ABAQUS 和MATLAB 進(jìn)行了4 種工況的數(shù)值模擬和有限元分析對(duì)比,主要得出以下結(jié)論:
(1)隨著系統(tǒng)非線性的增強(qiáng),幅值逐漸減小,而有限元解與理論解的差距則逐漸增大,說明Galerkin離散法對(duì)于弱非線性系統(tǒng)適用性更高,系統(tǒng)非線性越強(qiáng),誤差就越大。
(2)系統(tǒng)的平方非線性項(xiàng)會(huì)造成有限元解和本文方法所得結(jié)果分別產(chǎn)生向上和向下的偏移。
(3)隨著λ值的增大,兩種無量綱方法差距微小,總體而言經(jīng)過無量綱l方法得到的運(yùn)動(dòng)方程所得的系統(tǒng)時(shí)間歷程曲線更接近理論解;為了使結(jié)果更精確,當(dāng)系統(tǒng)非線性強(qiáng)時(shí)可采用無量綱d方法,系統(tǒng)非線性弱時(shí)則采用無量綱l方法較好。