典型橋梁斷面風致振動的氣動能量特征分析

劉祖軍，楊詠昕，周冰凌

（1.華北水利水電大學 土木與交通學院，河南 鄭州 450045；2.同濟大學 橋梁工程系，上海 200092；3.四川蜀泰建筑設計有限公司，四川 成都 610000）

0 引 言

顫振屬于自激振動，其物理機理可以從能量的角度加以深入的闡釋。處于氣流中的橋梁其能量反饋機制表現(xiàn)為氣流輸入到結構—氣流系統(tǒng)中的能量與結構阻尼耗散能量之間的平衡關系，當輸入到結構—氣流系統(tǒng)中的能量小于結構阻尼耗能時時，結構在初始擾動下將作衰減（阻尼）振動；而當輸入的能量大于結構阻尼耗能時結構在初始擾動下將作發(fā)散振動；兩者相等時結構在初始擾動下將作等幅簡諧振動。

在橋梁自激振動能量分析方面各國學者也進行了嘗試，Scanlan［1］最早建立了橋梁顫振的多模態(tài)分析方法，并從能量的觀點對橋梁的顫振穩(wěn)定性進行了很有價值的研究，給出了在一個振動周期內氣流沿橋梁斷面每延米輸入的總能量和結構耗能的表達式，但他僅給出了一個理論框架，如何從能量的角度對橋梁進行顫振分析，沒有給出具體的方法。

Larsen以CFD方法為基礎，根據(jù)離散渦計算中渦旋的運動規(guī)律提出了一個簡化分析模型［2］。這個模型描述了在橋梁斷面扭轉運動的一個周期里渦旋的運動情況，并通過積分估算由渦旋產(chǎn)生的氣動力對橋梁斷面所做的總功，通過用能量方法分析渦激力做攻與結構穩(wěn)定之間的關系。Larsen的研究具有很好的開創(chuàng)性，但是他的計算模型假定在旋渦沿橫截面移動時，旋渦升力保持不變并且結構阻尼為零，這和實際情況有很大的不同。

我國的劉高博士從結構與氣流系統(tǒng)內部能量平衡的觀點對系統(tǒng)的顫振進行研究［3］，發(fā)展了一種全橋多模態(tài)顫振分析方法－能量法，通過建立系統(tǒng)等效阻尼比與系統(tǒng)能量變化率之間的關系，推演了系統(tǒng)以及各階模態(tài)等效阻尼比的計算方法，并根據(jù)不同風速下系統(tǒng)能量變化率來判斷系統(tǒng)的顫振穩(wěn)定性。但這種方法不能很好的分析氣流與結構之間的細觀作用，特別是結構自激振動中系統(tǒng)能量的分布和轉化情況。同濟大學博士楊詠昕、劉祖軍基于分步分析的方法從能量角度對典型橋梁斷面顫振的機理進行了研究［4］。

在航空方面機翼的氣動彈性力學研究中，F(xiàn)razer早在1939年就分析了維持機翼強迫振動所需要的能量條件［5］，Nissim在20世紀70年代提出了氣動能量概念作為對機翼主動控制的理論基礎［6］，Jones分析了顫振發(fā)生的能量特征，上述工作構成機翼的氣動彈性力學研究應用的基礎［7］。在機械領域的葉輪機氣動彈性力學研究中Carta提出了用于顫振失穩(wěn)預測的能量穩(wěn)定判據(jù)［8］。后來Klose和Heinig通過研究證明能量法用于顫振失穩(wěn)的預測可被認為是特征值法在葉片高質量比下的特殊應用，并指出應用能量法預測的一些失效情況［9］。

本文通過風洞試驗研究了單箱的顫振性能，基于流固松耦合的計算策略，利用現(xiàn)有流體軟件的用戶自定義（UDF）功能，應用CFD數(shù)值方法，模擬了單箱顫振過程，根據(jù)分塊分析的思想研究了顫振過程中振動模型表面不同區(qū)域吸收氣流能量的特點，并分析了旋渦的非定常演化對模型表面不同區(qū)域壓力特性的影響。

1 單箱風洞測振試驗及結果分析

風洞試驗在同濟大學土木工程防災國家重點實驗室TJ－4邊界層風洞中進行，該風洞為低速回流式風洞，測振的試驗段尺寸為：寬0.814m，高0.8m，長2.0m，設計最大試驗風速為30m／s。單箱斷面依據(jù)大海帶東橋的主橋斷面為原型進行縮尺，其截面尺寸如圖1所示，縱向長度0.8m，采用鋁質芯棒，外覆泡沫制成。

圖1 單箱模型 （單位：mm）Fig.1 The model of box（unit：mm）

