朱進波， 鄭史雄， 唐 煜， 郭俊峰

(1.西南交通大學 土木工程學院 橋梁工程系，成都 610031；2.西南石油大學 土木工程與建筑學院，成都 610500)

橋梁顫振是一種自激發(fā)散振動，該振動現(xiàn)象會導致橋梁整體毀滅性的破壞。當前的顫振理論體系[1]已基本可以準確預測橋梁的顫振臨界風速。該風速預測方法是建立在結(jié)構(gòu)微幅振動的線性顫振理論框架[2-3]下，即橋梁顫振導數(shù)與斷面振幅無關(guān)，求得的顫振臨界風速具有明確的數(shù)值意義。在該風速后橋梁將直接作發(fā)散運動，在該體系之下并沒有研究顫振后狀態(tài)。

近年來隨著材料科學、結(jié)構(gòu)分析技術(shù)及施工方法的進步，橋梁結(jié)構(gòu)逐步向大跨、輕柔化發(fā)展，懸索橋逐漸成為大跨度橋型的首選。而橋梁跨度的增長使得橋梁剛度和阻尼急劇下降，橋梁結(jié)構(gòu)對風的敏感性增加，想要滿足顫振設(shè)防規(guī)定需要大幅增加建橋成本，可以預見風對大跨度橋梁結(jié)構(gòu)[4-6]的作用將逐漸成為提高橋梁跨徑的主要制約因素之一。

基于風洞試驗研究，Scanlan等將顫振導數(shù)表示的非定常氣動自激力引入振動平衡方程中，摒棄了流固耦合的直接求解，建立了經(jīng)典的顫振框架模型，采用半逆解法進行顫振求解。Wilde等[7-8]采用二自由度的狀態(tài)空間復模態(tài)特征值求解方法(Complex Eigenvalue Analysis, CEVA)，其基本思想是將二維顫振問題轉(zhuǎn)變?yōu)榍蠼鈴吞卣髦祮栴}，以某階模態(tài)下的阻尼比為零時，作為系統(tǒng)發(fā)散的依據(jù)，對應的風速為顫振臨界風速?；谠摽蚣芟碌慕饩哂忻鞔_的數(shù)學意義，當風速超過這一臨界點橋梁斷面的振幅將指數(shù)級增加而出現(xiàn)發(fā)散，可實際現(xiàn)象中卻并非如此。朱樂東等[9]在風洞試驗中發(fā)現(xiàn)，一些鈍體斷面顫振發(fā)生后并不會出現(xiàn)如線性理論預測的那樣，振幅按指數(shù)級迅速增長，而是由于自激力的非線性效應穩(wěn)定在一定幅值狀態(tài)。學者們開始關(guān)注橋梁進入顫振后的真實運動狀態(tài)，顫振研究逐步走向精細化。

想要弄清橋梁顫振后的特征并非易事，徐旭等[10]建立一個關(guān)于純扭轉(zhuǎn)顫振的非線性氣動力模型，研究橋梁顫振的非線性穩(wěn)定性問題；張朝貴[11]引入范德波爾方程來試圖解釋軟顫振現(xiàn)象；王騎等[12-13]將橋梁顫振形態(tài)轉(zhuǎn)變?yōu)閱巫杂啥扰まD(zhuǎn)方程，分別描述了極限環(huán)與軟顫振現(xiàn)象，反映了阻尼比與氣動力的非線性特性。盡管許多學者做了大量的研究工作，鑒于非線性問題的復雜性，對橋梁顫振后的特征狀態(tài)依然沒有定論。目前對于非線性自激力的探索主要集中在氣動力隨振幅的非線性變化規(guī)律與氣動本身的高次諧波特性兩個方面。

對于橋梁顫振后的振動狀態(tài)，最為關(guān)心的便是橋梁振動幅值與風速的關(guān)系。張彥等[14-15]分別從數(shù)值模擬和風洞實驗方面對不同振幅下橋梁斷面自激力進行研究。唐煜[16]基于傳統(tǒng)顫振耦合理論，發(fā)展出非線性二維兩自由度耦合顫振分析方法。

