桁架主梁橋梁顫振流固耦合數值仿真計算*

詹 昊 郭子會 閆斗平 廖海黎

(1.西南交通大學土木工程學院 成都 610031； 2.內蒙古伊泰準東鐵路有限責任公司 鄂爾多斯 010300；3.鄖縣錦宏路橋工程有限責任公司 十堰 442500； 4.中鐵大橋勘測設計院集團有限公司 武漢 430056)

顫振是一種空氣動力失穩現象，會引起災難性的后果。橋梁結構一般為非流線型斷面，當氣流流經時會發生渦旋和流動分離，形成復雜的空氣作用力，當空氣力受結構振動的影響較大時，同振動結構形成一個相互作用反饋機制，動力系統的空氣主要表現為一種自激力。當橋梁結構與產生自激氣動力的繞流氣流形成振動系統阻尼在不斷的相互反饋作用中由正值變為負值時，振動系統吸收的能量超越了自身耗能能力而造成系統振動發散，這種空氣動力失穩現象就是橋梁顫振。塔科瑪海峽橋風毀就是由顫振所引起。

隨著橋梁跨度的增加，其結構對風的敏感性也逐漸增加，因此，顫振穩定性往往是橋梁設計中首要考慮的因素。桁架主梁橋梁廣泛運用于公鐵兩用大橋、雙層公路橋梁和山區橋梁。對于基于流固耦合理論的顫振計算，國內外學者作了比較深入的研究，但由于桁架主梁構造復雜，一般簡化為平面模型計算。文獻[1-2]將鋼箱主梁簡化為平面模型，運用離散渦方法進行顫振計算，文獻[3]將鋼箱主梁簡化為平面模型，運用有限元方法進行顫振計算，文獻[4]將鋼桁主梁簡化為平面模型，運用無網格渦方法進行主梁靜力三分力系數計算，文獻[5]將桁架主梁簡化為平面模型，運用有限體積法計算主梁顫振導數，文獻[6-7]通過風洞實驗研究桁架主梁橋梁顫振穩定性。

黃岡長江大橋為主跨566 m的雙塔雙索面斜拉橋，主跨跨度、主桁桿件傾斜度、斜拉索破斷力和抗壓、抗拉支座均居世界已建成橋梁之首。黃岡公鐵兩用長江大橋效果圖見圖1，橫斷面圖見圖2。

圖1 黃岡公鐵兩用長江大橋效果圖

圖2 黃岡公鐵兩用長江大橋橫斷面布置圖(單位：m)

本文通過對FLUENT二次開發，建立了基于幾何三維的豎彎和扭轉流固耦合數值仿真計算模型，對黃岡長江大橋進行了顫振數值仿真計算，并與節段模型風洞試驗結果進行了比較。建立桁架主梁三維模型，運用流固耦合方法計算顫振臨界風速，國內外還未見相關文獻報道。該橋顫振檢驗風速為60.4 m/s，設計要求顫振臨界風速大于顫振檢驗風速。

1 數值仿真計算原理

將主梁作為質量、彈簧和阻尼系統，如圖3所示。研究對象為剛體，假設物體沿展向的豎向位移和扭轉位移相同。豎向振動方程和扭轉振動方程如式(1)和式(2)所示。

(1)

(2)

式中：m為單位長度的質量；Iθ為慣性矩；kh為豎向剛度；kθ為扭轉剛度；ch和cθ分別為結構豎彎阻尼和扭轉阻尼；Fh和Mθ為流體力；y為結構豎向位移；θ為扭轉角度。對于不可壓縮流體的連續方程和納維-斯脫克斯方程如式(3)和式(4)所示，梯度算子如式(5)所示。

U=0

(3)

Ut+(U·p+υ2U

(4)

(5)

流體質點速度向量U={u，υ，w}T；X，Y，Z為笛卡爾坐標系；ρ為流體密度，i，j，k為單位向量。假設空氣的密度ρ=1.225 kg/m3,氣溫20 ℃時，運動黏度υ為1.5×10-5m2/s。求解方程(3)、(4)，得到主梁周邊的壓力場和速度場，計算作用在主梁上的氣動力，這可以通過FLUENT直接完成。然后提取升力和力矩賦值給振動方程(1)、(2)，運用Newmark方法求解主梁振動方程。再將速度傳給主梁，通過FLUENT的動網格技術使主梁運動。然后開始下一步的流固耦合循環計算。以上流固耦合計算過程通過對FLUENT 二次開發完成。豎彎振動方程的Newmark方法的算法如下[9]。

1) 計算剛度k，質量m和阻尼矩c。

2.2.2 有效穗。考察結果詳見表5，分析可知施用磷肥的小麥有效穗平均為45.1萬/畝，比未施用磷肥處理有效穗39萬多6.1萬/畝，施用磷肥增加了小麥的有效穗。

3) 選擇步長Δt，參數γ和β，并計算下列有關參數。

4) 形成剛度

