薄平板斷面氣動自激力振幅效應數值模擬

劉志文，張瑞林，陳政清，2

（1.湖南大學風工程與橋梁工程湖南省重點實驗室，湖南長沙 410082；2.湖南大學土木工程學院，湖南長沙 410082）

1940年，舊塔科馬橋在19 m/s左右的風速下發生大幅振動，扭轉振幅大約為35°，并在長時間振動下因吊桿斷裂導致垮塌[1].此后，橋梁結構氣動穩定性問題受到廣泛關注，Scanlan等[2]將航空領域的顫振導數概念加以推廣，建立了適用于橋梁主梁斷面的線性顫振理論，為現代大跨度橋梁抗風設計提供指導.隨著橋梁跨度的進一步增加，大跨橋梁主梁在強風作用下可能會出現大幅振動現象，傳統的Scanlan線性顫振理論難以適用于橋梁主梁大幅振動情況.國內外許多學者開始關注振幅對橋梁主梁斷面的顫振性能影響.

橋梁主梁斷面氣動力在大振幅條件下存在高次諧波分量.陳政清等[3]通過強迫振動風洞試驗研究發現鈍體橋梁主梁斷面存在明顯的高次諧波分量.Diana等[4-5]通過墨西拿海峽大橋主梁斷面風洞試驗研究，發現氣動力存在高次諧波分量和氣動力的遲滯效應.王騎[6]通過風洞試驗研究發現大振幅下薄翼斷面與流線型箱梁斷面氣動力存在顯著的高次諧波分量.唐煜[7]采用計算流體動力學方法（Computational Fluid Dynamics，CFD）進行了薄平板斷面與流線型箱梁斷面強迫振動研究，結果表明：大振幅下氣動力存在高次諧波分量，并且振幅越大，高次諧波分量越明顯.熊龍等[8]、Lin等[9]分別進行了流線型箱梁斷面與B/H=5矩形斷面（B、H分別為矩形斷面的寬度與高度）強迫振動風洞試驗研究，結果表明：較大扭轉振幅下流線型箱梁斷面與B/H=5矩形斷面氣動力存在顯著的高次諧波分量.Gao等[10]、Zhang等[11]分別在雙邊肋主梁斷面與流線型箱梁斷面極限環顫振中發現氣動力存在二階與三階分量.

顫振導數是表征橋梁結構斷面氣動力的主要參數，為此部分學者進行了典型斷面顫振導數的振幅效應研究.Noda等[12]通過風洞試驗研究了寬高比B/H=13與150的薄矩形斷面在不同振幅下的顫振導數，研究表明：扭轉振幅對薄矩形斷面顫振導數影響較大.熊龍等[8]、Lin等[9]分別進行了流線型箱梁斷面與B/H=5矩形斷面大振幅的強迫振動風洞試驗研究，結果表明：扭轉振幅對流線型箱梁斷面與B/H=5矩形斷面顫振導數影響較大，其中對的影響最為明顯.Zhang等[11]采用CFD的方法進行了流線型箱梁斷面的強迫振動研究，結果表明：箱梁斷面顫振導數具有顯著的振幅相關性，而對豎向扭轉模態耦合因素不敏感.

綜上所述，國內外許多學者開展了主梁斷面氣動力的振幅效應的試驗與數值模擬研究，已有研究表明：大振幅下氣動力存在高階分量，而且隨振幅的增大顫振導數會發生改變.但目前對橋梁斷面氣動力幅值與遲滯相位隨振幅變化以及非線性氣動力作用下顫振位移響應演變規律的研究較少.

為此，本文選取薄平板斷面作為研究對象[3，7，13]，采用CFD方法進行薄平板斷面氣動自激力的振幅效應研究.首先采用強迫振動數值模擬方法，比較了不同振幅下的薄平板斷面顫振導數與氣動力遲滯相位，分析了薄平板斷面氣動力的頻譜特性，然后采用自由振動數值模擬方法，分析了薄平板斷面的顫振響應演變規律.此外，為驗證本文CFD數值模擬結果的可靠性，將顫振導數與顫振臨界風速結果與理想平板理論解以及已有文獻試驗結果進行了對比.

1 數值計算方法

1.1 流體力學控制方程

目前，在橋梁結構風工程領域，基于雷諾平均法（RANS）的數值模擬方法被廣泛采用[14-15]，其原理是將瞬態方程中的場變量分解為時均值和脈動值：

1.2 求解與計算

本文采用大型商業軟件ANSYS Fluent求解流場，采用剪應力輸運湍流模型（SST k-w）以使流體控制方程封閉，基于有限體積法求解控制方程.空間離散采用二階迎風格式，時間離散采用二階隱式時間積分.基于SIMPLEC算法處理壓力-速度耦合.

