劉祖軍，賈明曉，楊詠昕

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

箱梁是目前大跨度橋梁通常采用的截面形式。為提高其顫振臨界風速，往往可以通過改變截面的氣動外形來實現(xiàn)。但是箱梁的顫振與理想平板類似，多為彎扭耦合的顫振失穩(wěn)形態(tài)，其發(fā)生的細觀機理尚不夠清楚。

SCANLAN 等[1]引入航空領域顫振導數(shù)的概念，將作用在振動斷面上自激力表述為結構運動狀態(tài)與顫振導數(shù)的線性組合，并根據(jù)結構動力學方法建立了顫振運動控制方程，形成了完整的橋梁顫振分析理論。后來的研究者多以此為基礎，開展了橋梁的顫振性能研究。CHEN[2]提出的宏觀機理分析方法可以采用非迭代的手段研究顫振導數(shù)對系統(tǒng)顫振性能的影響。MATSUMOTO 等[3?4]提出了分步分析方法求解耦合的顫振控制方程，并以此為基礎，研究了多種斷面的顫振導數(shù)和氣動阻尼的變化規(guī)律。楊詠昕等[5]在Matsumoto 研究思路的基礎上，導出了二維三自由度耦合顫振分析方法，并用該方法研究了箱梁的顫振機理。丁泉順等[6]也采用類似方法，探討了扁平箱梁顫振的發(fā)生機理。這些對箱梁顫振機理的研究側重于分析顫振過程中顫振導數(shù)的變化規(guī)律。

除此之外，鮮榮等[7]采用風洞節(jié)段模型試驗研究了風嘴、欄桿、檢修軌道以及導流板對箱梁顫振穩(wěn)定性的影響。金挺等[8]通過對箱梁表面測壓數(shù)據(jù)進行積分的方法研究了箱梁斷面的雷諾數(shù)效應。任若松等[9]也采用表面測壓方法對準流線型橋梁斷面的表面壓力分布規(guī)律進行了研究。趙林[10]通過節(jié)段模型測壓試驗分析了在箱梁顫振過程中，模型表面壓力分布特性的變化規(guī)律。孟曉亮等[11]通過節(jié)段模型測振試驗和計算流體力學的方法，分析了不同風嘴角度對箱梁顫振和渦振性能的影響。張偉等[12]利用流場顯示技術研究了箱梁斷面顫振過程中，模型周圍流場的變化情況。這些研究分析了箱梁顫振過程中壓力場和流場的變化規(guī)律，但是沒有闡明模型表面壓力分布特性是如何影響箱梁的顫振導數(shù)以及顫振穩(wěn)定性的。ARGENTINI等[13]提出的分布式非定常氣動力表達方法則聚焦于表面風壓對顫振和抖振的影響。 LI 等[14]采用微觀機理分析方法，探討了表面壓力與顫振性能的聯(lián)系，并繪制了氣動阻尼的分布圖。

從研究現(xiàn)狀來看，目前從宏觀和細觀兩個層面上都對箱梁的顫振機理進行了研究，但兩個層面之間缺乏聯(lián)系，很少將體現(xiàn)自激作用的顫振導數(shù)與模型表面的分布壓力關聯(lián)起來，因此也就缺乏對箱梁顫振機理深入認識。本文推導了顫振導數(shù)與模型表面壓力之間的數(shù)學關系，研究了分布壓力的空間作用位置以及相位分布特性對箱梁顫振導數(shù)和顫振臨界風速的影響。通過對模型表面進行分區(qū)，研究了壓力的局部特性對顫振導數(shù)和系統(tǒng)振動能量的影響。本項研究加深了對箱梁顫振發(fā)生機理的認識，為有效實施箱梁的顫振控制提供理論基礎。

在空氣中振動的橋梁斷面，其運動狀態(tài)與作用在模型上的自激氣動力密切相關，而自激力是模型表面分布壓力的合力。因此，為了深入理解箱梁斷面顫振發(fā)生的本質就需要探討箱梁顫振過程中模型表面分布壓力的變化規(guī)律，并建立顫振導數(shù)與分布壓力之間的數(shù)學關系式。

