寬高比為4的矩形斷面渦激振動響應數值模擬

劉志文，周 帥，陳政清

（湖南大學 風工程試驗研究中心，長沙 410082）

渦激振動是由于結構尾流區旋渦脫落所引起的一種振動現象，具有強迫、自激、限幅等特點，是一種典型的非線性振動現象。實際工程結構如橋梁結構、天線、索結構、熱交換管及鉆井平臺立管等都有可能發生渦激共振現象，若設計不當則會在氣流作用下產生大幅渦激振動，從而引起結構疲勞或影響結構的使用性能，因此必須予以重視。

鈍體繞流和渦激振動在理論研究和工程實際中都有著很重的意義，現有的研究成果大多是通過試驗得到的［1，2］，鈍體結構渦激振動的數值模擬研究相對較少。Nomura［3，4］采用基于任意拉格朗日 － 歐拉法計算了雷諾數分別為Re=100、Re=140時圓柱和H型柱體的渦激振動響應，該方法假定柱體和流體交界面網格點的運動速度和柱體表面的運動速度相同，在遠離柱體的邊界上假定網格點固定。該方法最大的缺點是柱體的運動容易引起計算網格的畸變。曹豐產等［5］采用動網格法分析圓柱和流體的耦合作用，采用二階投影法和多重網格法求解N－S方程，柱體振動方程采用Newmark－β法求解，實現了圓柱渦激振動響應的數值模擬。Sun等［6］采用LES模型和分塊迭代方法對雷諾數為Re=200時圓柱渦激振動響應進行了數值模擬。李廣望［7］采用任意拉格朗日－歐拉法求解圓柱的渦激振動，結合分塊耦合法提高圓柱渦激振動計算效率，但仍對大振幅條件下的網格畸變和纏繞等不能很好解決。方平治、顧明［8］針對二維方柱在高雷諾數條件下的渦激振動進行了數值模擬，計算得到了方柱的渦激振動“鎖定”和“拍”現象。徐楓、歐進萍［9］基于動網格技術和滑移網格技術，對二維方柱渦激振動響應進行了數值模擬。Pan等［10］采用k－ωSST湍流模型對低質量－阻尼比參數的二維圓柱渦激振動響應進行了數值模擬，計算得到的流體力、圓柱振動響應以及圓柱尾流渦結構均與試驗結果吻合良好。值得注意的是，氣流在流經圓柱、方柱時一般會發生“分離”現象，而很少出現“再附”現象。在實際橋梁結構工程中，氣流流經主梁斷面后一般會同時出現“分離”和“再附”現象。對出現“分離”和“再附”現象的結構進行渦激振動數值模擬對工程實際而言具有重要的意義。

為此本文采用“剛性運動區域+動網格區域+靜止網格區域”的思路建立網格［11］，來解決結構振動過程中可能出現的網格畸變與負體積的問題，并將求解結構振動響應的Newmark－β通過UDF嵌入到Fluent中，對寬高比為4的矩形斷面結構渦激振動響應進行了數值模擬，該方法區別于文獻［9］中的滑移網格技術。

1 計算模型

1.1 控制方程

粘性不可壓縮流體的控制方程有質量守恒方程（連續性方程）、動量守恒方程（N－S方程）和能量守恒方程。在進行渦激振動響應計算時不考慮能量方程，在直角坐標系下，基于雷諾平均的連續性方程和N－S方程分別為：

式中，i，j=1，2；ρ為空氣密度，即ρ=1.225 kg/m3，μ為動力粘性系數，即μ=1.789 4 ×10－5kg/m·s。

柱體振動控制方程為：

柱體周圍流場采用Fluent6.2中的k－ωSST湍流模型進行計算。具體求解設置為：采用SIMPLEC（Semi implicit method for pressure linked equation consistent）算法求解動量方程中速度分量和壓力的耦合問題；采用PRESTO（Pressure staggering option）求解速度分量；對流項采用QUICK（Quadratic upwind interpolation for convection kinetics）求解。在進行柱體渦激振動計算時，首先針對靜止柱體繞流進行計算，計算物理時間步長為0.005 s，計算殘差控制在 5 ×10－4，計算進行到充分繞流為止；然后將柱體突然釋放，進行柱體渦激振動響應計算。進行渦激振動響應計算的每個時間步內，先求解流體控制方程，得到速度場、壓力場以及作用在柱體上的升力FL（t），通過嵌入Fluent中的UDF程序來提取作用于柱體上的升力，并將該升力代入柱體振動方程（3）的右端，采用Newmark－β法求解柱體振動響應，并將柱體速度賦予柱體周圍隨柱體一起做剛體運動的隨動流體，然后將柱體的速度通過Fluent中的動網格宏DEFINE_CG_MOTION進行傳遞，使柱體周圍的動網格獲得速度，從而更新動網格位置，待網格迭代收斂后，整個流場更新完成從而進行下一個時間步的計算，直到計算結束為止。渦激振動計算時間步設置為0.005 s，計算殘差控制在2 ×10－3。