模型的基本參數(shù)如下：m＝1.3625kg／m，Im＝0.01277kg·m，豎向頻率fh＝2.825Hz，扭轉頻率fa＝7.451Hz。測振時采用激光位移計記錄模型試驗風速下的位移信號。

當風速小于18m／s時結構的豎向位移振幅很小，但是存在較小幅度的扭轉振動。當風速超過18m／s后結構振動位移極大值和方差就突然增大，豎向振動的參與程度也不斷增加，直到風速到達20.8m／s時結構出現(xiàn)了顫振失穩(wěn)現(xiàn)象。圖2和圖3給出了顫振臨界風速下的位移時程曲線。

2 顫振過程的數(shù)值模擬

2.1 單箱顫振數(shù)值模擬方法

圖2 豎向振動位移時程響應 （20.8m／s）Fig.2 Vertical vibration displacement vs.time curve of box section（20.8m／s）

圖3 扭轉振動位移時程響應（20.8m／s）Fig.3 Torsional vibration displacement vs.time curve of box section（20.8m／s）

橋梁顫振發(fā)散是流體與結構相互動力作用的結果。處理流固耦合的問題可以采用強耦合和松耦合這兩種方法。強耦合求解方法是以求解流體運動的迭代過程為主，將物體運動的求解耦合在流體運動求解過程中.松耦合方法在每一個時間步內分別對流體和結構依次求解，然后通過中間平臺交換數(shù)據(jù)信息，從而實現(xiàn)兩個場的耦合求解。這種方法對計算機硬件的要求不高，并且可以充分發(fā)揮CFD（Computational Fluid Dynamics）計算和CSD（Computational Structural Dynamics）計算的各自優(yōu)點，能實現(xiàn)數(shù)值計算的模塊化。在分析橋梁結構風致振動問題時多采用松耦合的計算策略。

本文采用松耦合的計算方法，利用現(xiàn)有商業(yè)軟件Fluent提供的基于ALE（Arbitrary Lagrange Euler）方法的動網(wǎng)格技術，來實現(xiàn)顫振的數(shù)值模擬，通過軟件的用戶自定義函數(shù)描述模型及周圍氣流的剛體運動，具體的計算流程圖如圖4所示。

圖4 顫振數(shù)值模擬流程圖Fig.4 The process of flutter numerical simulation

顫振數(shù)值模擬過程的具體步驟為：

（1）首先進行流體域內固定斷面的非定常繞流計算，將計算獲得的三分力系數(shù)通過Fluent軟件的用戶自定義功能輸出到指定的文本中，然后再將其換算為作用在模型表面上的升力和升力矩，從而獲取t時刻的氣動力。

（2）根據(jù)計算所獲得的升力和升力矩，通過數(shù)值方法對振動結構進行結 構動力學求解。根據(jù)newmark－β方法求解豎向和扭轉振動方程，獲得t＋Δt時刻結構運動的豎向速度和扭轉速度。

（3）將模型振動的速度賦值給包圍模型斷面并隨模型一起做剛體運動的部分流體區(qū)域，通過Fluent提供的用法自定以函數(shù)（UDF）描述模型及周圍流體的運動，采用動網(wǎng)格方法進行流體域的求解，獲得振動斷面的氣動升力和升力矩。

（4）重復（2）－（3）步，獲得橋梁斷面振動響應的時程。

（5）根據(jù)斷面的響應時程曲線判斷振動是否已經(jīng)發(fā)散，如果已經(jīng)發(fā)散則該風速即為顫振臨界風速，否則增加風速按照（1）－（4）步重新計算直到達到振動發(fā)散。

2.2 顫振數(shù)值模擬結果分析

數(shù)值模擬時采用商業(yè)軟件Fluent提供的RANS方法的k－ωSST兩方程模型［10－11］，計算域的大小參考了同濟大學土木工程防災國家重點試驗室TJ－4風洞的中段試驗端設置，計算域沿流線長度為3m（其中上游1m，下游2m），橫向寬度為0.8m。計算時壁面附近最小網(wǎng)格尺度為0.0004m，計算域采用分塊結構化網(wǎng)格。網(wǎng)格數(shù)量為15.6萬。

計算參數(shù)設置為：動量，湍動能和能量耗散均采用兩階迎風格式進行離散，壓力速度耦合采用SIMPLE算法，求解器采用分離式，計算模式選用兩階隱式。邊界條件的設定為：速度入口，入口初始風速15m／s，湍流強度0.5%，壓力出口，計算域的上端和下端設為對稱邊界條件，表面采用無滑移的壁面條件。