將推導模型表面分布壓力與顫振導數(shù)之間的聯(lián)系，通常采用8 個顫振導數(shù)來描述的氣動升力和升力矩[15?17]，如式(1)所示：

1 顫振導數(shù)與表面壓力的聯(lián)系

1.1 顫振導數(shù)與表面壓力的數(shù)學關系

式中：φα為氣動力與模型扭轉運動之間的相位差；x為沿模型寬度方向上的任意一點與模型中心之間的距離。

聯(lián)立式(3)～式(4)、式(10)～式(11)，并通過比較系數(shù)可以求出H2?、H3?、A?2、A?3這4 個顫振導數(shù)：

需要特別說的是本文在進行上述公式推導時假設橋梁顫振臨界狀態(tài)的振動形式為等幅的簡諧振動。其主要是基于以下幾點考慮：

1)目前對于橋梁―氣流系統(tǒng)線性顫振臨界狀態(tài)的分析通常是認為當處于顫振臨界點時，系統(tǒng)的阻尼為零，系統(tǒng)的振動為等幅簡諧振動，氣流輸入到系統(tǒng)的能量剛好等于系統(tǒng)的阻尼耗能。例如：塔科馬大橋在墜毀之前就經(jīng)歷的較長時間的扭轉振動，較多學者在對該橋的顫振機理進行分析時也大多假設結構的振動為等幅的簡諧振動。另外，在采用強迫振動試驗確定顫振導數(shù)時，也通常是強迫斷面做等幅的簡諧振動，從而獲得作用在斷面上的氣動力，進而識別顫振導數(shù)。因此本文基于現(xiàn)有的方法，仍假設橋梁―氣流系統(tǒng)線性顫振臨界狀態(tài)的運動形式為等幅的簡諧運動。

2)這一假設僅適用于橋梁發(fā)生小幅振動，并且不考慮系統(tǒng)非線性影響因素的條件下。對于線性系統(tǒng)而言，如果系統(tǒng)的振動位移為等幅的簡諧振動，那么輸入到系統(tǒng)的氣動力也應該是簡諧振動形 式，氣動力與振動位移之間只存在相位差。由于氣動力是模型表面壓力積分之和，因此如果氣動力是按照簡諧振動方程的形式發(fā)生變化，則模型表面的分布壓力也可以近似認為是按照同樣的規(guī)律發(fā)生變化。但是如果強烈的流固耦合作用使得氣動力具有較強的非線性時(氣動力含有高階成份)，這些假設條件就不再適用了。

1.2 數(shù)值實現(xiàn)途徑

通過上述的理論分析可知，采用模型表面壓力來求解顫振導數(shù)的關鍵是如何獲取模型表面壓力系數(shù)的幅值Cp?和相位差φ，結合數(shù)值方法研究具體的實現(xiàn)過程，步驟如下(圖1)：

圖1 數(shù)值實現(xiàn)途徑Fig. 1 The process of numerical simulation

1)設定振動頻率ω，振動幅值及來流風速，采用數(shù)值方法結合動網(wǎng)格技術模擬結構的強迫振動。

2)數(shù)值模擬時當結構運動狀態(tài)達到平穩(wěn)后，提取模型表面任意一點N個時刻的壓力P，設x為n個 時刻，y為 與之相對應的壓力P， 即：x=[t1,t2,···,tn]，y=[P1,P2,···,Pn]。

3)針對模型表面的任意一點，采用函數(shù)y=asin(ωx)+bcos(ωx)+c來擬合步驟2 中提取的N個時刻的壓力值P，獲得擬合函數(shù)的系數(shù)a、b、c。

上述建立了顫振導數(shù)與模型表面壓力之間的聯(lián)系，為分析模型表面壓力分布特性對箱梁顫振穩(wěn)定性的影響提供了基礎。需要說明的是采用表面壓力來描述顫振導數(shù)時，需要獲得壓力幅值Cp?以及壓力與振動位移之間的相位差φ。這里的相位差定義為箱梁表面壓力的峰值與模型振動最大位移之間的相位延遲。

2 箱梁斷面風致振動的數(shù)值模擬