1.2 計算區域與網格劃分

考慮到工程實際中，氣流流經橋梁主梁斷面時一般會發生“分離”和“再附”現象，這區別于氣流流經方柱和圓柱時的流動現象；同時為驗證本文方法的精度，選擇文獻［12］中寬高比為4的矩形斷面進行渦激振動響應數值模擬。矩形斷面寬為B=0.3 m，矩形斷面高為D=0.075 m，將矩形斷面豎向渦激振動系統簡化為沿豎向振動的彈簧－質量－阻尼系統，如圖1所示。渦激振動計算參數取自文獻［12］，自振頻率fn=3.905 Hz，單位長度矩形柱體質量為1.362 kg/m，阻尼比為0.177%，對應的折算速度Vred=V∞/fnD（其中V∞為來流風速（m/s），fn為結構振動頻率（Hz），D為矩形斷面高（m））分別為 4、6、8、9、10、12、14 及 16，對應的雷諾數（以矩形斷面高度D為參考尺寸）范圍為6 015～24 060。

計算區域為：上游邊界距矩形斷面為10B，上下側邊界距矩形斷面中心距離為10B，下游邊界距矩形斷面為20B。邊界條件設置如下：上游為速度入口，即Velocity-inlet邊界條件，湍流強度為Iu=0.5%，湍流粘性比為 10%；出口邊界采用Outflow邊界條件，即完全發展的出流邊界條件；上下兩側采用對稱邊界條件，即Symmetry邊界條件；矩形斷面四周采用無滑移邊界條件，即No-Slip Wall邊界條件。

圖1 結構渦激振動模型Fig.1 Vortex-induced vibration model of structure

圖2 矩形斷面網格分塊劃分及邊界條件Fig.2 Grid blocks and boundary conditions of the rectangular section

為了解決動網格模擬中網格的畸變、纏繞等問題，采用“剛性運動區域+動網格區域+靜止網格區域”的方法建立矩形斷面渦激振動網格劃分區域。在剛性運動區域，采用結構化網格對結構周圍的區域網格進行加密，這部分網格在結構振動過程中隨結構一起運動，不需要重新劃分網格。動網格區域采用三角形網格進行劃分，在結構振動過程中，該部分網格進行重新劃分，由于該部分網格遠離結構，故網格間距較大，網格不容易出現負體積，從而有效解決結構大幅渦激振動可能引起的網格畸變和負體積現象。為了提高計算效率，最外側網格采用靜止的結構化網格進行劃分。在劃分網格時，首先對矩形斷面附近建立邊界層網格，矩形斷面附近第一層網格中心至矩形斷面表面的距離為δ=0.001B，然后分塊進行網格劃分，整個計算區域網格總數為75 848。圖3分別給出了矩形斷面周圍流場區域網格總體圖及局部示意圖。

圖3 矩形斷面網格劃分Fig.3 Computational grid of the rectangular section

2 數值模擬結果

2.1 靜止矩形斷面繞流結果

限于篇幅，圖4僅給出了折算風速為Vred=V∞/fnH=9.0時靜止矩形斷面升力系數時程曲線及升力系數FFT變換幅值頻譜曲線。從圖4中可以看出氣流流經該矩形斷面時所產生的升力系數出現明顯的振幅不變的振蕩現象，表明氣流流經該斷面時所產生的漩渦已達到穩定狀態。表1給出了寬高比為4的矩形斷面阻力系數、升力系數根方差以及斯托羅哈數計算結果。從表1中可以看出，本文阻力系數、斯托羅哈數結果與文獻［14］試驗結果比較接近，但升力系數根方差比文獻［13］的結果要偏小。

圖4 靜止矩形斷面繞流數值模擬結果Fig.4 Numerical simulation results of flow past the stationary rectangular section

表1 寬高比為4的矩形斷面繞流數值模擬與試驗結果Tab.1 Numerical and experimental results of flow past the fixed rectangular section with aspect ratio 4

2.2 彈性支承矩形斷面渦激振動模擬結果

根據文獻［12］的寬高比為4的矩形斷面渦激振動振幅與速度圖（即A－V圖），分別對折算風速為4、6、8、9、10、12、14 及 16 時進行了渦激振動響應計算。限于篇幅，圖5僅給出折算風速Vred=V/fnD分別為4、10及16時，矩形斷面升力系數CL、渦激振動響應時程曲線。從圖5中可以看出，不同折算風速條件下，彈性支承矩形斷面的渦激振動響應振幅是不同的。