由于RANS方法平均的結果是將流場時空變化的細節(jié)一概抹平，丟失了流場脈動運動中的大量信息，并且為了使N－S方程封閉，提出了Reynolds應力模型和渦粘模型12－14］，這些模型存在著對經(jīng)驗數(shù)據(jù)過分依賴和預報精度相對較低的局限性。

由于受到計算機硬件條件的限制，本文仍然采用了RANS方法模擬箱梁顫振過程，并提取了模型表面壓力進行能量分析，需要說明的是計算獲得壓力是平均后的結果，因此本文計算的能量可以反映結構振動過程中的總體變化規(guī)律，但忽略了瞬時脈動壓力對結構振動能量的貢獻。

圖5 豎向振動位移時程響應 （20.8m／s，CFD數(shù)值模擬）Fig.5 Vertical vibration displacement vs time curve of box section（20.8m／s，CFD）

該箱梁的顫振臨界風速為20.8m／s。圖5、圖6是平板顫振時刻的位移時程，數(shù)值模擬所獲的顫振頻率為4.92Hz，實測顫振頻率為4.59Hz。

圖6 扭轉振動位移時程響應（20.8m／s，CFD數(shù)值模擬）Fig.6 Torsional vibration displacement vs time curve of box section（20.8m／s，CFD）

3 基于分塊分析思路的顫振過程能量分析

3.1 分塊分析能量計算方法

分塊分析主要是為了比較詳細的捕捉表面壓力的變化，尋找引起顫振發(fā)散的主要的能量輸入特點，同時能夠考慮尾部旋渦對氣流能量輸入方式的影響和影響范圍。根據(jù)對斷面進行分區(qū)，然后按照下式分別計算不同區(qū)域的能量輸入特點.單位區(qū)域的能量計算公式如下：

單位面積上一個周期內氣動力輸入到系統(tǒng)的能量w：

則面積為S（圖7的網(wǎng)格線部分）的區(qū)域氣動力輸入的能量w（s，t）：

P（x，t）為結構表面的壓力，v（x，t）為結構運動速度，n（x，t）為壓力與速度之間的夾角。

圖7 分塊分析示意圖Fig.7 The Block Analysis

根據(jù)本文提出的分塊分析思路，編制了相應的計算程序，通過提取數(shù)值計算獲得的模型表面壓力，詳細的研究了在一個周期內氣流向結構輸送能量的過程，模型表面進行了如圖8所示的分區(qū)。

圖8 模型表面分區(qū)Fig.8 The partition of the model

圖9～圖14是模型扭轉振動時不同分區(qū)輸入到系統(tǒng)的能量。第1個區(qū)域在前半個周期內上表面消耗能量，下表面吸入能量，后半個周期則反之。一個完整周期內輸入到系統(tǒng)的能量為正，但數(shù)值較小。上下表面輸入到系統(tǒng)能量具有較好的周期性特征，具體如圖9所示。第4個區(qū)域的輸入系統(tǒng)的能量特點和1區(qū)基本相似，能量的輸入也具有周期性的特點，二者的不同之處在于第4個區(qū)域是消耗系統(tǒng)的能量。第4區(qū)域受到尾部旋渦的直接作用。而尾部旋渦運動具有一定的規(guī)律性，并且尾端風嘴處上下側交替出現(xiàn)的旋渦聚集了較大的能量，控制了模型尾端的自由振動，使得模型尾端的運動和尾部旋渦的運動趨勢基本同步，因此造成了該區(qū)域氣流輸入能量具有明顯的周期性（圖15）。

圖9 1區(qū)扭轉運動的能量隨時間變化關系Fig.9 The energy of torsion in the first partition vs.time curve

圖10 2區(qū)扭轉運動的能量隨時間變化關系Fig.10 The energy of torsion in the second partition vs.time curve

圖11 3區(qū)扭轉運動的能量隨時間變化關系Fig.11 The energy of torsion in the third partition vs.time curve

圖12 4區(qū)扭轉運動的能量隨時間變化關系Fig.12 The energy of torsion in the fourth partition vs.time curve

圖13 各分區(qū)扭轉振動能量隨時間變化關系Fig.13 The energy of torsion in every partition vs.time curve

第2個區(qū)域為系統(tǒng)提供了較多的能量，在整個振動過程中上表面在前1／4周期消耗能量，后3／4周期內吸入能量，下表面的情況與之相反，總能量為正，數(shù)值較大，且增加較快，如圖10。第3個區(qū)域也就是風嘴部位是扭轉系統(tǒng)的主要能量來源，在一個完整周期內上下表面均吸入系統(tǒng)能量并且數(shù)值較大（圖11）。因此在顫振過程中模型的主要吸能部位分布在迎風端風嘴及其附近的區(qū)域。從第23個區(qū)域的能量隨時間的變化關系曲線上可以看出這兩個部位的能量輸入的周期性不明顯，受尾部旋渦運動的影響較小。

