拱橋鋼箱吊桿馳振穩定性數值計算研究*

詹 昊 廖海黎

(西南交通大學土木工程學院1) 成都 610031) (中鐵大橋勘測設計院集團有限公司2) 武漢 430056)

0 引 言

通常馳振是細長結構在橫風向作用下的一種不穩定的單自由度發散振動現象，會導致結構破壞.對于橋梁結構，馳振易發生在拱橋吊桿、橋塔這類阻尼較小的細長結構上.由于馳振具有很大的破壞性，應采取措施盡力避免馳振發生，使馳振臨界風速高于馳振檢驗風速.對于不同切角形狀的正方形和長方形截面的氣動力性能研究相對較少，且多為靜態繞流研究，可以參見文獻[1-4]的風洞實驗結果及文獻[5]的數值仿真計算結果.對于馳振問題的數值計算主要有兩種方法.方法一，首先通過計算流體方法求出不同風攻角下的阻力系數和升力系數，再根據葛勞渥-登哈托原理判斷馳振是否發生，用相關公式計算馳振臨界風速.方法二，建立流固耦合數值計算模型直接吹出馳振臨界風速.從相關文獻來看，多用方法一進行研究，運用方法二進行研究相對較少[6].

本文通過兩種方法對橫橋向來風時，南京大勝關大橋吊桿進行馳振數值仿真計算分析，并對其截面進行優化.兩種方法相互驗證，與風洞實驗結果進行對比分析.并展示了馳振中的極限環振動現象.

1 工程概況

南京大勝關橋為雙線高速和雙線I級干線鐵路橋梁，主橋采用6跨連續鋼桁拱橋，三桁承重結構.大橋吊桿采用鋼箱型截面，最長吊桿長度超過50 m，柔度大，發生馳振的可能性大.對于拱橋吊桿的風致振動制振措施一般是改變吊桿氣動外形的空氣動力學措施和設置阻尼器的機械措施.機械措施所需成本及后期維護費用很高；截面形狀的微小變化往往會改變結構的氣動性能，在保證吊桿安全的情況下，改變吊桿氣動外形的氣動力學措施幾乎無需增加成本，因此，吊桿截面選型顯得十分重要.南京大勝關大橋三種截面方案見圖1，計算參數為表1.

圖1 三種方案截面圖 (單位：mm)

參數實橋長方形切角長方形A切角長方形B 質量/(kg·m-1)1 09910.9910.9910.99 順橋向豎彎自振頻率/Hz2.5442.5442.5444.10阻尼比0.0020.00280.00280.0012 桿件迎風面尺寸H/m1.120.1120.1120.1032

2 方法一

2.1 數值仿真計算原理

風軸坐標系下的阻力、升力和升力矩系數分別為

(1)

(2)

(3)

式中：0.5ρU2為氣流動壓；L為節段模型的長度；FD(α)，FL(α)和MZ(α)分別為攻角α情況下采用風軸坐標系時的阻力、升力和俯仰力矩；H為結構高度.

馳振臨界風速為

(4)

式中：m，ζ和ω分別為結構的質量，阻尼比和角頻率；H為結構高度；ρ為流體密度，空氣密度常溫時取1.225 kg/m3.

2.2 數值仿真計算模型

長方切角B截面網格劃分圖見圖2，其他兩種截面網格劃分與長方形截面相似.整個流場網格由三部分構成，包圍截面的邊界層結構網格區域，非結構網格過渡區域和最外面的結構網格區域.動網格變形采用動態分層法.風向從左至右，左側設定為速度入口，右側設定為自由出流.上下邊界為無滑移固壁邊界.矩形計算區域為3.3 m×3 m，模型上游來流區域為0.8 m，下游尾流區為2.5 m，離上下邊界各為1.5 m.數值計算中，采用有限體積法求解，其中對流項采用二階迎風差分格式，壓力和速度的耦合采用SMPLEC算法.RNGk-ε湍流模型對于有較大速度梯度的流場及分離流計算較準確，計算時間適中，本文計算采用RNGk-ε湍流模型.

圖2 計算網格劃分

2.3 數值仿真計算結果

橫橋向來風時，阻力系數CD和升力系數CL隨風攻角α變化曲線圖見圖3～4.

圖3 阻力系數隨風攻角變化曲線(U=10 m/s)

圖4 升力系數隨風攻角變化曲線(U=10 m/s)

由圖3～4可知，CD大致以0°風攻角為中心對稱分布，CL大致以0°風攻角為中心反對稱分布.長方形截面CD計算值大于風洞試驗值，CL計算值和風洞試驗值吻合較好；切角長方形截面A的CL和CD計算值和風洞試驗值吻合較好；切角長方形截面BCD和風洞試驗值吻合較好，CL數值計算值和風洞試驗值有所差別，但變化趨勢基本一致.根據馳振臨界風速公式計算的結果見表2.

由表2可知，馳振力系數和風洞試驗基本吻合，通過式(4)計算得到的馳振臨界風速兩者基本吻合.橫橋向來風時，實際橋梁最長吊桿中間高度的馳振檢驗風速為50.4 m/s，長方形截面和長方形截面A不能滿足馳振要求.長方形截面B由于馳振力系數為正數，不會發生馳振.

表2 方法1計算結果

3 方法二

3.1 數值仿真計算原理

