賈明曉，劉祖軍，，楊詠昕

(1.華北水利水電學(xué)院 土木與交通學(xué)院，鄭州 450011;2.同濟大學(xué) 橋梁工程系，上海 200092)

數(shù)值方法與計算機技術(shù)進步提升了對顫振本質(zhì)的認(rèn)識，已能較準(zhǔn)確預(yù)測橋梁顫振的臨界風(fēng)速、顫振頻率及顫振形態(tài)。顫振屬自激振動，其物理機理可從能量角度進行闡釋。處于氣流中橋梁的能量反饋機制表現(xiàn)為氣流輸入結(jié)構(gòu)—氣流系統(tǒng)中能量與結(jié)構(gòu)阻尼耗散能量之間的平衡關(guān)系。當(dāng)輸入結(jié)構(gòu)—氣流系統(tǒng)中能量小于結(jié)構(gòu)阻尼耗能時，結(jié)構(gòu)在初始擾動下將作衰減(阻尼)振動;而當(dāng)輸入能量大于結(jié)構(gòu)阻尼耗能時結(jié)構(gòu)在初始擾動下將作發(fā)散振動;兩者相等時結(jié)構(gòu)在初始擾動下將作等幅簡諧振動。

對橋梁自激振動能量分析，Scanlan［1］最早提出橋梁顫振的多模態(tài)分析方法，并從能量觀點對橋梁顫振穩(wěn)定性進行研究，給出一個振動周期內(nèi)氣流沿橋梁斷面每延米輸入的總能量與結(jié)構(gòu)耗能表達式，但此僅為理論框架，如何從能量角度對橋梁進行顫振分析，并未給出具體方法。Larsen［2］以CFD方法為基礎(chǔ)，據(jù)離散渦計算中渦旋運動規(guī)律提出簡化分析模型。該模型描述橋梁斷面扭轉(zhuǎn)運動一個周期內(nèi)渦旋的運動情況，并通過積分估算由渦旋產(chǎn)生的氣動力對橋梁斷面所做總功，用能量方法分析渦激力做功與結(jié)構(gòu)穩(wěn)定之關(guān)系。文獻［2］頗具開創(chuàng)性，但其計算模型假定在旋渦沿橫截面移動時，旋渦升力保持不變且結(jié)構(gòu)阻尼為零，此與實際情況有很大不同。劉高［3］從結(jié)構(gòu)與氣流系統(tǒng)內(nèi)部能量平衡觀點對系統(tǒng)顫振進行研究，發(fā)展一種全橋多模態(tài)顫振分析方法-能量法，通過建立系統(tǒng)等效阻尼比與系統(tǒng)能量變化率之間的關(guān)系，推導(dǎo)出系統(tǒng)及各階模態(tài)等效阻尼比計算方法，并據(jù)不同風(fēng)速下系統(tǒng)能量變化率判斷系統(tǒng)顫振穩(wěn)定性。但該方法不能分析氣流與結(jié)構(gòu)間的細觀作用，尤其結(jié)構(gòu)自激振動中系統(tǒng)能量分布及轉(zhuǎn)化情況。

在機翼氣動彈性力學(xué)研究中，F(xiàn)razer［4］分析了維持機翼強迫振動所需的能量條件;Nissim［5］提出氣動能量概念作為對機翼主動控制的理論基礎(chǔ);Jones［6］分析顫振發(fā)生的能量特征;Carta［7］提出用于顫振失穩(wěn)預(yù)測的能量穩(wěn)定判據(jù);Klose等［8］通過研究證明能量法用于顫振失穩(wěn)的預(yù)測被認(rèn)為特征值法在葉片高質(zhì)量比下的特殊應(yīng)用，并指出用能量法預(yù)測的失效情況。

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

1 平板斷面風(fēng)洞測振試驗

1.1 平板風(fēng)洞測振試驗及結(jié)果分析

風(fēng)洞試驗在同濟大學(xué)土木工程防災(zāi)國家重點實驗室TJ-4邊界層風(fēng)洞中進行，該風(fēng)洞為低速回流式風(fēng)洞，測振試驗段尺寸為:寬 0.814 m，高 0.8 m，長 2.0 m，設(shè)計最大試驗風(fēng)速30 m/s。平板斷面風(fēng)洞測振采用模型如圖1所示，縱向長度0.8 m。為保證輕質(zhì)高強，材料采用碳纖維材質(zhì)。

圖1 平板模型(單位:mm)Fig.1 The model of plate(unit:mm)

模型的基本參數(shù):m=1.525 kg/m，Im=0.010 28 kg/m，豎向頻率fh=2.588 Hz，扭轉(zhuǎn)頻率fa=5.026 Hz。用激光位移計記錄模型試驗風(fēng)速下位移信號。