各分區(qū)對扭轉振動的能量貢獻如圖13所示，從單箱扭轉振動主要能量輸入曲線可以發(fā)現(xiàn)氣流輸入到單箱的主要能量在一個周期內是不斷增加的，所以單箱的顫振多為突然性失穩(wěn)。豎向振動的主要能量輸入則是主要由第2個區(qū)域引起，具體如圖14。上述分析表明單箱顫振時扭轉能量的主要吸入部分是迎風端的風嘴和靠近風嘴附近的區(qū)域如圖15所示。

圖14 各分區(qū)豎向振動能量隨時間變化關系Fig.14 The energy of vertical in every partition vs.time curve

圖15 氣流能量輸入圖Fig.15 The energy input by air

3.2 顫振過程旋渦演化對模型不同區(qū)域壓力特性的影響

通過CFD數(shù)值方法模擬了單箱處于顫振臨界狀態(tài)的繞流特征。由于振動模型周圍的流場是千變萬花的，即使模型振動到同一位置時，其周圍對應的流場也完全不同，因此采用了相位平均的方法［15］探討了單箱尾部風嘴處的旋渦的演化規(guī)律。

從流場的相位平均的分析結果（圖16）可以看出箱梁振動過程中尾部風嘴處的旋渦存在這樣的演化過程：當結構處于平衡位置時，在風嘴下側存在回流區(qū)，將要形成旋渦；當結構由平衡位置振動到負最大位移時在風嘴的上下側同時存在旋渦，但是上側旋渦尺度較小，下部旋渦尺度很大，結構處于下旋渦的控制下；隨后結構負最大位移回復到平衡位置時，上側旋渦消失，下部旋渦的尺度減小，由平衡位置振動到正最大位移時，在在風嘴的上下側又產(chǎn)生了一對尺度較大，接近橢圓形的旋渦，結構處于上下旋渦的共同控制下。

圖16 單箱斷面的旋渦驅動結構運動的流（顫振風速20.8m／s，數(shù)值模擬）Fig.16 Vortex driven structural movement of box section（20.8m／s，CFD ）

圖17 箱梁表面4個分區(qū)上下表面壓力隨時間變化關系Fig.17 Four regional of box girder surface pressure versus time

從4個區(qū)域輸入到系統(tǒng)扭轉振動的能量隨時間的變化關系，我們可以推測出尾部旋渦對氣流輸入到系統(tǒng)能量的影響。模型尾端風嘴附近的區(qū)域受到旋渦的直接作用，輸入到系統(tǒng)的能量具有很明顯的周期性特點。圖17給出了4個區(qū)域上下表面壓力分布隨時間的變化關系。為了清楚說明尾部旋渦對模型不同部位的影響定義了上下表面壓力差的相對變化量：＝（Cpu－Cpd）／Cpu，其中Cpu、Cpd為模型上下表面最大壓力。第4個區(qū)域的（4）＝4.28（1）＝1.43（2）＝1.31，（3）＝0（圖18），因此可以看出越靠近尾流區(qū)域，上下表面壓力差值的相對值越大，第4區(qū)域最大，而距尾流區(qū)最遠的3區(qū)上下表面之間的壓力相對差很小，兩個曲線基本重合。從流場機制上說箱梁尾部風嘴附近受到旋渦的直接推動作用，其表面壓力的變化特點與旋渦的運動規(guī)律很接近，具有明顯的周期性。尾部風嘴上下側的旋渦交替控制結構的運動，因此造成尾部風嘴區(qū)域的上下表面的壓力相對波動較大，而遠離尾流的區(qū)域受此影響較小，其上下表面的壓力波動的程度相對較小。

圖18 箱梁表面4個分區(qū)上下表面相對壓力差的波動系數(shù)Fig.18 The volatility coefficient of relative pressure difference on the four model surface regional

4 結 論

本文通過風洞試驗測試了單箱的顫振性能，采用CFD數(shù)值計算方法，并基于流固松耦合的計算策略模擬了單箱的顫振過程，結合分塊分析的思路定量研究氣流對振動模型表面不同區(qū)域的能量輸入特點及其對模型不同部位表面壓力的影響。通過上述研究得出以下幾點結論：

（1）分塊分析的計算結果表明迎風端的風嘴及其附近區(qū)是扭轉振動能量的主要輸入部位，對系統(tǒng)的振動的穩(wěn)定性起重要作用。

