寬高比為4的矩形斷面渦激振動響應(yīng)數(shù)值模擬

劉志文，周 帥，陳政清

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

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

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

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

1 計算模型

1.1 控制方程

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

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

柱體振動控制方程為：

柱體周圍流場采用Fluent6.2中的k－ωSST湍流模型進行計算。具體求解設(shè)置為：采用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，計算進行到充分繞流為止；然后將柱體突然釋放，進行柱體渦激振動響應(yīng)計算。進行渦激振動響應(yīng)計算的每個時間步內(nèi)，先求解流體控制方程，得到速度場、壓力場以及作用在柱體上的升力FL（t），通過嵌入Fluent中的UDF程序來提取作用于柱體上的升力，并將該升力代入柱體振動方程（3）的右端，采用Newmark－β法求解柱體振動響應(yīng)，并將柱體速度賦予柱體周圍隨柱體一起做剛體運動的隨動流體，然后將柱體的速度通過Fluent中的動網(wǎng)格宏DEFINE_CG_MOTION進行傳遞，使柱體周圍的動網(wǎng)格獲得速度，從而更新動網(wǎng)格位置，待網(wǎng)格迭代收斂后，整個流場更新完成從而進行下一個時間步的計算，直到計算結(jié)束為止。渦激振動計算時間步設(shè)置為0.005 s，計算殘差控制在2 ×10－3。

1.2 計算區(qū)域與網(wǎng)格劃分

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

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

圖1 結(jié)構(gòu)渦激振動模型Fig.1 Vortex-induced vibration model of structure

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

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

圖3 矩形斷面網(wǎng)格劃分Fig.3 Computational grid of the rectangular section

2 數(shù)值模擬結(jié)果

2.1 靜止矩形斷面繞流結(jié)果

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

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

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

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

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

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

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

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

3 結(jié)論

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

（1）采用“剛性運動區(qū)域+動網(wǎng)格區(qū)域+靜止網(wǎng)格區(qū)域”的思路建立網(wǎng)格，可以有效解決結(jié)構(gòu)渦激振動響應(yīng)計算中可能產(chǎn)生的網(wǎng)格畸變、負體積等問題；

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

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

［1］鮮 榮，廖海黎，李明水.大跨主梁沿跨向渦振Scanlan非線性模型應(yīng)用［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］曹豐產(chǎn)，項海帆.圓柱非定常繞流及渦激振動的數(shù)值計算［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］方平治，顧 明.高雷諾數(shù)條件下二維方柱渦激振動的數(shù)值模擬［J］.同濟大學學報（自然科學版），2008，36（2）：161－165.

［9］徐 楓，歐進萍.方柱非定常繞流與渦激振動的數(shù)值模擬［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）.