基于CFD 流線型橋梁斷面二維渦激振動(dòng)分析

宋紅紅，楊 剛，姜亞麗

（大連海事大學(xué) 交通運(yùn)輸工程學(xué)院，遼寧 大連 116026）

0 引言

渦激振動(dòng)是由風(fēng)荷載作用在橋梁斷面氣流分離出旋渦及脫落產(chǎn)生，當(dāng)旋渦脫落頻率與橋梁自振頻率相等時(shí)易產(chǎn)生渦激共振現(xiàn)象.流線型箱梁具有較好的抗顫抖振功能，但是大跨度的流線型箱梁在較低風(fēng)速下易產(chǎn)生渦激共振現(xiàn)象[1].渦激共振現(xiàn)象雖然不會(huì)像顫振和馳振那樣對(duì)橋梁產(chǎn)生動(dòng)力失穩(wěn)破壞，但是由于渦激振動(dòng)發(fā)生的風(fēng)速低，頻率高，長時(shí)間的渦激共振現(xiàn)象會(huì)造成橋梁構(gòu)件疲勞損壞，影響行車舒適性和行車安全[2].近幾年國內(nèi)外橋梁渦激振動(dòng)現(xiàn)象并不少見.如英國Kessock 大橋因渦激振動(dòng)引起110 mm 的豎向位移；2012 年，江西九江長江公鐵兩用鋼拱橋發(fā)生吊桿渦激共振現(xiàn)象；2013年，上海楊浦大橋發(fā)生拉索共振、風(fēng)雨振現(xiàn)象[3]；2020 年，武漢鸚鵡洲大橋、廣州虎門大橋都因渦激振動(dòng)作用產(chǎn)生較大起伏，嚴(yán)重影響行車安全.

橋梁渦激振動(dòng)的研究方法主要有風(fēng)洞實(shí)驗(yàn)、數(shù)值模擬、現(xiàn)場實(shí)測(cè)、理論分析.通過這幾種方法對(duì)渦激振動(dòng)的產(chǎn)生機(jī)理、渦激振動(dòng)的影響因素、渦激振動(dòng)的計(jì)算方法、渦振控制進(jìn)行研究[4].基于CFD（Computational Fluid Dynamics）建立橋梁數(shù)值風(fēng)洞技術(shù)已經(jīng)逐步成熟[5-10].文獻(xiàn)[1]中采用數(shù)值風(fēng)洞的方式，避免了實(shí)際風(fēng)洞尺寸的限制，研究了渦激振動(dòng)阻塞比效應(yīng).文獻(xiàn)[11]中發(fā)現(xiàn)了高寬比較大的流線型箱梁斷面渦振的發(fā)生機(jī)制.文獻(xiàn)[12]中將高寬比為4 的矩形斷面渦激振動(dòng)響應(yīng)結(jié)果與風(fēng)洞試驗(yàn)比較，結(jié)果吻合較好，證明采用CFD 技術(shù)進(jìn)行渦激振動(dòng)研究的可行性.還有學(xué)者將CFD 與風(fēng)洞試驗(yàn)相結(jié)合提出一種新的研究方法[13-14]，通過這種方法對(duì)半開口分離雙箱梁的渦激性能及氣動(dòng)措施進(jìn)行了探討[2].盡管橋梁渦激振動(dòng)方面的研究已經(jīng)取得一些重要成果，但是由于不同橋梁產(chǎn)生渦激振動(dòng)的機(jī)理不同，影響渦激振動(dòng)的因素繁多，因此，對(duì)于受風(fēng)環(huán)境影響嚴(yán)重的橋梁有必要分析其渦激振動(dòng)的機(jī)制，找到高效模擬風(fēng)環(huán)境和進(jìn)行渦激振動(dòng)模擬的求解方法，為預(yù)防橋梁渦激共振和渦激振動(dòng)應(yīng)急處理打下基礎(chǔ).