（2）從模型表面不同區(qū)域的能量輸入特點可以發(fā)現(xiàn)由于模型尾部風嘴附近上下側的旋渦交替作用具有明顯的規(guī)律性，因此造成模型靠近尾流區(qū)域的氣流能量輸入方式具有較明顯的周期性，而模型迎風端的氣流輸入規(guī)律性較差，周期性不明顯。

（3）分塊分析結果表明由于箱梁尾部風嘴上下側的旋渦交替作用對結構振動的影響很大，從而造成了尾部風嘴區(qū)域的上下表面壓力相對波動較大。而遠離尾部旋渦的模型區(qū)域受旋渦交替作用影響較小，因此使得該模型區(qū)域上下表面的壓力波動的程度相對較小。

（4）CFD數(shù)值模擬的結果表明單箱顫振過程中模型的周期性振動與模型尾部風嘴附近的上下側旋渦的交替作用存在著一定的聯(lián)系。

［1］SCANLAN R H.The action of flexible bridges under wind，I：flutter theory［J］.J of Sound and Vibration，1978，60（2）：187－199.

［2］LARSEN A.Aerodynamics of the Tacoma narrows bridgc－60years later［J］.Joumd of Structural Engineering International，2000，（10）：243.

［3］LIU Gao，WANG Xiuwei，QIANG Shizhong，et al.Flutter analysis of long－span suspension bridges by energy method［J］.China Journal of Highway and Transport，2000，12（3）：20－24.（in Chinese）劉高，王秀偉，強士中，等.大跨度懸索橋顫振分析的能量法［J］.中國公路學報，2000，12（3）：20－24.

［4］LIU Zujun，GE Yaojun，YANG Yongxin.Energy transformation mechanism of coupled bending－torsional flutter［J］.Jouranl of Tongji University（Natural Science），2011，39（7）：949－954.（in Chinese）劉祖高，葛耀君，楊詠昕.彎扭耦合顫振過程中的能量轉換機理［J］.同濟大學學報（自然科學版），2011，39（7）：949－954.

［5］FRAZER R A.On the power input required to maintain forced oscillations of an aeroplane wing in flight［J］.ARCR＆M，1939.

［6］NISSIM E.Flutter suppression using activecontrols based on the concept of aerodynamic energy［R］.NASATND 6199，1997.

［7］JONES J G.On the energy characteristics of the aerodynamic matrix and the relationship to possible flutter［J］.The Aeronautical Quqrterly，1983，3：212－225.

［8］CARTA F O.Coupled blade－disk－shroud flutter instabilities in turbojet engine rotors［J］.Journal of Engineering for Power，1967：419－426.

［9］KLOSE A，HEINIG K.A comparison of flutter calculations based on Eigenvalue and energy method［R］.AD－A2119741，1989.

［10］LUBCKE H，SCHMIDT S，RUNG T，et al.Comparison of les and rans in bluff－body flows［J］.Journal of Wind Engineering and Industrial Aerodynamics，2001，89（14－15）：1474－1485.

［11］TETSURO TAMURA，YOSHIYUKI ONO.LES analysis on aeroelastic instability of prisms in turbulent flow［J］.Journal of WInd Engineering and Industrial Aerodynamics，2003，（91）：51－64.

［12］LIU Zujun，GE Yaojun，YANG Yongxin.Identification of aerodynamic parameters of a flat plate with multimoving grid technique based on large eddy simulation［J］.Journal of Vibration and Shock，2011，30（4）：156－160.（in Chinese）劉祖軍，葛耀君，楊詠昕.基于大渦模擬方法的多層動網(wǎng)格技術識別平板氣動參數(shù)［J］.振動與沖擊，2011，30（4）：156－160.

［13］ZHOU Zhiyong，YANG Likun.Numerical study on the mechanism of torsional flutter forΠ－shaped section［J］.Acta Aerodynamics Sinica，2009，27（6）：683－689.（in Chinese周志勇，楊立坤.Π形板梁分離流扭轉顫振機理數(shù)值研究［J］.空氣動力學學報，2009，27（6）：683－689.

［14］ZHOU Zhiyong.Numerical simulation study on Reynolds number effect on bridge decks［J］.Acta Aerodynamics Sinica，2009，27（6）：664－670.（in Chinese）周志勇.橋梁斷面雷諾數(shù)效應數(shù)值模擬研究［J］.空氣動力學學報，2009，27（6）：664－670.

［15］FUJISAWA N，TAKEDA G，IKE N.Phase－averaged characteristics of flow around a circular cylinder under acoustic excitation control［J］.Journal of Fluids and Structures，2004，19：159－170.