本文基于氣動自激力隨振幅響應變化這一思想，依托顫振理論框架和求解思路，提出更加規(guī)范、精細的考慮顫振幅值因素的二自由度復模態(tài)特征值求解方法。基于此方法，提出非線性顫振幅值響應搜索程序，一定程度上揭示了橋梁顫振后的特征狀態(tài)。

1 流線型箱梁斷面的顫振研究方法

隨著大跨度橋梁不斷的投入建設(shè)，為滿足其氣動穩(wěn)定性需求，流線型箱梁斷面型式逐漸成為大跨懸索橋主梁斷面型式的首選。目前主要采用直接風洞試驗法和基于風洞試驗識別參數(shù)的理論分析法來檢測橋梁的顫振穩(wěn)定性能。如今伴隨著流體力學與計算機硬件的發(fā)展，數(shù)值風洞識別氣動參數(shù)的方法逐漸被接受與認可，其純理論的計算方法有望代替物理風洞試驗方法。該方法可分為三個部分，即自激氣動力模型、顫振導數(shù)數(shù)值識別和二維耦合顫振分析方法，全部過程都是數(shù)值計算。

選取Scanlan等的經(jīng)典顫振理論框架作為自激氣動模型

(1)

(2)

基于此方法，本文以某大跨度懸索橋為研究對象，研究其顫振幅值特性，該橋梁結(jié)構(gòu)基本參數(shù)如表1，流線型箱梁斷面尺寸如圖1所示。

圖1 主梁斷面圖Fig.1 Section diagram of the main beam

表1 橋梁結(jié)構(gòu)基本參數(shù)Tab.1 Basic parameters of the bridge structure

2 氣動參數(shù)識別

基于Fluent軟件，建立二維數(shù)值風洞模型，通過著名的N-S方程實現(xiàn)對流體時間域和空間域的控制，采用SSTk-ω湍流模型，實現(xiàn)數(shù)值模擬。應用強迫振動法[17]精確控制橋梁斷面運動，采用具備非線性彈簧效果的多變形子區(qū)域動網(wǎng)格方法更新網(wǎng)格。

2.1 自激氣動力分析

橋梁斷面按縮尺比1∶40建模，豎向振動幅值選取區(qū)域為0.001B～0.14B(B=0.867 5 m)，扭轉(zhuǎn)振動幅值選取區(qū)域為0.1°～30°。選定一振動幅值，分別對斷面賦予單頻單一自由度強迫振動，每個振幅下的折算風速選取區(qū)間為2～10，來流風速取為20 m/s，計算出相應的自激氣動力，具體的工況設(shè)置如表2所示。

表2 工況設(shè)置Tab.2 Operation setting

通過對各折算風速狀態(tài)下的自激氣動力進行頻譜分析，可以研究氣動力隨幅值變化的氣動規(guī)律。由于算例橋主梁的顫振臨界風速約為66.3 m/s，振動頻率約為0.285 Hz，相應的折算風速約為6.7左右，限于篇幅，本文以折算風速為6的狀態(tài)為例，分析其氣動頻譜特性。

首先列出流線型箱梁斷面0°攻角時在單頻豎向振動下的部分氣動力頻譜分析圖(去除靜力分量)。

由圖2和圖3可以看出，當箱梁斷面做單頻單一豎向振動時，產(chǎn)生的自激力頻譜中除了包含與振動頻率同頻的基頻分量外，還會有些許頻率是斷面振動頻率整數(shù)倍的高次諧波分量。在氣動力成分占比中，基頻諧波成分明顯，盡管隨著振幅增大高次諧波占比從無到有，但是與基頻成分相比很小，幾乎可以忽略。

圖2 單頻豎向不同幅值振動下的氣動升力頻譜圖Fig.2 Aerodynamic lift spectrum under different vertical amplitude vibrations of the single frequency

圖3 單頻豎向不同幅值振動下的氣動扭矩頻譜圖Fig.3 Aerodynamic torsional spectrum under different vertical amplitude vibrations of the single frequency