數值模擬計算在DELL Precision 7920 Tower工作站進行，采用4個核并行計算，網格數量約為8.57萬，計算效率約為1 h/千步.

1.3 強迫振動與自由振動數值模擬方法

采用大型流體力學分析軟件ANSYS Fluent動網格技術，通過動網格驅動宏（DEFINE-CG-MOTION）實現結構斷面運動.對于結構強迫振動，指定結構斷面相應速度即可實現結構運動狀態更新.對于結構自由振動，通過求解結構振動方程得到結構振動的位移響應，再通過動網格實現結構運動狀態的更新.結構在升力和升力矩作用下的振動方程為：

式中：m為單位長度質量；ξh0為豎向運動阻尼比；wh0為豎向運動圓頻率；I為單位長度慣性矩；ξa0為扭轉運動阻尼比；wa0為扭轉運動圓頻率；Lse、Mse分別為升力與升力矩.結合Fluent動網格技術，通過用戶定義函數（UDF）將Newmark-β 數值算法嵌入Fluent中以求解式（4）（5），計算薄平板斷面自由振動響應.

1.4 薄平板斷面及計算域

薄平板斷面寬度B為0.45 m，高度H為0.001 m，寬高比B/H為450 ∶1.計算域設為40B×20B（B為薄平板斷面寬度）.圖1所示為薄平板斷面形狀以及計算域與邊界條件.計算域與邊界條件如下：計算域左側為速度入口邊界（Velocity inlet），計算域右側為壓力出口邊界（Pressure outlet），計算域上下側為對稱邊界（Symmetry），薄平板斷面為固定壁面邊界（Wall）.

圖1 薄平板斷面形狀、計算域及邊界條件（單位：mm）Fig.1 Thin plate section，computational domain and boundary conditions（unit：mm）

1.5 網格無關性與時間步無關性驗證

以扭轉振幅α0=12°、折算風速U/fB=17.78的強迫振動為例進行網格無關性與時間步無關性驗證.不同網格劃分方案如表1所示，圖2（a）所示為氣動升力與力矩時程計算結果，其中時間步Δt均為0.000 5 s.由圖2（a）可知不同網格方案計算得到的氣動升力與力矩時程曲線幾乎完全一致.對網格2進行時間步無關性驗證，時間步Δt分別為0.000 3 s、0.000 5 s與0.000 8 s，圖2（b）所示為氣動升力與力矩時程計算結果，由圖2（b）可知不同時間步計算得到的氣動升力與力矩時程曲線幾乎完全一致.

表1 薄平板斷面網格參數Tab.1 Grid parameters of the thin plate section

圖2 氣動力時程曲線（α0=12°、U/fB=17.78）Fig.2 Aerodynamic time history curve（α0=12°，U/fB=17.78）

為確保計算結果可靠，同時考慮計算效率，確定網格2作為計算模型，即網格數量為8.57萬，首層網格高度為10-5B，時間步Δt為0.000 5 s，近壁面網格情況如圖3所示.計算結果顯示y+＜1，如圖4所示.

圖3 近壁面網格Fig.3 Near-wall grids

圖4 薄平板斷面y+分布Fig.4 y+distribution of thin plate section

2 振幅對氣動自激力的影響

2.1 不同振幅下顫振導數

顫振導數是表征結構斷面氣動自激力的重要氣動參數，其與結構斷面運動狀態線性組合表示氣動力的線性部分.在顫振導數識別方面有強迫振動方法與自由振動方法，考慮到強迫振動方法識別顫振導數具有重復性好、折算風速范圍廣的優點[16]，本文采用分狀態單自由度強迫振動方法識別顫振導數.

風速U分別為4 m/s、6 m/s、8 m/s、10 m/s、12 m/s、16 m/s和18 m/s；雷諾數Re為1.2×105～5.5×105；扭轉運動振幅α0分別取為1°、2°、4°、6°、8°、10°和12°；豎向運動振幅h0分別取為0.02B、0.08B、0.16B和0.20B；振動頻率f均為2.0 Hz.薄平板斷面線性氣動自激力表達式為：

當薄平板斷面分別作扭轉強迫振動與豎向強迫振動時，對應的扭轉與豎向位移分別為：

式中：α0為扭轉運動振幅；h0為豎向運動振幅.

根據不同振幅條件下薄平板斷面氣動自激力時程，采用最小二乘法進行薄平板斷面顫振導數識別，不同振幅對應的薄平板斷面的顫振導數隨折算風速的變化曲線如圖5所示.

由圖5可知，薄平板斷面在小振幅下（α0=1°、h0=0.02B）的顫振導數均與理想平板顫振導數理論解[17]吻合較好，表明本文強迫振動數值模擬結果具有良好的精度.