大連長山大橋?yàn)榇罂缍阮A(yù)應(yīng)力混凝土矮塔斜拉橋，主梁采用流線型扁平箱梁，主橋長度1.79 km，主跨長540 m，橫跨長山海峽，由于其跨徑大和近海的地理位置使其受風(fēng)荷載的影響非常明顯.此外，在大橋的設(shè)計(jì)過程和現(xiàn)有的研究中并沒有考慮渦激振動(dòng)的影響.因此，對(duì)大連長山大橋進(jìn)行渦激振動(dòng)分析非常有必要.采用CFD 技術(shù)對(duì)長山大橋主梁截面進(jìn)行二維渦激振動(dòng)研究，探究其脫渦過程及渦激共振的可能性，為評(píng)價(jià)本橋的抗風(fēng)性能提供一定的參考依據(jù).

1 數(shù)值方法與計(jì)算模型

1.1 渦激振動(dòng)計(jì)算方法

將橋梁截面渦激振動(dòng)模型簡化為豎向單自由度彈簧-質(zhì)量-阻尼系統(tǒng)，見圖1.

圖1 單自由度彈性支承柱體結(jié)構(gòu)渦激振動(dòng)模型Fig.1 viv model of a two-degree-of-freedom elastic supported column structure

單自由度彈性支撐體結(jié)構(gòu)渦激振動(dòng)控制微分方程為

式中，y(t)為t時(shí)刻的振動(dòng)位移，m；y˙(t) 為t時(shí)刻的振動(dòng)速度，m/s；˙y˙(t) 為t時(shí)刻的振動(dòng)加速度，m/s2；m為單位長度柱體的質(zhì)量為彈性支承柱體的振動(dòng)系統(tǒng)的阻尼，N/(m·s-1)；ξ為彈性支撐體系阻尼比；為彈性支承柱體的振動(dòng)系統(tǒng)的剛度，N/m；ωy為彈性支承柱體的振動(dòng)系統(tǒng)的振動(dòng)圓頻，Hz；Fy(t)為t時(shí)刻單位長度柱體所受到的豎向力.

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

式中，ρ為空氣密度，取1.225 kg/m3；μ為動(dòng)力粘性系數(shù)，取1.789 4×10-5kg/m·s；u、p分別為流場速度和壓力的時(shí)間平均值，為雷諾應(yīng)力.

1.2 計(jì)算模型

（1）二維流場及網(wǎng)格劃分

以長山大橋典型主梁截面為研究對(duì)象，流場計(jì)算域布置及橋梁截面分別見圖2、圖3.計(jì)算域采用圓形流場，可避免漩渦打到邊界上再折射回來，并保證流體在流場內(nèi)充分繞流.計(jì)算域邊界距離截面形心一般為20 倍的特征長度，因此設(shè)置圓形流場外緣距截面形心為20 倍的截面寬度；入口的邊界條件設(shè)置為速度入口（Velocity inlet），通過給定初始速度來模擬來流風(fēng)速，出口邊界設(shè)置為壓力出口（Pressure outlet），將其與大氣壓的壓差設(shè)置為0 來模擬出口邊界與外界相通.

圖2 計(jì)算域及邊界條件布置Fig. 2 clculation field and boundary condition layout

圖3 主梁截面Fig. 3 schematic of main beam section

由于橋梁截面形狀相對(duì)復(fù)雜，該模型采用三角形非結(jié)構(gòu)化網(wǎng)格，見圖4、圖5.網(wǎng)格劃分采用嵌套網(wǎng)格方式，將橋梁截面視為剛體，截面外圍分別為剛體運(yùn)動(dòng)區(qū)域、動(dòng)網(wǎng)格區(qū)域、靜網(wǎng)格區(qū)域.在靠近截面的地方網(wǎng)格密度較大，最小網(wǎng)格面積約為0.052 m2，流域外圍網(wǎng)格密度較小，最大網(wǎng)格面積約為2.03 m2，這樣既能保證計(jì)算精度，又能通過減少網(wǎng)格數(shù)量來提高計(jì)算效率.

圖4 整體網(wǎng)格示意Fig. 4 overall grid diagram

圖5 局部網(wǎng)格示意Fig.5 ocal grid diagram

（2）求解設(shè)置及計(jì)算流程