接著列出流線型箱梁斷面0°攻角時在單頻扭轉(zhuǎn)振動下的部分氣動力頻譜分析圖。從圖4中可以看出，與單頻單一豎向振動相同，單頻單一扭轉(zhuǎn)振動下也會產(chǎn)生以基頻為主、整數(shù)倍頻摻雜的多頻氣動力。當振幅增大時，二次倍頻占比開始明顯，三次倍頻諧波及更高次諧波占比可以忽略。隨著幅值再次增大，更高次波成分逐漸顯露。

圖4 單頻扭轉(zhuǎn)不同幅值振動下的氣動升力頻譜圖Fig.4 Aerodynamic lift spectrum under different torsional amplitude vibrations of the single frequency

為分析高次諧波的成分，調(diào)整豎向坐標為對數(shù)坐標，給出升力幅值頻譜圖。如圖5所示當振幅從2°-5°-18°時，三次諧波分量愈發(fā)明顯，其占比逐漸超越了二次諧波，三次諧波的重要性在大振幅振動中顯現(xiàn)出來。這與王騎的試驗結(jié)果是吻合的。對于產(chǎn)生的升力或者力矩，隨振幅的變化規(guī)律是同步的。

圖5 單頻扭轉(zhuǎn)不同幅值振動下的氣動升力頻譜圖(對數(shù)坐標)Fig.5 Aerodynamic lift spectrum under different torsional amplitude vibrations of the single frequency

綜上分析，在單頻單自由度振動時，箱梁斷面會產(chǎn)生非線性氣動力。隨著振幅的增大，倍頻氣動力占比逐漸增大，其復雜的特性規(guī)律難以捕捉。同樣，氣動力隨幅值的變化規(guī)律再也不能用“線性關(guān)系”簡單描述， 相應的通過氣動力識別的顫振導數(shù)則不再隨振幅微幅變動。

2.2 非線性顫振導數(shù)

用最小二乘擬合[18]識別出隨折算風速與振幅變化的非線性顫振導數(shù)，通過插值得到的三維曲面圖，見圖6。

圖6 顫振導數(shù)曲面Fig.6 Flutter derivative surface

那么當實際風速跨過傳統(tǒng)線性理論預測的顫振發(fā)散風速時，振動幅值指數(shù)倍增大，因斷面的氣動參數(shù)的非線性變化，橋梁斷面的系統(tǒng)振動阻尼可能由負轉(zhuǎn)正，振動又趨于平穩(wěn)。

3 復模態(tài)特征值求解方法

為了簡化公式的形式且更直觀地說明，將自激力進行符號表示簡化

式中，Hi，Ai(i=1,2,3,4)是有量綱的顫振導數(shù)，它們是因子與顫振導數(shù)合并的結(jié)果，同樣其也僅和折算風速與振幅有關(guān)。代入式(1)和式(2)中

(3)

用矩陣記號，將式(3)改為

(4)

(5)

那么就轉(zhuǎn)變成著名的狀態(tài)方程，向量(x,y,h,α)T稱為狀態(tài)向量，因為它能唯一確定二自由度的運動狀態(tài)。式中每個子矩陣都是2×2矩陣，0表示零矩陣。

因為顫振臨界狀態(tài)時結(jié)構(gòu)處于等頻運動狀態(tài)，因此可以設(shè)

(6)

式中，ω應為實數(shù)。將式(6)代入式(5)得

(7)

注意對任意的eiωt，式(7)都是滿足的，于是應有

(8)

為了使式(8)有非零解，那么其系數(shù)行列式必須等于零。式(8)的系數(shù)行列式是一個關(guān)于ω的四階方程，稱為特征方程。該方程中有兩個未知數(shù)ω與U，如果不考慮顫振振幅因素的影響，不妨代入一組微振幅下的顫振導數(shù)，通過將試探的風速值代入系數(shù)矩陣中，并調(diào)用復系數(shù)矩陣特征值算法，得到ω的四個根，一旦某一個根是實數(shù)，那么此ω為顫振臨界狀態(tài)下結(jié)構(gòu)的振動頻率，對應的風速為臨界風速，如上所述，求解此特征方程的方法是一個逐步搜索的過程。

則系統(tǒng)在顫振臨界狀態(tài)之前的任意狀態(tài)、任意時刻的運動方程可表示為

(9)

當系統(tǒng)以某階復頻率做主振動時(r=1 or 2)，振動系統(tǒng)的周期表達式為