圖5 不同振幅下薄平板斷面顫振導數隨折算風速的變化Fig.5 Flutter derivatives of the thin plate section at different amplitudes vs.reduced wind velocity

由圖5（a）～（d）可知，扭轉振幅對薄平板斷面的顫振導數影響較大，其中對顫振導數的影響最為明顯.當扭轉振幅α0=6°，且折算風速較大時，顫振導數與小振幅下（α0=1°）顫振導數開始出現差異；在高折算風速下，隨著振幅增大（α0=10°、12°）顫振導數均出現了由負轉正的變化，這與文獻[12]薄矩形板（B/H=150）的風洞試驗結果一致.當扭轉振幅為α0=6°、8°、10°、12°時，隨折算風速的增大顫振導數與小振幅下（α0=1°）顫振導數的差別逐漸增大.扭轉振幅對顫振導數的影響較小，但在高折算風速下也有一定影響.

由圖5（e）～（h）可知，豎向振幅對薄平板斷面的顫振導數影響總體較小.

2.2 氣動自激力的幅值與遲滯相位

薄平板斷面的線性氣動自激力表示為正弦函數形式，即：

式中：L0、M0為氣動力幅值；φL、φM為氣動力相對于位移的遲滯相位.比較式（6）～（9）與式（10）（11），可以將顫振導數用氣動力幅值與氣動力遲滯相位表示，如式（12）～（19）所示.由式（12）～（19）可知，顫振導數的改變可由L0/α0、M0/α0、L0/h0、M0/h0與φL、φM表征.

根據薄平板斷面強迫振動模擬所得的氣動力，分別對式（10）（11）進行最小二乘擬合，可得到薄平板斷面的氣動力幅值與遲滯相位，具體結果如圖6所示.

由圖6（a）～（f）可知，對于扭轉運動，L0/α0、M0/α0隨振幅變化較小.升力與力矩的遲滯相位的正弦值sin φLα、sin φΜα隨振幅變化較大，其中sin φMα在折算風速U/fB為20.0時，其值在1°振幅時為0.28，在12°振幅時為-0.26，這也引起了顫振導數在高折算風速下，隨著振幅的增大發生了由負轉正的變化.升力與力矩的遲滯相位余弦值cos φLα、cos φΜα隨振幅變化不明顯，且其絕對值大多大于0.94，所以遲滯相位并沒有引起顫振導數發生較大變化.

圖6 氣動力幅值與遲滯相位正弦值與余弦值Fig.6 Amplitude and hysteresis phase sinusoidal and cosine values of aerodynamic forces

由圖6（g）～（l）可知，對于豎向運動，L0/h0、M0/h0隨振幅變化較小.升力與力矩的遲滯相位的正弦值sin φLh、sin φΜh變化亦較小，且其絕對值大多大于0.96，所以顫振導數隨振幅變化較小.升力與力矩的遲滯相位的余弦值cos φLh、cos φΜh在低折算風速下略有變化，高折算風速下變化不明顯，這引起顫振導數在低折算風速下有較小變化.

2.3 氣動自激力頻譜特性

采用快速傅里葉變換方法（Fast Fourier Transform，FFT）對薄平板斷面單自由度強迫振動得到的氣動自激力進行頻譜分析，圖7為折算風速U/fB為17.78的薄平板斷面氣動力頻譜圖.

由圖7（a）（b）可知，當薄平板斷面的扭轉振幅為12°時，升力與力矩均存在明顯高階分量，主要為三階、五階等奇數倍頻分量；由圖7（c）（d）可知，豎向振幅為0.20B時，氣動力同樣存在高階分量，主要為三階分量，與扭轉運動相比，豎向運動引起的薄平板斷面的氣動力高階分量則相對較小.

圖7 氣動力頻譜特性（U/fB=17.78）Fig.7 Spectrum characteristics of aerodynamic forces（U/fB=17.78）

為進一步分析扭轉振幅對薄平板斷面氣動力分量的影響，定義氣動力分量幅值比Ri：

式中：Fi為i階氣動升力或力矩幅值；f為頻率；Ff為頻率f對應的氣動力幅值.由式（20）計算不同扭轉振幅下薄平板斷面氣動力一階、三階與五階分量的幅值比，結果如圖8所示.

圖8 不同扭轉振幅下氣動力分量幅值比（U/fB=17.78）Fig.8 Amplitude ratio of aerodynamic components at different torsional amplitude（U/fB=17.78）

由圖8可知，隨著扭轉振幅的增大，薄平板斷面氣動力一階分量所占比例不斷下降，當扭轉振幅大于8°時，薄平板斷面氣動力出現較為明顯的三階與五階分量.