根據(jù)上述建立的模型表面壓力與顫振導數(shù)之間的數(shù)學關系，采用計算流體動力學方法(CFD)模擬了箱梁的強迫振動，提取了模型表面不同時刻的壓力值，結合本文1.2 節(jié)中提出的分析方法，研究了分布壓力對顫振導數(shù)的影響。

箱梁斷面以大海帶橋的主梁為原型進行了縮尺，其幾何尺寸如圖2 所示。應用商業(yè)軟件Fluent進行數(shù)值模擬計算，計算域的選取如圖3 所示，計算時壁面附近最小網(wǎng)格尺寸為0.0004 m，計算域采用分塊結構化網(wǎng)格。網(wǎng)格數(shù)量為15.6 萬。湍流模型采用基于RANS 方法的k?ωS S T兩方程模型，壓力速度耦合采用SIMPLE 算法，求解器采用分離式?；诹鞴倘躐詈系挠嬎悴呗院蛣泳W(wǎng)格技術分別模擬了箱梁的豎向和扭轉強迫振動。

圖2 箱梁斷面 /mmFig. 2 The box girder section

圖3 二維計算區(qū)域及邊界設置示意Fig. 3 Two-dimensional computational domain and boundary

結合該模型的顫振風洞試驗結果，設置計算風速為20.8 m/s，豎向振幅為1.8 cm，扭轉幅值為2.4°，振動頻率為6.05 Hz。計算域的上端和下端設為壁面邊界條件。計算時間步長取 ?t=0.004 s，每次迭代計算的殘差小于 10?5認為計算結果收斂。調(diào)整近壁面的第一層網(wǎng)格尺寸，使y+<5(圖4)。

圖4 模型斷面近壁面網(wǎng)格的y+值Fig. 4 The grid value of y+ on the model surface

3 箱梁表面壓力分布特性與顫振導數(shù)之間的聯(lián)系

3.1 模型表面壓力分布特點與顫振導數(shù)的關系

根據(jù)1.2 節(jié)中的方法，研究了箱梁表面壓力分布特點與顫振導數(shù)之間的關系。模型運動的相位定義如圖5 所示。為了說明數(shù)值擬合的精度，圖6給出了模型分別處于平衡位置(T/2 時刻)，此時模型具有最大速度，以及模型運動到最大位移時(3T/4 時刻)，模型上下表面壓力擬合的殘差沿模型寬度方向的分布。從圖6 可以看出，數(shù)值擬合的殘差絕大部分在5%以內(nèi)，只是在個別點處擬合的殘差超過了5%，因此擬合的精度滿足要求。另外，圖7 給出了T/2時刻，模型上下表面壓力的計算值與擬合值，從圖中可以看出，二者之間吻合較好。

圖5 模型的振動狀態(tài)Fig. 5 The model vibration state

圖6 壓力數(shù)值擬合的殘差在模型表面的分布Fig. 6 The pressure numerical fitting residuals' distribution on the model surface

圖7 模型表面壓力的計算值與擬合值( T/2)Fig. 7 The model surface pressure calculation result and fitting values ( T/2)

上述分析說明了本文方法的有效性。根據(jù)CFD 計算獲得的模型運動狀態(tài)下的氣動力，結合本文1.2 節(jié)中的方法，計算了擬合系數(shù)c沿箱梁表面的分布情況，如圖8 所示。擬合系數(shù)c的分布情況代表了箱梁表面平均壓力的分布特征。從圖8可以看出箱梁迎風側風嘴附近的平均壓力較大，上下表面壓力的分布規(guī)律類似，但力的作用方向相反。

圖8 擬合系數(shù)c 沿模型表面分布圖Fig. 8 The distribution of fitting coefficient c along the model surface

圖10 給出了擬合系數(shù)b沿模型表面的分布情況，同圖9 相比，二者的分布規(guī)律相似，也是在模型迎風側風嘴所在區(qū)域數(shù)值較大，而在箱梁尾部所在區(qū)域數(shù)值較小。由于系數(shù)b的分布與氣動導數(shù)A?2直接相關聯(lián)，因此箱梁迎風端的外形對顫振導數(shù)A?2影響很大，而顫振導數(shù)A?2對系統(tǒng)振動的穩(wěn)定起到了非常重要的作用。