(10)

根據(jù)歐拉公式變換為

(11)

豎向與扭轉(zhuǎn)振動相位角

(12)

豎向與扭轉(zhuǎn)振動幅值比為

(13)

此處絕對值為模的概念。

值得注意的是，選取任意一組微幅振動下的顫振導數(shù)代入計算后，在達到了顫振臨界狀態(tài)時，會得到式(13)中幅值比一值。由于考慮了顫振振幅響應因素，應滿足代入的顫振導數(shù)對應的幅值比與臨界狀態(tài)的幅值比相同。故之前代入的一組顫振導數(shù)只能稱作是一種試探，因此不斷調(diào)整試探振幅使得其滿足上一結(jié)果的幅值比，如此往復，直到代入的顫振導數(shù)對應的幅值比與求解的幅值比相同為止，即為顫振的最終狀態(tài)。上述方法稱為考慮顫振振幅因素的二自由度復模態(tài)特征值求解方法。其求解流程圖如圖7所示。

圖7 改進的復模態(tài)特征值求解方法流程圖Fig.7 Improved flow chart of complex modal eigenvalue solution method

需要說明的是，由于無法預知顫振時起振的具體幅值，故而依舊只能選取微幅概念中的某一幅值，保持某一自由度振幅不變，不斷改變另一振幅來滿足彎扭幅值比的要求。本文研究對象為流線型箱梁斷面，并且其扭彎頻率比為2.13，很有可能發(fā)生顫振時為扭轉(zhuǎn)驅(qū)動機制。為研究方便，故而將扭轉(zhuǎn)振幅作為不變量。需要強調(diào)的是，對于求得的最終狀態(tài)對應的振幅，無論假設(shè)豎向或扭轉(zhuǎn)幅值為不變量，選取另一變量的任意振幅值作為初始試探值，最終狀態(tài)的結(jié)果都是趨于一致的，即初值的改變只會影響搜索路徑。

4 顫振后狀態(tài)探索研究

4.1 顫振臨界風速求解

假設(shè)顫振時扭轉(zhuǎn)振幅為0.1°，結(jié)合上述求得的顫振導數(shù)曲面圖，基于考慮顫振振幅因素的二自由度復模態(tài)特征值求解方法，得到顫振臨界風速為66.34 m/s。為了使本文的方法更具有說服力，選用斷面結(jié)構(gòu)類似的南京四橋主梁斷面進行驗證。將唐煜研究中運用同種方法求得的三維顫振導數(shù)代入計算，得到扭轉(zhuǎn)角1°微幅振動下的顫振臨界風速為70.4 m/s，這與文獻[19]中用直接風洞實驗法測得的風速70.8 m/s十分接近，除了反映顫振導數(shù)計算的合理性、基于風洞試驗識別參數(shù)的理論分析法適用性外，也充分說明了考慮顫振振幅因素的二自由度復模態(tài)特征值求解方法是準確的。

4.2 顫振后的振動狀態(tài)可能性假設(shè)

一旦來流風速超過臨界風速，系統(tǒng)的氣動扭轉(zhuǎn)阻尼為負，橋梁斷面作發(fā)散運動，振幅隨之增大，相應的氣動自激力的高階成分逐漸顯現(xiàn)，氣動參數(shù)的非線性特性愈發(fā)顯著。橋梁斷面周圍的氣動環(huán)境隨著振幅變化而改變，那么就極有可能出現(xiàn)再次穩(wěn)定的現(xiàn)象。如風洞試驗中的軟顫振現(xiàn)象，就可以從某方面說明氣動參數(shù)隨振幅的變化導致起振后斷面仍然陷入環(huán)振蕩現(xiàn)象。