結(jié)構(gòu)渦激振動(dòng)求解過程采用Fluent 里的SST k-ω湍流模型進(jìn)行模擬.具體求解設(shè)置為：采用SIMPLEC（Semi Implicit Method for Pressure Linked Equation Consistent）算法求解動(dòng)量方程中速度分量和壓力耦合問題；采用PRESTO（Pressure Staggering Option）求解速度分量；對(duì)流項(xiàng)采用QUICK（Quadratic Upwind Interpolation For Convection Kinetics）求解.單自由度彈性支撐結(jié)構(gòu)渦激共振計(jì)算流程見圖6.

圖6 渦激振動(dòng)計(jì)算流程Fig. 6 viv calculation flow chart

進(jìn)行結(jié)構(gòu)渦激共振計(jì)算時(shí)，首先對(duì)結(jié)構(gòu)進(jìn)行靜態(tài)繞流計(jì)算，然后釋放截面某方向的自由度，進(jìn)行瞬態(tài)動(dòng)力分析，得到每個(gè)時(shí)間步內(nèi)結(jié)構(gòu)上的升力、升力矩、壓力分布等，并將該升力代入式（1），采用Runge-kutta 法求解柱體振動(dòng)響應(yīng)，再將結(jié)構(gòu)動(dòng)力響應(yīng)賦給結(jié)構(gòu)周圍的流體，然后通過Fluent 中動(dòng)網(wǎng)格宏DEFFINE_CG_ MOYION 進(jìn)行傳遞，使結(jié)構(gòu)周圍的動(dòng)網(wǎng)格獲得速度，從而更新網(wǎng)格位置.待網(wǎng)格迭代收斂后，整個(gè)流場更新完成，然后進(jìn)行下一個(gè)時(shí)間步的計(jì)算，直到計(jì)算結(jié)束.

（3）動(dòng)力求解方法

通過UDF 將迭代求解方法嵌入到Fluent 中求解式（1），采用Runge-kutta 迭代求解方法，其表達(dá)式為

式中，

采用四階Runge-kutta 對(duì)二階結(jié)構(gòu)振動(dòng)微分方程進(jìn)行求解，首先要把二階微分方程進(jìn)行降階處理，因此可設(shè)為

將式（6）代入式（1）得

即

經(jīng)計(jì)算，可得結(jié)構(gòu)在下一時(shí)間步的位移為

速度為

式（9）～式（10）中，

2 長山大橋主梁截面渦激振動(dòng)分析

2.1 靜態(tài)繞流分析

在靜態(tài)繞流過程中，流場入口采用速度入口的邊界條件：設(shè)置湍流強(qiáng)度為0.5%，湍流粘性比為10%，風(fēng)速為10 m/s，即折算風(fēng)速為

式中，V為來流風(fēng)速，m/s；fh為結(jié)構(gòu)自振頻率，Hz；D為截面高度，m；流場出口采用壓力出口的邊界條件.設(shè)計(jì)算殘差為5×10-4，步數(shù)為5 000，得到靜態(tài)繞流三分力時(shí)程，見圖7.從圖7 中可以看出，在計(jì)算2 500 步后三分力系數(shù)呈周期性變化，表征流場達(dá)到穩(wěn)定狀態(tài).

圖7 三分力系數(shù)時(shí)程Fig.7 time history of three component force coefficient

2.2 渦激振動(dòng)分析

渦激振動(dòng)分析采用1.2 節(jié)中的求解設(shè)置，選擇與靜態(tài)繞流相同邊界條件，關(guān)閉靜態(tài)分析，開啟瞬態(tài)分析進(jìn)行求解，導(dǎo)入U(xiǎn)DF 程序完成控制方程求解及網(wǎng)格更新.長山大橋模型一階豎彎頻率為0.43Hz，風(fēng)攻角為0°，主梁斷面質(zhì)量為57 980 kg/m，設(shè)置時(shí)間步長為0.01.

通過Fluent 動(dòng)態(tài)捕捉，得到一個(gè)周期內(nèi)的脫渦過程，見圖8.

圖8 脫渦過程Fig. 8 vortex shedding process