圖9 擬合系數(shù)a 沿模型表面分布圖Fig. 9 The distribution of fitting coefficient a along the model surface

圖10 擬合系數(shù)b 沿模型表面分布圖Fig. 10 The distribution of fitting coefficient b along the model surface

3.2 箱梁表面壓力的局部特性對顫振導數(shù)的影響

通過3.1 節(jié)的分析可知在箱梁振動過程中模型表面壓力的分布是不均勻的，模型表面的不同區(qū)域的分布壓力對顫振導數(shù)的影響也是不相同，因此模型表面壓力的局部特性對箱梁顫振穩(wěn)定性產(chǎn)生的影響也不盡相同。因此，本節(jié)提出通過分塊分析的方法來定量研究模型表面壓力的局部特性對顫振導數(shù)的影響，模型表面分區(qū)如圖11 所示。

圖11 箱梁斷面的分區(qū)示意圖Fig. 11 The partition of H-shape section

根據(jù)CFD 數(shù)值計算結果，提取了箱梁斷面在一個振動的周期內(nèi)模型表面4 個分區(qū)上下表面壓力差隨計算時間變化關系，如圖12 所示。從圖12可以看出，在一個振動周期內(nèi)，位于模型的迎風側的3 區(qū)和2 區(qū)的壓力差的變化隨時間呈現(xiàn)出較明顯的波動性，這與模型的周期性運動相對應。而位于模型背風側的1 區(qū)和4 區(qū)的壓力差基本上不隨時間發(fā)生較大的變化。

圖12 箱梁斷面不同區(qū)域的壓力差Fig. 12 The pressure lag of different area on the box girder surface

結合建立的顫振導數(shù)與模型表面壓力之間關系，討論模型表面不同區(qū)域的分布壓力對顫振導數(shù)的貢獻。

圖13 模型表面各分區(qū)對顫振導數(shù)的貢獻Fig. 13 Each partition contribution of the model to aerodynamic derivatives

通過計算模型表面不同區(qū)域輸送給主梁的氣動能量來更加清楚的說明模型表面分布壓力的局部特性對主梁振動產(chǎn)生的影響。

一個周期內(nèi)單位面積上的氣動力輸送給振動系統(tǒng)的非定常平均積累功w(x,t) 和功率N(t)的定義為：

式中：P(x,t) 為 結構表面的壓力；v(x,t)為結構運動速度；n(x,t)為壓力與速度之間的夾角。

根據(jù)式(21)計算了在一個振動周期內(nèi)箱梁表面的不同區(qū)域輸送給主梁的氣動能量(圖14)。迎風側風嘴所在的3 區(qū)對模型輸入的能量最大，2 區(qū)次之，相比之下1 區(qū)輸入的氣動能量則較小，而位于模型背風側風嘴處的4 區(qū)輸入給主梁的氣動能量會由正變?yōu)樨?，因此有利于系統(tǒng)振動的穩(wěn)定性。從這一分析結果可以看出迎風側風嘴附近的分布壓力對運動斷面輸入正能量，而背風側風嘴附近的分布壓力消耗運動斷面的能量，因此在模型振動過程中，如果模型表面的主要分布壓力向模型迎風側移動，則氣動力輸入到系統(tǒng)的正能量逐漸增加，最終會導致模型振動穩(wěn)定性的喪失。

圖14 箱梁斷面不同區(qū)域輸送給振動系統(tǒng)的能量Fig. 14 The energy input to the system by the box girder different areas

4 探討箱梁表面壓力分布特性對顫振穩(wěn)定性的影響

以上分析了箱梁表面壓力分布特點與顫振導數(shù)的之間聯(lián)系，以及分布壓力的局部特性對顫振導數(shù)的影響，以此為基礎來研究壓力的分布特性是如何對箱梁顫振穩(wěn)定性產(chǎn)生影響。為了便于說明問題，建立如圖15 所示的箱梁斷面的坐標系，靠近來流風場所在的箱梁區(qū)域為迎風側，而靠近箱梁尾部的區(qū)域為背風側。

圖15 箱梁斷面的坐標系Fig. 15 The coordinate system of box girder

4.1 提高箱梁顫振穩(wěn)定性所需要的壓力分布特性