圖6給出了折算風速Vred=V/fnD分別為4、10及16時矩形斷面升力系數CL的FFT變化幅值頻譜曲線。從圖6中可以看出，當折算風速為Vred=V/fnD=4時，矩形斷面升力系數CL的FFT變化幅值頻譜對應的卓越頻率為fs=1.532 6 Hz，低于矩形斷面結構振動頻率，未發生鎖定現象；當折算風速為Vred=V/fnD=10時，矩形斷面升力系數CL的FFT變化幅值頻譜對應的卓越頻率為fs=3.759 1 Hz，與矩形斷面的結構振動頻率接近，即出現“鎖定”現象；隨著折算風速的增加，矩形斷面的漩渦脫落頻率逐漸增加，“鎖定”現象消失。

圖7 渦激振動幅值及漩渦脫落頻率比隨折算風速變化Fig.7 VIV amplitude and ratio of shedding frequency to natural frequency versus reduced wind velocity

圖7分別給出了矩形斷面渦激振動振幅隨折算風速的變化曲線、漩渦脫落頻率與結構振動頻率之比隨折算風速的變化曲線。從圖7（a）中可以看出，寬高比為4的矩形斷面渦激振動響應數值模擬結果（如鎖定區間、最大振幅）均與文獻［10］中的風洞試驗結果吻合較好。由圖7（b）可知，當折算風速為Vred=V/fnD=8～12時，矩形斷面漩渦脫落頻率與結構振動頻率之比接近1，即發生了“鎖定”現象。

3 結論

本文針對寬高比為4的矩形斷面的渦激振動進行了數值模擬，并與已有試驗結果進行了比較，主要結論如下：

（1）采用“剛性運動區域+動網格區域+靜止網格區域”的思路建立網格，可以有效解決結構渦激振動響應計算中可能產生的網格畸變、負體積等問題；

（2）通過對寬高比為4的矩形斷面的渦激振動數值模擬可知，采用本文所建立的方法進行具有“分離”、“再附”現象的鈍體斷面的渦激振動是可行的。

需要說明的是，本文僅僅針對形狀較為簡單的鈍體斷面進行了渦激振動數值模擬，對于斷面形狀較為復雜的鈍體斷面（如考慮防撞護欄、檢修車軌道等附屬設施的橋梁主梁斷面），則需要進一步的研究與試驗驗證。

［1］鮮 榮，廖海黎，李明水.大跨主梁沿跨向渦振Scanlan非線性模型應用［J］.振動與沖擊，2009，28（4）：54－58.

［2］歐陽克儉，陳政清，韓 艷，等.橋面中央開口懸索橋渦激共振與制渦試驗研究［J］.振動與沖擊，2009，28（7）：199－202.

［3］Nomura T. Finite elementanalysis ofvortex-induced vibrations of bluff cylinders［J］.Journal of Wind Engineering and Industrial Aerodynamics，1993，46，47：587 －594.

［4］ Nomura T.A numerical study on vortex-induced oscillations of bluff cylinders［J］.Journal of Wind Engineering and Industrial Aerodynamics，1993，50：75－83.

［5］曹豐產，項海帆.圓柱非定常繞流及渦激振動的數值計算［J］.水動力學研究與進展（A 輯），2001，16（1）：111－118.

［6］李廣望，任安祿，陳文曲.ALE方法求解圓柱的渦激振動［J］.空氣動力學報，2004，22（3）：283－288.

［7］ Sun D，Owen J S，Wright N G，et al.Fluid-structure interaction of prismatic line-like structures using LES and Block-Iterative coupling［C］//Aerodynamics Laboratory National Research Council Canada.5th International Colloquium on Bluff Body Aerodynamics and Applications，Ottawa，Canada：2004.133－136.

［8］方平治，顧 明.高雷諾數條件下二維方柱渦激振動的數值模擬［J］.同濟大學學報（自然科學版），2008，36（2）：161－165.

［9］徐 楓，歐進萍.方柱非定常繞流與渦激振動的數值模擬［J］.東南大學學報（自然科學版），2005，35（S1）：35－39.

［10］ Pan Z Y，Cui W C，Miao Q M.Numerical simulation of vortex-induced vibration of a circular cylinder at low massdamping using RANS code［J］.Journal of Fluids and structures，2007，23：23 －37.

［11］ Huang L，Liao H L，Wang B，et al.Numerical simulation for aerodynamic derivatives of bridge deck［J］.Simulation Modelling Practice and Theory，2009，17：719－729.

［12］ Masaru M，Tomomi Y，Hitoshi T，et al.Vortex-induced vibration and its effect on torsional flutter instability in the case of B/D=4 rectangular cylinder［J］.Journal of wind Engineering and Industrial Aerodynamics，2008，96：971－983.

［13］ Sun D，Owen J S，Wright N G.Application of thek-ωturbulence model for a wind-induced vibration study of 2D bluff bodies［J］.Journal of Wind Engineering and Industrial Aerodynamics，2009，7，77 －87.

［14］ Nakaguchi H，Hashimoto K，Muto S，An experimental study on aerodynamic drag of rectangular cylinders［J］.J.Japan Soc.Aeronaut.Space Sci，1968，16：1 －5（in Japanese）.