流線型箱梁斷面的非線性顫振幅值特性研究

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

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

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

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

基于風洞試驗研究，Scanlan等將顫振導數表示的非定常氣動自激力引入振動平衡方程中，摒棄了流固耦合的直接求解，建立了經典的顫振框架模型，采用半逆解法進行顫振求解。Wilde等[7-8]采用二自由度的狀態空間復模態特征值求解方法(Complex Eigenvalue Analysis, CEVA)，其基本思想是將二維顫振問題轉變為求解復特征值問題，以某階模態下的阻尼比為零時，作為系統發散的依據，對應的風速為顫振臨界風速。基于該框架下的解具有明確的數學意義，當風速超過這一臨界點橋梁斷面的振幅將指數級增加而出現發散，可實際現象中卻并非如此。朱樂東等[9]在風洞試驗中發現，一些鈍體斷面顫振發生后并不會出現如線性理論預測的那樣，振幅按指數級迅速增長，而是由于自激力的非線性效應穩定在一定幅值狀態。學者們開始關注橋梁進入顫振后的真實運動狀態，顫振研究逐步走向精細化。

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

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

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

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

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

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

(1)

(2)

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

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

表1 橋梁結構基本參數Tab.1 Basic parameters of the bridge structure

2 氣動參數識別

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

2.1 自激氣動力分析

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

表2 工況設置Tab.2 Operation setting

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

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

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

圖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°攻角時在單頻扭轉振動下的部分氣動力頻譜分析圖。從圖4中可以看出，與單頻單一豎向振動相同，單頻單一扭轉振動下也會產生以基頻為主、整數倍頻摻雜的多頻氣動力。當振幅增大時，二次倍頻占比開始明顯，三次倍頻諧波及更高次諧波占比可以忽略。隨著幅值再次增大，更高次波成分逐漸顯露。

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

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

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

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

2.2 非線性顫振導數

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

圖6 顫振導數曲面Fig.6 Flutter derivative surface

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

3 復模態特征值求解方法

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

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

(3)

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

(4)

(5)

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

因為顫振臨界狀態時結構處于等頻運動狀態，因此可以設

(6)

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

(7)

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

(8)

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

則系統在顫振臨界狀態之前的任意狀態、任意時刻的運動方程可表示為

(9)

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

(10)

根據歐拉公式變換為

(11)

豎向與扭轉振動相位角

(12)

豎向與扭轉振動幅值比為

(13)

此處絕對值為模的概念。

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

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

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

4 顫振后狀態探索研究

4.1 顫振臨界風速求解

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

4.2 顫振后的振動狀態可能性假設

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

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

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

4.3 顫振幅值響應特性分析

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

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

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

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

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

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

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

步驟5當內部風速搜索值小于或等于整體風速值時，則繼續增大振幅轉入步驟4計算；

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

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

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

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

圖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

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

5 結 論

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

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

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

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

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

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