圖16 箱梁表面壓力和相位差的分布Fig. 16 The distribution of surface pressure and phase lag long the box girder

4.2 壓力分布特性對箱梁顫振臨界風速的影響

4.1 節(jié)初步探討了為提高箱梁顫振臨界風速所需的壓力分布規(guī)律。結合箱梁模型的具體參數(shù)，根據(jù)假設的模型表面壓力分布，應用本文的式(12)、式(13)計算顫振導數(shù)，通過二維耦合顫振分析方法，計算當分布壓力處在模型表面不同位置，并且處于不同的相位差分布時箱梁的顫振臨界風速。箱梁截面仍然采用如圖2，縱向長度0.8 m。模型的基本參數(shù)如下：每延米的質量m=1.3625 kg/m，每延米的質量慣性矩Im=0.01277 kg·m2/m，豎向振動頻率fh=2.825 Hz ， 扭轉振動頻率fa=7.452 Hz，豎彎阻比 ξh=0.0127 ， 扭轉阻比 ξa=0.047。

首先分析分布壓力的空間位置對模型顫振臨界風速的影響。假設模型表面壓力的分布范圍是0.2b，壓力的幅值Cp?=10.0，壓力中心距截面形心的距離xu分別為 ? 0.9b、 ? 0.4b、 0 .4b、 0 .9b(圖17(a))。相位差呈對稱分布，如圖17(b)所示，在靠近迎風側φa=π/3、φh=π/4 ； 而在背風側φa=?π/3、φh=?π/4。根據(jù)上述壓力分布規(guī)律，應用本文推導的顫振導數(shù)與分布壓力之間的關系，可以計算出上述4 種壓力分布情況下的顫振導數(shù)，應用二維耦合顫振分析，求得了與上述4 種壓力分布狀況相對應的顫振臨界風速(圖17(c))。從圖中可以看出當分布壓力位于模型的迎風側時，箱梁的顫振臨界風速最小，數(shù)值為16 m/s， 而當分布壓力位于模型的背風側時，箱梁的顫振臨界風速最高，數(shù)值為36.5 m/s。因此，當箱梁表面峰值壓力逐漸向模型尾部移動時，箱梁的顫振臨界風速會逐漸提高。這也意味著如果箱梁尾部存在較強的壓力作用則有利于箱梁顫振的穩(wěn)定性。

圖17 箱梁表面的壓力分布特性和顫振臨界風速Fig. 17 The pressure distribution characteristics and flutter critical wind speed of the box girder

分析相位差的分布對箱梁顫振臨界風速的影響，這里進行了兩種情況下的計算：第1 種情況下，xu分別取為? 0.9b、? 0.4b、0 .4b、0 .9b，相位差呈反對稱分布，在靠近迎風側φa=π/3、φa=?π/4；在背風側φa=?π/3、φh=+π/4，如圖18所示。通過計算發(fā)現(xiàn)在這種相位差分布情況下，上述4 種壓力分布狀況對應的顫振臨界風速高達288.8 m/s，也就意味著很難發(fā)生顫振失穩(wěn)。

圖18 相位差沿箱梁表面成反對稱分布Fig. 18 The anti-symmetrically distributed phase lag along the box girder

第2 種情況下的計算條件為：?0.9b≤xu1≤?0.4b、0.4b≤xu2≤0.9b；相位差φa與φh分布規(guī)律相同，在迎風側 π/6≤φa=φh≤π/2，在背風側?π/2≤φa=φh≤?π/6。通過二維耦合顫振分析方法獲得箱梁顫振臨界風速如圖19 所示。

圖19(a)給出了當分布壓力位于模型迎風側時(?0.9b≤xu1≤?0.4b) 相位差由 π/6 變 化到 π/2的過程中，箱梁顫振臨界風速變化規(guī)律。從圖19(a)可以看出當迎風側的分布壓力向模型下游移動時，顫振臨界風速逐漸增加。而當分布壓力空間位置固定時，相位差由 π/6 變化到 π/2的過程中，箱梁顫振臨界風速先增加后減小。