在平板測振試驗中，記錄的各級風(fēng)速下結(jié)構(gòu)位移響應(yīng)結(jié)果見圖2、圖3。風(fēng)速小于18 m/s時結(jié)構(gòu)位移振幅較小，基本處于靜止?fàn)顟B(tài)，風(fēng)速超過18 m/s后結(jié)構(gòu)振動位移值及方差突然增大，豎向振動參與程度不斷增加，直之風(fēng)速達20.4 m/s時結(jié)構(gòu)出現(xiàn)顫振失穩(wěn)現(xiàn)象。

圖2 豎向位移隨風(fēng)速變化關(guān)系Fig.2 Vertical vibration displacements vs wind velocity curve

圖3 扭轉(zhuǎn)振動位移響應(yīng)-風(fēng)速曲線Fig.3 Torsional vibration displacement vs wind velocity curve of plate section

圖4 位移幅值譜(20.4 m/s)Fig.4 Displacement amplitude spectrum(20.4 m/s)

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

2.1 平板顫振數(shù)值模擬方法

氣動彈性問題具有重要研究價值，當(dāng)橋梁結(jié)構(gòu)與空氣存在相對運動時，橋梁結(jié)構(gòu)會受到氣動力作用，氣動力會導(dǎo)致結(jié)構(gòu)位置及形狀發(fā)生改變。結(jié)構(gòu)位形的改變引起流體占據(jù)實際空間發(fā)生變化，反過來又導(dǎo)致作用在結(jié)構(gòu)上氣動力改變形成流體與固體的相互耦合作用。橋梁的顫振發(fā)散也是流固耦合作用的結(jié)果。

處理流固耦合問題可采用強耦合與松耦合兩種方法。強耦合求解方法以求解流體運動的迭代過程為主，將物體運動的求解耦合在流體運動求解中;松耦合方法為在每個時間步內(nèi)分別對流體與結(jié)構(gòu)依次求解，通過中間平臺交換數(shù)據(jù)信息，實現(xiàn)兩個場的耦合求解。此方法對計算機硬件要求不高，且可充分發(fā)揮CFD(Computational Fluid Dynamics)計算與CSD(Computational Structural Dynamics)計算的各自優(yōu)點，能實現(xiàn)數(shù)值計算模塊化［9］。分析橋梁結(jié)構(gòu)風(fēng)致振動問題時多采用松耦合的計算策略。

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

(1)進行流體域內(nèi)固定斷面的非定常繞流計算，提取靜止?fàn)顟B(tài)下初始流場及作用在斷面上的三分力，獲得t時刻氣動力。

(2)據(jù)計算所得升力、升力矩，通過數(shù)值方法對振動結(jié)構(gòu)進行結(jié)構(gòu)動力學(xué)求解。據(jù)newmarkk-β方法求解豎向及扭轉(zhuǎn)振動方程，獲得t+Δt時刻結(jié)構(gòu)運動的豎向速度與扭轉(zhuǎn)速度:

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

(4)重復(fù)(2)、(3)步，獲得橋梁斷面振動響應(yīng)時程。

(5)據(jù)斷面響應(yīng)時程曲線判斷振動是否已發(fā)散，若已發(fā)散則該風(fēng)速即為顫振臨界風(fēng)速，否則增加風(fēng)速按(1)～(4)步重新計算直到振動發(fā)散。

具體計算流程如圖5所示。

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

數(shù)值模擬用商業(yè)軟件Fluent提供的RANS方法的k-ωSST兩方程模型，計算域大小參考同濟大學(xué)土木工程防災(zāi)國家重點試驗室TJ-4風(fēng)洞中段試驗端設(shè)置，計算域沿流線長度為3 m(上游1 m，下游2 m)，橫向?qū)挾?.8 m。計算時壁面附近最小網(wǎng)格尺度為0.000 4 m，計算域用分塊結(jié)構(gòu)化網(wǎng)格。網(wǎng)格數(shù)量16.2萬。

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

采用數(shù)值方法對顫振過程進行模擬及預(yù)測，數(shù)值模擬獲得該平板顫振臨界風(fēng)速為21.2 m/s與試驗結(jié)果20.4 m/s尚有一定差距。圖6～圖8為平板運動時程及幅值譜，數(shù)值模擬所得顫振頻率為4.73 Hz，實測顫振頻率為 4.59 Hz。

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

圖6 豎向振動位移時程響應(yīng)(20.4 m/s，CFD數(shù)值模擬)Fig.6 Vertical vibration displacement vs time curve of plate section(20.5 m/s，CFD)

圖7 扭轉(zhuǎn)振動位移時程響應(yīng)(20.4 m/s，CFD 數(shù)值模擬)Fig.7 Torsional vibration displacement vs time curve of plate section(20.4 m/s，CFD)