合理的氣動力與振動方程搭接方法、合適的氣動參數(shù)識別方法都使得建立非線性顫振框架模型困難重重。為了在現(xiàn)有的顫振框架下進行顫振后狀態(tài)的有益探索，故而結(jié)合塔克馬橋風毀工程實際事件、軟硬顫振風洞實驗現(xiàn)象，假設(shè)橋梁顫振的最終形態(tài)都為同一頻率的簡諧扭轉(zhuǎn)和簡諧豎向的耦合運動。那么在該假設(shè)下，對顫振后橋梁斷面振動形態(tài)進行穩(wěn)態(tài)響應設(shè)想：當風速達到顫振臨界風速時，橋梁斷面作簡諧平穩(wěn)運動；一旦越過臨界風速，系統(tǒng)阻尼由零轉(zhuǎn)負，橋梁斷面做發(fā)散運動；隨著振動幅值增大，橋梁斷面氣動環(huán)境發(fā)生改變，在某一振動幅值下，系統(tǒng)阻尼由負轉(zhuǎn)正之際，橋梁振動形態(tài)又趨于穩(wěn)定；當風速繼續(xù)增長到某一值，系統(tǒng)阻尼由正轉(zhuǎn)負，在上一穩(wěn)定振幅下又會出現(xiàn)簡諧振動穩(wěn)定狀態(tài)；再次越過上一穩(wěn)定風速，斷面又將做發(fā)散運動；當斷面幅值不斷增長，系統(tǒng)阻尼再也不能由正轉(zhuǎn)負時，橋梁斷面也就一直發(fā)散振動直至破壞。

盡管上述假設(shè)將顫振發(fā)生過程中的非線性行為描述弱化為簡諧環(huán)振動現(xiàn)象，但是簡諧環(huán)振動假設(shè)可以幫助我們認知非線性氣動彈性響應的一些規(guī)律，挖掘出主要相關(guān)的氣動因素。

4.3 顫振幅值響應特性分析

根據(jù)上述設(shè)想，研究顫振后的振動狀態(tài)，基于考慮顫振幅值因素的二自由度復模態(tài)特征值求解方法，給出非線性顫振幅值響應搜索程序流程圖，如圖8所示。

圖8 非線性顫振幅值響應搜索程序流程圖Fig.8 Flow chart of the search program for nonlinear flutter amplitude response

在顫振幅值響應風速搜索過程中分為整體風速與內(nèi)部風速，設(shè)置內(nèi)部風速的目的是便于程序整體編譯，整個內(nèi)部風速搜索可以視為一個子程序，整個風速搜索過程需遵循以下步驟：

步驟1設(shè)定基本參量如初始風速、初始頻率、起振微小扭轉(zhuǎn)與豎向振幅和最大限定扭轉(zhuǎn)角等，選定整體風速初始值；

步驟2進入內(nèi)部風速搜索，微幅等距增長風速值，迭代計算指定振幅下的系統(tǒng)頻率和阻尼，每階風速后都需判定該階下的系統(tǒng)阻尼，若系統(tǒng)阻尼由正轉(zhuǎn)負且彎扭比合適，則停止內(nèi)部風速搜索，每一個完整的子程序搜索過程的最終風速作為本次的內(nèi)部風速值，記錄相應參數(shù)值；

步驟3轉(zhuǎn)入整體風速搜索，若內(nèi)部風速搜索值小于或等于整體風速值，則認為顫振發(fā)生，并記錄當前風速，否則增長整體風速，轉(zhuǎn)入步驟2，若整體風速達到限定最大風速值則停止搜索；

步驟4一旦顫振發(fā)生，保持整體搜索風速不變，由于此時的系統(tǒng)阻尼為負，斷面作發(fā)散振動，等距微幅增長扭轉(zhuǎn)角，并結(jié)合上一扭轉(zhuǎn)角于顫振時的彎扭比計算出當下的豎向振幅試探值，轉(zhuǎn)入內(nèi)部搜索子程序中，記錄相關(guān)參數(shù)，轉(zhuǎn)入整體風速搜索，若扭轉(zhuǎn)角大于限定的扭轉(zhuǎn)角上限值則停止搜索；

步驟5當內(nèi)部風速搜索值小于或等于整體風速值時，則繼續(xù)增大振幅轉(zhuǎn)入步驟4計算；

步驟6當內(nèi)部風速搜索值大于整體風速值時，記錄此時的振幅參數(shù)，保持振幅不變，增長整體風速，進入內(nèi)部搜索程序，當整體風速不小于內(nèi)部風速時，轉(zhuǎn)步驟3。