圖19(b)給出了當分布壓力位于模型背風側時( 0.4b≤xu2≤0.9b) 相位差由 ?π/2 變 化到 ?π/6的過程中，箱梁顫振臨界風速變化規(guī)律。從圖19(b)可以看出當模型背風側的分布壓力向模型尾部移動時，箱梁的顫振臨界風速增加；而當分布壓力空間位置固定時，在相位差由 ?π/2 變化到?π/6的過程中，箱梁顫振臨界風速風速也不斷增加。

圖19 顫振臨界風速隨箱梁表面壓力分布特性的變化規(guī)律Fig. 19 The change of flutter critical wind speed with the distribution characteristics of the box girder surface pressures

通過探討表面分布壓力的空間位置和相位差對箱梁顫振臨界風速的影響?？梢缘贸鲆恍﹩⑹荆?/p>

1)顫振導數(shù)與模型表面壓力的空間位置和相位差密切相關。在峰值壓力向箱梁尾部移動過程中，顫振導數(shù)A?1和A?2的絕對值會不斷減小，但是壓力分布的空間位置對顫振導數(shù)H2?和H3?卻沒有影響。箱梁迎風側的分布壓力對模型的顫振穩(wěn)定性起了不利的作用，而當迎風側分布壓力向模型尾部移動時，箱梁的顫振臨界風速將會不斷的提高。這也說明了如果箱梁尾部存在較強分布壓力，則會有利于箱梁顫振的穩(wěn)定性。

2)相位差的分布對所有的顫振導數(shù)都產(chǎn)生影響。相位差φa、φh呈反對稱分布則能夠對振動系統(tǒng)的穩(wěn)定性產(chǎn)生有利作用，提高了系統(tǒng)的顫振臨界風速。如果φa與φh的空間分布規(guī)律相同，那么當分布壓力位于模型迎風側時，顫振臨界風速隨相位差的增加先增大后減?。划敺植級毫ξ挥谀Ｐ捅筹L側時，顫振臨界風速隨相位差的增加而增大。

5 結論

本文建立了箱梁表面壓力與顫振導數(shù)之間的數(shù)學關系式，應用CFD 數(shù)值模擬方法獲得了箱梁振動狀態(tài)下的表面壓力，研究了箱梁表壓力與顫振導數(shù)的聯(lián)系，探討了表面壓力分布的空間位置和相位差的分布規(guī)律對箱梁顫振穩(wěn)定性的影響。通過本項研究得出以下幾點結論：

(1)對主要顫振導數(shù)A?1、A?2、H3?貢獻較大的是箱梁迎風側風嘴所在的區(qū)域，并且決定A?1和H3?的表面壓力的分布規(guī)律類似。當豎向振動與扭轉振動的相位差不斷減小時，決定A?1和H3?的表面壓力的分布規(guī)律也就基本趨于相同，從而使得A?1和H3?的耦合作用更為強烈，因此，氣動耦合效應與模型表面壓力的分布特性密切相關。

(2)箱梁迎風側風嘴及其附近區(qū)域的壓力特性決定了箱梁主要顫振導數(shù)的性質，這對有效實施風振控制措施具有指導意義，通過合適的氣動控制措施改變這些區(qū)域的流場特性，則能影響壓力在這些區(qū)域的分布規(guī)律，從而提顫振臨界風速。

(3)箱梁表面壓力分布的空間位置對箱梁的顫振導數(shù)和顫振臨界風速都有重要的影響。表面分布壓力位于箱梁迎風側時，則對箱梁的顫振穩(wěn)定起了不利的作用，而當箱梁尾部存在較強的分布壓力時，箱梁的顫振臨界風速將會提高。當迎風側的分布壓力向箱梁尾部移動時，顫振導數(shù)A?1和A?2的絕對值會不斷減小，而顫振導數(shù)H2?和H3?不受此影響，箱梁的顫振臨界風速會逐漸增加。

(4)箱梁表面壓力的相位差沿模型寬度的分布特性會對箱梁的所用顫振導數(shù)都產(chǎn)生影響。當表面壓力與扭轉和豎向振動最大位移之間的相位差φa、φh呈反對稱分布時，系統(tǒng)的顫振臨界風速較高。而當φa與φh的分布特性相同時，箱梁顫振臨界風速隨相位差的變化規(guī)律與模型表面壓力所處的空間位置有較大關系。