圖8 位移幅值譜(20.4 m/s)Fig.8 Displacement amplitude spectrum(20.4 m/s)

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

3.1 分塊分析能量計算方法

分塊分析主要為詳細捕捉表面壓力變化，尋找引起顫振發(fā)散的主要能量輸入特點，考慮尾部旋渦對氣流能量輸入方式影響及影響范圍。對斷面分區(qū)，分別計算不同區(qū)域的能量輸入特點。單位區(qū)域能量為單位面積上一個周期內(nèi)氣動力輸入到系統(tǒng)的能量，其計算式為:

面積為S(圖9的網(wǎng)格線部分)的區(qū)域氣動力輸入能量w(s，t)計算式為:

式中:P(x，t)為結(jié)構(gòu)表面壓力，v(x，t)為結(jié)構(gòu)運動速度，n(x，t)為壓力與速度夾角。

圖9 分塊分析示意圖Fig.9 The block analysis

據(jù)本文提出的分塊分析思路編制計算程序，通過提取數(shù)值計算獲得的模型表面壓力，詳細研究一個周期內(nèi)氣流向結(jié)構(gòu)輸送能量過程，對模型表面進行分區(qū)見圖10。

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

圖11～圖14為平板扭轉(zhuǎn)振動時不同分區(qū)輸入到系統(tǒng)的能量。第1區(qū)域在前半個周期內(nèi)上表面消耗能量，下表面吸入能量，后半個周期則反之。一個完整周期內(nèi)輸入到系統(tǒng)的能量為正，但數(shù)值較小。從上下表面輸入到系統(tǒng)能量隨時間變化關(guān)系曲線上可看出，輸入的能量周期性明顯，如圖11所示。第4區(qū)域輸入系統(tǒng)的能量特點與1區(qū)基本相似，能量輸入也具有周期性特點，二者不同之處在于第4區(qū)域消耗系統(tǒng)能量。由于第4區(qū)及第1區(qū)靠近尾渦，因此受到尾部旋渦的周期性運動影響，該兩區(qū)域的能量輸入具有明顯的周期性特征(圖14)。

第2區(qū)域為系統(tǒng)提供了較多能量，在整個振動過程中下表面始終輸入正能量，上表面在前1/4周期消耗能量，后3/4周期內(nèi)吸入能量，總能量為正，數(shù)值較大，且增加較快，見圖12。第3區(qū)域即風(fēng)嘴部位為扭轉(zhuǎn)系統(tǒng)最主要能量來源，在一個完整周期內(nèi)上下表面均吸入系統(tǒng)能量(圖13)，且該區(qū)域吸入的能量遠大于其它區(qū)域?qū)ο到y(tǒng)能量的貢獻，為引起結(jié)構(gòu)振動發(fā)散的主要能量源。從第2、3區(qū)域能量隨時間變化關(guān)系曲線上可看出，該兩部位能量輸入周期性不明顯，受尾部旋渦運動影響較小。

豎向振動主要能量輸入則主要由第2區(qū)域引起，見圖15，各分區(qū)對扭轉(zhuǎn)振動能量的貢獻如圖16所示。通過上述分析知，風(fēng)嘴為扭轉(zhuǎn)振動的主要能量輸入部位，對系統(tǒng)扭轉(zhuǎn)振動穩(wěn)定性起重要作用，見圖17。

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

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

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

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

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

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

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

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

通過對CFD數(shù)值模擬計算結(jié)果研究分析發(fā)現(xiàn)平板流線型好的斷面周圍旋渦較難識別，由于振動模型周圍流場千變?nèi)f化，即使模型振動到同一位置，其周圍流場也完全不同。因此先用相位平均方法［10］得出模型振動4分點相位的流場圖，后用Snapshot POD旋渦識別方法［11］重建流場，所得流場特征見圖18。

由數(shù)值模擬結(jié)果看出，平板因流線型較好，振動模型周圍流場中旋渦識別較困難。據(jù)POD流場分解結(jié)果，平板處于顫振臨界狀態(tài)時，其尾部旋渦形狀為較規(guī)則圓形，旋渦渦心位于同一高度，呈直線排列，渦街結(jié)構(gòu)較穩(wěn)定，隨著模型振動而發(fā)生擺動，但渦心連線基本保持在同一條直線上。