由圖8 中可見，旋渦首先產(chǎn)生在橋梁上截面風(fēng)速入口處，然后順著風(fēng)速風(fēng)向慢慢向右推移至橋梁上截面中間位置，隨后轉(zhuǎn)移到橋梁截面尾部，最后旋渦脫落在橋梁截面尾部流場中.從旋渦脫落的整個(gè)過程中可以發(fā)現(xiàn)，在橋梁上截面和下截面及截面尾部渦量較大.圖9 為折算風(fēng)速為12 的渦量圖，對(duì)比發(fā)現(xiàn)不同的風(fēng)速下產(chǎn)生的渦量和脫渦頻率不同，風(fēng)速為12 時(shí)橋梁截面尾部流場產(chǎn)生的旋渦較多.

圖9 折算風(fēng)速12 渦量Fig. 9 translate wind speed 12 vorticity

Strouhal 通過實(shí)驗(yàn)發(fā)現(xiàn)當(dāng)流體流過圓柱后，在尾流中將出現(xiàn)交替脫落的旋渦，并且旋渦脫落的頻率、風(fēng)速及圓柱體直徑之間存在以下關(guān)系

式中，fs為旋渦脫落頻率，Hz；D為圓柱體直徑，mm；U為來流速度，m/s2；St為Strouhal 數(shù).后經(jīng)學(xué)者研究發(fā)現(xiàn)，式（13）也可用于橋梁，其中D為橋梁迎風(fēng)面主梁高度.

按照以上方法，分別對(duì)折算風(fēng)速4、6、8、10、12、14、16、18、20、22 這10 種風(fēng)速下進(jìn)行渦激振動(dòng)計(jì)算.限于篇幅，僅給出8、10、12 風(fēng)速下升力時(shí)程和升力頻譜分析結(jié)果，見圖10～圖12.從圖10～圖12 可以看出，隨著風(fēng)速的增大脫渦頻率逐漸增大，符合式（13）中風(fēng)速與脫渦頻率呈線性增加的關(guān)系.

圖10 折算風(fēng)速8 渦激振動(dòng)響應(yīng)Fig. 10 translate wind speed 8 viv response

圖11 折算風(fēng)速10 渦激振動(dòng)響應(yīng)Fig. 11 Translate wind speed 10 viv response

圖12 折算風(fēng)速12 渦激振動(dòng)Fig. 12 translate wind speed 12 viv response

圖13為豎向位移與折算風(fēng)速關(guān)系，從圖中可以看出，渦振鎖定折算風(fēng)速區(qū)間為8～20，引起的豎向無量綱位移為0.023，此時(shí)脫渦頻率為1.43 Hz，與本橋豎向自振頻率0.43 Hz 相差較大，不易出現(xiàn)渦激共振現(xiàn)象.

圖13 位移隨風(fēng)速的變化Fig.13 displacement changes with wind speed

3 結(jié)論

基于Fluent 動(dòng)網(wǎng)格技術(shù)對(duì)長山大橋典型主梁截面進(jìn)行二維流場渦激振動(dòng)模擬，得到以下結(jié)論.

（1）采用“橋梁截面+運(yùn)動(dòng)區(qū)域+動(dòng)網(wǎng)格區(qū)域+靜止網(wǎng)格區(qū)域+流場邊界”的二維流場建模方式能夠有效防止負(fù)網(wǎng)格的產(chǎn)生，并且能夠減少網(wǎng)格數(shù)量，提高計(jì)算效率.

（2）對(duì)橋梁截面進(jìn)行了折算風(fēng)速為4、6、8、10、12、14、16、18、20、22 這10 種風(fēng)速下單自由度渦激振動(dòng)模擬，分析得到易引起渦激振動(dòng)的折算風(fēng)速為8～20 之間，豎向無量綱振幅最大為0.023，此時(shí)的脫渦頻率為1.43 Hz，與本橋豎向自振頻率相差較大，不易出現(xiàn)渦激共振現(xiàn)象.

（3）流線型箱梁的漩渦產(chǎn)生脫落過程為：產(chǎn)生于橋梁截面入風(fēng)端，經(jīng)橋梁上截面，移至截面尾部，最后脫落在截面尾部尾流中.因此，今后在改變流線型箱梁前端的截面形狀時(shí)，如大型風(fēng)屏障或者大型貨車停靠時(shí)，應(yīng)注意橋梁的渦激振動(dòng)情況，避免發(fā)生共振.