將吊桿簡化成質量彈簧阻尼系統，建立豎向彎曲流固耦合模型，數值仿真計算原理示意圖見圖5[7-8].

圖5 數值仿真計算原理示意

(5)

(6)

式中：m為結構的質量；kh為豎彎慣性矩；kh為豎向剛度；ch為結構阻尼；Fh為豎向氣動力；y為物體豎向位移.

對于不可壓縮流體的連續方程和納維-斯脫克斯方程為

(7)

用FLUENT軟件求解流體控制方程式(6)～(7)，得到吊桿周圍流體的速度場和壓力場以及作用在柱體上的升力.將Newmark方法的代碼嵌入用戶自定義函數(UDF)同軟件連接，通過UDF來提取升力代入吊桿振動方程(2)，求解吊桿的動力響應.然后將吊桿的速度通過FLUENT的動網格技術進行傳遞使網格獲得速度.最后用速度與時間步相乘來得到網格位置的更新，開始下一個時間步計算.如此循環得到各時間步的振動位移.計算網格劃分同方法1.

3.2 數值仿真計算結果

各截面方案旋渦脫落見圖6，不同風速下長方形截面A位移時程曲線見圖7.

圖6 旋渦脫落圖(U=10 m/s)

由圖6可知，長方形截面切角后，旋渦脫落強度減小，其中切角長方形截面B的渦脫落強度最小.

圖7 不同風速下吊桿位移時程曲線(切角長方形截面A)

由圖7可知，U=7 m/s，豎向振幅很小，不到1 mm；U=8.0m/s，豎向振幅突然增大，達到了65 mm；隨著風速的增大，豎向振幅繼續增大，最后穩定在一個極限振幅上.與渦激共振不同，馳振沒有鎖定風速區間，其振幅隨著風速的增大而增大，且比渦激共振振幅大得多，達到桿件截面尺寸數倍.不同風速下長方形截面和切角長方形截面A最大振幅見表3～4.

表3 不同風速下最大振幅(長方形截面)

H=11.2 cm，為桿件迎風面高度.

表4 不同風速下最大振幅(切角長方形截面A)

以結構振幅開始迅速增大時的風速作為馳振臨界風速，由表3～4可知，長方形截面馳振臨界風速為7.5～8.5 m/s, 切角長方形截面A馳振臨界風速為7～8 m/s，切角長方形截面B在各級風速下振幅沒有超過2 mm，沒有發生馳振.方法二計算結果與風洞試驗值對比見表5，節段模型風洞試驗采用直接吹出馳振臨界風速.方法一計算結果見表6，兩種方法計算結果和風洞實驗基本吻合.

表5 方法二計算結果

表6 方法一計算結果

風速U=16m/s時，切角長方形截面A不同時間段的升力系數時程曲線及對應的功率譜密度函數分別見圖8～9，升力系數功率譜密度曲線峰值對應的橫坐標為旋渦脫落頻率值.

0～35 s是桿件靜止到振動到達穩態的整個過程：其中0～5 s是桿件從靜止到開始振動；17.5～35 s是物體振動達到穩定狀態，其中30～32 s是穩態中的一段放大圖像.由圖8可知，0～5 s時，渦脫落頻率f=27 Hz，這與截面靜止繞流時的渦脫落頻率相同.17.5～35 s時物體振動達到穩定狀態，由于物體大幅度振動，對渦脫落頻率有調整作用，渦脫落頻率和物體的固有頻率趨向一致，這與文獻[9-10]中的結論一致.升力系數曲線大的波動周期與物體的固有頻律相近，桿件的固有頻率是2.544 Hz，此時主要旋渦脫落頻率f=2.5 Hz.同時還存在3倍、5倍、10倍和15倍的次要頻率，即存在f=7.5，12.5，22.5，27.5 Hz的渦脫落頻率.方法二的計算結果見表5. 當V=24 m/s，切角長方形截面A相軌跡見圖10.

圖8 不同時間段的升力系數時程曲線(切角長方形截面A,V=16 m/s)

圖9 不同時間段升力系數功率譜密度曲線(切角長方形截面A,V=16m/s)

圖10 相平面圖(切角長方形截面A,V=24 m/s)

由圖10可知，相軌跡逐漸擴大，最后穩定在一個極限環內.

4 結 論

1) 運用FLUENT計算吊桿的阻力系數和升力系數隨風攻角變化值和風洞實驗值基本吻合，通過式(4)計算得到的馳振臨界風速兩者基本吻合； 通過對FLUENT二次開發，運用流固耦合方法計算得到的馳振臨界風速和節段模型風洞實驗吹出的馳振臨界風速基本吻合.可見這兩種數值計算方法具有較高的精度，為今后馳振研究提供了有力工具.

2) 長方形截面切角后，旋渦脫落強度減小，其中切角較大的長方形截面B的渦脫落強度最小.流動形態的改變導致結構受力性能的改變，影響結構的馳振穩定性.

3) 長方形截面切角能否改善馳振性能取決于切角的形狀和大小，對切角的形狀和大小非常敏感.當順橋向來風時，長方形截面和切角長方形截面A在較高風速時會發生弛振，表現為隨風速增大振幅逐漸增大的極限環振蕩.切角長方形截面B在不會發生弛振.為今后工程設計提供了參考.