3 顫振位移響應演變規律

3.1 顫振臨界風速

采用自由振動方法進行薄平板斷面顫振穩定性的直接計算[18-19]，薄平板斷面自由振動計算特征參數采用文獻[1]中的相關參數，如表2所示.

表2 薄平板斷面自由振動計算特征參數Tab.2 Calculating dynamic parameters of free vibration of the thin plate section

顫振臨界狀態附近不同計算風速下薄平板斷面位移時程曲線如圖9所示.薄平板斷面顫振臨界風速與頻率的計算結果及其與頻域理論解如表3所示.由表3可知，CFD計算誤差較小，表明本文自由振動數值模擬結果具有良好的精度.

圖9 位移時程曲線（U=16.20、16.34、16.50 m/s）Fig.9 Displacement time history curve（U=16.20，16.34，16.50 m/s）

表3 薄平板斷面顫振臨界風速與頻率Tab.3 Critical flutter wind velocity and freguency of the thin plate section

3.2 顫振位移響應演變規律

非線性氣動力作用下，結構振動的頻率和阻尼比不是一個常數，而是隨振幅變化，可根據薄平板斷面自由振位移響應時程計算得到：

式中：Q（ti）為第i個位移峰值；ti為第i個位移峰值對應的時間；ξ（t）為瞬時阻尼比；ω（t）瞬時圓頻率；φhα（t）為豎向與扭轉振動位移相位差.另外，瞬時阻尼比為系統阻尼與氣動阻尼兩者共同作用的結果.

圖10為薄平板斷面在風速U為17.0 m/s（U/fa0B=12.49）時自由振動位移響應時程曲線.根據位移響應時程由式（21）～（23）計算得到薄平板斷面豎向與扭轉振動對應的瞬時頻率、瞬時阻尼比以及豎向與扭轉振動位移相位差隨時間變化曲線，如圖11所示.

圖10 位移時程曲線（U=17.0 m/s）Fig.10 Displacement time history curve（U=17.0 m/s）

圖11 瞬時頻率、瞬時阻尼比與豎向扭轉位移相位差（U=17.0 m/s）Fig.11 Instantaneous frequency，instantaneous damping ratio and phase difference between the vertical displacement and torsional displacement（U=17.0 m/s）

由圖11（a）可知，豎向與扭轉瞬時頻率在小振幅（α0＜1.44°）下不隨振幅變化，且兩者為同一數值.隨著振幅的增大，兩者均呈現了減小的趨勢，在減小過程中豎向瞬時頻率略大于扭轉瞬時頻率，換言之，是扭轉瞬時頻率下降略早于豎向瞬時頻率.由于薄平板斷面為彎扭耦合顫振，豎向與扭轉頻率變化趨勢具有一致性，但在變化過程中兩者頻率不是同一個數值.值得說明的是瞬時頻率存在一定跳躍性，這是由時間離散引起的數值誤差.

由圖11（b）可知，豎向與扭轉瞬時阻尼比絕對值隨著振幅增大而不斷增大，說明顫振發散速率越來越大.這與顫振導數隨振幅變化有關，隨著扭轉振幅的增大顫振導數逐漸變小，再變為負值，導致隨著扭轉振幅的增大顫振導數提供的氣動正阻尼逐漸減小，再變為提供氣動負阻尼.薄平板斷面顫振為豎向與扭轉耦合的振動形式，所以豎向與扭轉瞬時阻尼比的絕對值均隨振幅的增大而增大，且豎向與扭轉瞬時阻尼比基本保持一致.

由圖11（c）可知，扭轉振幅α0＞1.44°時，豎向與扭轉位移相位差隨著振幅增大不斷增大.豎向位移相位超前于扭轉相位，又由于豎向瞬時頻率略大于扭轉瞬時頻率，導致豎向與扭轉位移相位差進一步增大.

4 結論

以薄平板為研究對象，采用計算流體動力學方法進行強迫振動與自由振動研究，分析了氣動自激力的振幅效應，得到如下主要結論：

2）當扭轉振幅大于8°時，薄平板斷面氣動力存在較為明顯高次諧波分量，主要為三階與五階分量；豎向振幅引起的薄平板斷面氣動力高階分量成分相對較小.

3）薄平板斷面顫振發散過程中，豎向與扭轉瞬時頻率在小振幅（α0＜1.44°）下基本不隨振幅變化，隨著振幅進一步增大兩者均呈現減小的趨勢；豎向與扭轉瞬時阻尼比絕對值隨著振幅增大而不斷增大；扭轉振幅α0＞1.44°時，豎向與扭轉位移相位差隨著振幅增大不斷增大.