需要注意的是，對一個參數(shù)的計算和迭代求解都是在指定的振幅條件下進行的。本文中增加了彎扭幅值比的概念，為了便于程序的整體編譯，設(shè)置了內(nèi)部風速搜索。內(nèi)部風速僅具有比較作用，當內(nèi)部風速大于整體風速時，其意義為在當前整體風速下，橋梁斷面于指定幅值條件振動下系統(tǒng)阻尼為正；當內(nèi)部風速小于整體風速時，橋梁斷面的系統(tǒng)阻尼為負。故而運行上述程序進行求解，結(jié)果如圖9所示，其前段以紅色圈標注的細部圖如圖10所示。

圖9 顫振振幅響應總體圖Fig.9 Whole graph of flutter amplitude response

從圖10中可以看出，設(shè)定的起振扭轉(zhuǎn)角為0.1°，當越過顫振臨界風速66.34 m/s后，橋梁斷面的振幅響應并不是直接發(fā)散的，而是隨著風速的增長爬坡式的增大，扭轉(zhuǎn)角從0.1°，-0.6°，-0.9°，-1.5°，-1.8°，-2.1°，-2.4°，-2.7°緩緩增大。當越過風速66.42 m/s時，橋梁斷面振幅響應出現(xiàn)了階躍性變化，扭轉(zhuǎn)角從2.7°陡增到19.6°，兩個極限穩(wěn)態(tài)間的風速間隔約0.1 m/s，一旦越過風速66.49 m/s，如圖9所示，顫振振幅響應又陷入了爬坡式的增長，直到風速越過75.61 m/s后，隨著振幅的增長，斷面再也達不到穩(wěn)態(tài)響應。

圖10 顫振振幅響應局部圖Fig.10 Local graph of flutter amplitude response

圖11 顫振響應彎扭幅值比Fig.11 Crankle amplitude ratio of flutter response

圖12 顫振響應相位角Fig.12 Phase angle of flutter response

該橋型斷面與大多流線型斷面類似，都是以扭轉(zhuǎn)形態(tài)為主的顫振機制，在顫振振幅響應搜索中，如圖11、圖12所示，起始的彎扭振幅比幾乎區(qū)別不大，當越過風速66.49 m/s時，彎扭比陡增至14.2，之后相位角呈現(xiàn)先增加后減小的趨勢，扭轉(zhuǎn)運動的幅度比例增大，扭轉(zhuǎn)與豎向運動的相位角不斷縮小，橋梁斷面的運動形態(tài)逐漸向單一自由度逼近。

5 結(jié) 論

本文基于風洞試驗識別參數(shù)的理論分析法，計算了流線型斷面的三維顫振導數(shù)曲面、改進了傳統(tǒng)的復模態(tài)特征值求解流程和建立了非線性顫振振幅響應搜索程序。

(1)通過控制振動幅值識別出不同振幅下的氣動力，單頻單自由度振動下橋梁斷面氣動力會有非線性成分，隨著幅值增大，非線性成份比重擴大，單頻單一扭轉(zhuǎn)運動下尤其明顯。值得注意的是在扭轉(zhuǎn)運動下，當振幅達到某一值時，三次諧波力占比超越了二次諧波力。

(2)通過最小二乘擬合識別出隨振幅與折算風速變化的顫振導數(shù)，并通過插值繪制成三維曲面。對于流線型箱梁斷面而言，顫振導數(shù)H1,H4,A1,A4變化不大，而A2,H3則變化顯著，一定程度上反映了氣動力隨振幅改變的非線性性質(zhì)。

(3)發(fā)展傳統(tǒng)的顫振理論，引入彎扭幅值比概念，擴展出考慮顫振振幅因素的二自由度復模態(tài)特征值求解方法，并通過實橋算例證明了該方法的準確性與適用性。

(4)編譯非線性顫振振幅響應搜索程序，給出了微幅顫振后的狀態(tài)。流線型箱梁發(fā)生顫振后，振幅會出現(xiàn)階躍性增長，達到某一大幅值后，又基本處于振幅緩增的振動狀態(tài)。

本文從顫振響應振幅入手，試圖揭示顫振后的非線性動力學行為，對顫振后狀態(tài)進行了嘗試性的描述與分析。