對于每個時間步長的計算步驟如下。

①計算t+Δt時刻的荷載

③計算t+Δt時刻的加速度、速度

同理可求扭轉振動方程的解。

2 數值仿真計算模型

數值仿真計算模擬節段模型風洞實驗，但按照實際橋梁尺寸建模，克服了節段風洞實驗模型由于縮小尺寸而帶來的雷諾數效應誤差，主梁取56 m長的節間。計算參數如表1所示[9]，以最不利的第一階正對稱豎彎和第一階正對稱扭轉模態組合計算，結構阻尼比取0.005，節段風洞試驗模型見圖4。

表1 計算參數

圖4 節段風洞試驗模型

數值仿真計算模型見圖5，數值仿真計算區域見圖6，長寬高分別為420 m、72 m和370 m。計算網格劃分見圖7。

圖5 數值計算模型

圖6 數值計算區域

圖7 數值計算網格

如圖7所示，風向從左至右，左側設定為速度入口，右側設定為自由出流。上下邊界為無滑移固壁邊界，數值計算中，靠近主梁網格加密,遠離主梁網格逐漸稀疏。網格總數約100萬，采用有限體積法求解，其中對流項采用中心差分格式，壓力和速度的耦合采用SMPLEC算法，計算采用LES湍流模型，時間步長取0.005 s。由于流固耦合計算費時，再加上網格數量大，因此本文采用了并行計算技術，運用多核工作站進行計算。

3 數值仿真計算

旋渦脫落圖見圖8。

圖8 旋渦脫落圖

由圖8可知，桁架主梁構造復雜，由于弦桿的存在，旋渦脫落表現為復雜的形態，若簡化為平面模型計算，將與實際結構存在較大差別。0°風攻角時，主梁位移時程曲線見圖9、圖10。

圖9 主梁豎向位移時程曲線(0°風攻角)

圖10 主梁扭轉位移時程曲線(0°風攻角)

由圖9和圖10可見，當風速為99 m/s時，主梁豎向位移幅值和扭轉位移幅值基本保持恒定；當風速為101 m/s時，主梁豎向位移幅值和扭轉位移幅值隨時間逐漸增大，呈現發散振動的趨勢。以開始發散振動時的風速作為顫振臨界風速，可知0°風攻角時，顫振臨界風速為99～101 m/s。同理可以計算得到+3°風攻角時，顫振臨界風速為78～80 m/s，顫振臨界風速見表2。

表2 顫振臨界風速

由表2可知，數值仿真計算值和節段模型風洞實驗值基本吻合，顫振臨界風速大于顫振檢驗風速60.4 m/s,滿足顫振穩定性要求。

4 結論

1) 桁架主梁構造復雜，由于弦桿的存在，旋渦脫落表現為復雜的形態，若簡化為平面模型，將與實際結構存在較大差別，難以反映真實的流場情況。

2) 本文數值仿真計算結果和風洞試驗結果基本吻合，說明了方法的有效性和準確性。本文模擬節段模型風洞實驗，同時可以按照實際橋梁建模，能夠克服節段風洞實驗模型由于縮小尺寸而帶來的雷諾數效應誤差。

3) 本方法可以建立結構的三維幾何模型進行抗風計算，適用于任何橋梁主梁結構，為橋梁結構抗風計算提供了新的思路。

[1] LARSEN A，WALTHER J H. Aeroelastic analysis of bridge girder sections based on discrete vortex simulations[J].Journal of Wind Engineering and Industrial Aerodynamic，1997，67/68：253-265.

[2] LARSEN.Advances in aeroelastic analyses of suspension and cable-stayed bridges[J]. Journal of Wind Engineering and Industrial Aerodynamics，1998,74：73-90.

[3] FRANDSEN J B. Numerical bridge deck studies using finite elements.Part I:flutter[J].Journal of Fluids and Structures,2004,19(2):171-191.

[4] MOLHOLM M, JOHANNES H, RASMUSEN J T, LARSEN A, et al.On estimating the aerodynamic admittance of bridge sections by a mesh-free vortex method [J].Journal of Wind Engineering and Industrial Aerodynamics,2015,146:117-127.

[5] 陳艾榮，艾輝林.計算橋梁空氣動力學：大渦模擬[M]北京：人民交通出版社，2010.

[6] 白樺，李宇，李加武，等.鋼桁架懸索橋顫振穩定性能研究[J].振動與沖擊，2013(4):90-95.

[7] 鄧育林,郭慶康,何雄君,等.考慮流固耦合作用的橋梁深水矩形空心高墩振動特性分析[J].武漢理工大學學報(交通科學與工程版),2016,40(6):968-972.

[8] 劉慶寬，盛永青，馬文勇，等.小寬高比鋼桁架懸索橋顫振穩定氣動措施的試驗研究[J].實驗流體力學，2012(3)：26-30.

[9] 王元漢，李麗娟，李銀平.有限元法基礎與程序設計[M].廣州:華南理工大學出版社，2002.