由4區(qū)域輸入至系統(tǒng)扭轉(zhuǎn)振動能量隨時間的變化關(guān)系可推測出尾部旋渦對氣流輸入系統(tǒng)能量的影響。靠近尾流區(qū)域受尾部旋渦運動規(guī)律作用，輸入系統(tǒng)的能量周期性特點明顯。為更詳細研究尾流運動對結(jié)構(gòu)振動的影響，圖19給出4區(qū)域上下表面壓力分布隨時間的變化關(guān)系。為清楚說明尾部旋渦對模型不同部位影響，定義上下表面壓力差的相對變化量為:=(CPu-CPd)/CPu，其中CPu，CPd為模型上下表面壓力最大值。4區(qū)域上下表面的相對壓力差為:(4)=2.21，(1)=1.61，(2)=1.11，(3)=0(圖20)，由此看出，越靠近尾流區(qū)域，上下表面壓力差值的相對值越大，第4分區(qū)最大，而距尾流區(qū)最遠的第3分區(qū)上下表面之間的壓力相對差較小，兩曲線基本重合。即越靠近尾流區(qū)域，其表面壓力變化特點與尾流旋渦運動規(guī)律越接近。平板振動時尾部旋渦成直線排列并隨平板一起做大幅振動，會造成接近尾流區(qū)域的表面壓力相對變化較大，而遠離尾流區(qū)域影響較小。

圖18 平板斷面旋渦驅(qū)動結(jié)構(gòu)運動流(顫振風(fēng)速20.4 m/s，數(shù)值模擬)Fig.18 Vortex driven structural movement of plate section(20.4m/s，CFD)

圖19 平板表面4個分區(qū)上下表面壓力隨時間變化關(guān)系Fig.19 Four regional of plate girder surface pressure vs.time

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

4 結(jié)論

本文通過風(fēng)洞試驗測試了流線型較好平板斷面的顫振性能，采用CFD數(shù)值計算方法，并基于流固松耦合計算策略模擬了平板顫振過程，結(jié)合分塊分析思路定量研究氣流對振動模型表面不同區(qū)域的能量輸入特點及其對模型不同部位表面壓力影響，結(jié)論如下:

(1)分塊分析計算結(jié)果表明迎風(fēng)端風(fēng)嘴是扭轉(zhuǎn)振動能量的主要輸入部位，對系統(tǒng)扭轉(zhuǎn)振動穩(wěn)定性起重要作用;氣流輸入到平板的主要能量在一個周期內(nèi)不斷增加，導(dǎo)致平板顫振多為突然性失穩(wěn)。

(2)由分塊分析結(jié)果知，模型尾部風(fēng)嘴處旋渦對氣流能量輸入方式影響較大。平板處于顫振狀態(tài)時，尾部受旋渦的直接作用。尾端風(fēng)嘴處穩(wěn)定的渦街可控制模型尾端自由振動，使模型尾端運動與渦街?jǐn)[動趨勢基本同步;尾部渦街?jǐn)[動的規(guī)律性，造成模型靠近尾流區(qū)域能量輸入方式周期性明顯;而迎風(fēng)端的規(guī)律性較差，周期性不太明顯。

(3)分塊分析的結(jié)果表明平板尾部規(guī)律性擺動的渦街控制結(jié)構(gòu)運動，造成尾部風(fēng)嘴區(qū)域上下表面壓力相對波動較大，而遠離尾流區(qū)域影響較小，其上下表面壓力波動程度相對較小。

(4)平板顫振過程中在模型尾部產(chǎn)生的旋渦尺度較小，形態(tài)基本為圓形，氣流的能量分布在該尺度相當(dāng)?shù)男郎u之間，當(dāng)模型振動到平衡位置時，渦街位于同一高度，然后隨著模型的振動渦街發(fā)生相應(yīng)擺動，而渦心連線基本保持直線。

［1］Scanlan R H.The action of flexible bridges under wind，Ⅰ:flutter theory［J］.J.of Sound and Vibration，1978，60(2):187-199.

［2］ Larsen A.Aerodynamics of the tacoma narrows bridgc-60 years later［J］.Joumd of Structural Engineering International，2000，10(4):243-248.

［3］劉 高.大跨度懸索橋顫振分析的能量法［J］.中國公路學(xué)報，2000，12(3):20-24.

LIU Gao.Flutter analysis of long-span suspension bridges by energy method［J］.China Journal of Highway and Transport，2000，12(3):20-24.

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

［5］Nissim E.Flutter suppression using active controls based on the concept of aerodynamic energy［M］.NASA TND 6199，1997.

［6］Jones J G.On the energy characteristics of the aerodynamic matrix and the relationship to possible flutter［R］.The Aeronautical Quarterly，Part 3，1983:212-225.

［7］ Carta F O.Coupled blade-disk-shroud flutter instabilities in turbojet engine rotors［J］.Journal of Engineering for Power，1967，89(3):419-426.

［8］Klose A，Heinig K.A comparison of flutter calculations based on Eigenvalue and energy method［R］.AD-A2119741，1989.

［9］ 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):1471-1485.

［10］ 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(2):159-170.

［11］ Berkooz G，Holmes P，Lumley J L.Proper orthogonal decomposition in the analysis of turbulent flows［J］.Annual Review of Fluid Mechanics，1993，25:537-550.