祝志文， 石亞光， 顏 爽

(1. 汕頭大學(xué) 土木與環(huán)境工程系，廣東 汕頭 515063； 2.中冶南方工程技術(shù)有限公司，武漢 430000)

現(xiàn)代橋梁跨度的不斷增大使得結(jié)構(gòu)逐步向輕柔化和低阻尼方向發(fā)展，導(dǎo)致其對(duì)風(fēng)致作用更加敏感。使得其在常遇風(fēng)速下存在發(fā)生渦激振動(dòng)的可能。渦激振動(dòng)是一種具有自激和自限幅特性的風(fēng)致振動(dòng)，已建成的丹麥大帶東橋[1]，日本東京灣大橋[2]，巴西Rio-Niteroi橋[3]和西堠門大橋[4]均發(fā)生過大幅度的渦激振動(dòng)。2020年5月5日廣東虎門大橋多次發(fā)生渦激振動(dòng)，引起社會(huì)廣泛關(guān)注。當(dāng)橋梁跨度大于1 000 m后，結(jié)構(gòu)基頻很小且模態(tài)密集，在常遇風(fēng)速內(nèi)可依次激發(fā)橋梁主梁的多個(gè)模態(tài)，也即可激發(fā)基頻以上的橋梁高階模態(tài)振動(dòng)[5-6]。橋梁渦激振動(dòng)雖不如顫振一樣導(dǎo)致橋梁的災(zāi)難性破壞，但會(huì)影響橋面行車的舒適性和安全性，過大的渦振所產(chǎn)生的高應(yīng)力幅可能會(huì)在橋梁主體結(jié)構(gòu)的焊縫部位產(chǎn)生疲勞破壞。因此，準(zhǔn)確預(yù)測大跨度橋梁的渦激振動(dòng)并采取有效的預(yù)防措施顯得尤為重要。

現(xiàn)階段橋梁渦激振動(dòng)研究主要借助節(jié)段模型風(fēng)洞試驗(yàn)[7-9]，但節(jié)段模型試驗(yàn)獲得的主梁最大渦激振動(dòng)幅值可能與實(shí)橋存在差別，因而可能需要修正，其中的途徑如通過節(jié)段模型結(jié)果與多點(diǎn)彈性支承模型的對(duì)比得到[10-11]?，F(xiàn)有渦激振動(dòng)數(shù)值模擬的研究也主要集中在形狀規(guī)則的斷面上，如圓柱、方柱和橢圓等。Daniels等[12]采用大渦模擬(large eddy simulation,LES)方法結(jié)合動(dòng)網(wǎng)格技術(shù)對(duì)長寬比為4的矩形柱進(jìn)行了渦激振動(dòng)數(shù)值模擬。徐楓等[13]采用二維RANS(Reynolds-averages Navier-Stokes)方法模擬了不同截面形狀彈性支撐柱體的流致振動(dòng)，結(jié)果表明柱體的振動(dòng)以橫向振動(dòng)為主。Seyyed等[14]針對(duì)低雷諾數(shù)下的橢圓渦激振動(dòng)進(jìn)行了二維模擬，并分析了不同攻角和不同橢圓外形對(duì)渦激振動(dòng)的影響。對(duì)橋梁渦激振動(dòng)的計(jì)算流體動(dòng)力學(xué)(computational fluid dynamics，CFD)模擬，特別是包含防撞欄桿等附屬設(shè)施的成橋狀態(tài)斷面開展得非常少，且鮮見針對(duì)帶防撞欄橋梁主梁高階模態(tài)渦激振動(dòng)的CFD研究報(bào)道。

本文以大跨度橋梁較為普遍采用的閉口扁平箱梁為例，采用分區(qū)網(wǎng)格劃分策略解決橋梁振動(dòng)過程中容易出現(xiàn)的網(wǎng)格畸變問題，并將求解橋梁渦激振動(dòng)響應(yīng)的Newmark-β法通過自開發(fā)的UDF程序嵌入到Fluent中，開展了大帶東橋扁平箱梁斷面的高階模態(tài)渦激振動(dòng)模擬，得到了高階模態(tài)渦振鎖定曲線并在此基礎(chǔ)上對(duì)渦激振動(dòng)的影響因素進(jìn)行了分析。

1 基準(zhǔn)模型CFD模擬的數(shù)值實(shí)現(xiàn)

在1999年6月24日召開的第10屆國際風(fēng)工程會(huì)議上，國際風(fēng)工程協(xié)會(huì)確定丹麥大帶東橋主橋加勁梁作為橋梁CFD模擬的基準(zhǔn)模型[15]。該加勁梁為扁平六邊形閉口箱梁，梁寬31 m，高4.4 m，橋面布設(shè)鋼制中央分隔欄和上下游側(cè)欄桿，如圖1所示，其中S.C.代表加勁梁剪切中心，在橋軸線距離底板0.535倍加勁梁梁高位置。需要指出，早期關(guān)于大帶東橋CFD研究的諸多文獻(xiàn)，由于方案的變更導(dǎo)致一些研究的加勁梁氣動(dòng)外形與成橋不一致的情況，包括梁寬和欄桿外形等。需要指出，該橋在4～12 m/s的風(fēng)速下就曾發(fā)生了顯著的多階豎向渦激振動(dòng)，包括高階模態(tài)渦激振動(dòng)。

圖1 大帶東橋主跨加勁梁斷面(m)Fig.1 Stiffening girder cross section of main span of Great Belt East Bridge(m)

1.1 控制方程

計(jì)算風(fēng)工程中，橋梁繞流問題的二維黏性不可壓雷諾時(shí)均Navier-Stokes方程表示為

(1)

(2)

SST(shear stress transport)k-ω湍流模型，結(jié)合了標(biāo)準(zhǔn)k-ε模型和k-ω模型的各自優(yōu)點(diǎn)，在外部自由流中采用k-ε模型，在近壁區(qū)采用k-ω模型。一般認(rèn)為SSTk-ω湍流模型是雷諾時(shí)均模型中求解精度最好的湍流模型之一，計(jì)算量也顯著小于大渦模擬。由于本文渦激振動(dòng)模擬涉及大量時(shí)間步上的非定常計(jì)算，綜合計(jì)算量和計(jì)算精度，本文數(shù)值模擬采用SSTk-ω湍流模型。

為模擬高階模態(tài)渦激振動(dòng)，本文針對(duì)大帶東橋?qū)嶋H觀測到的第八階豎彎渦激振動(dòng)開展研究。為此將加勁梁斷面豎向渦激振動(dòng)系統(tǒng)簡化為單自由度的彈簧-質(zhì)量-阻尼系統(tǒng)，如圖2所示，其中α表示來流風(fēng)攻角。其豎彎運(yùn)動(dòng)方程表示為

(3)

利用Newmark-β法求解式(3)的運(yùn)動(dòng)微分方程，對(duì)下一時(shí)刻t+Δt，振動(dòng)控制方程可以寫為

FL(t+Δt)/M

(4)

式中,Δt為時(shí)間步長。因此，在t+Δt時(shí)刻的位移、速度和加速度可以表示為

(5)

(6)

其中,

(8)

(9)

式(5)～式(9)中的a0,a1,a2,…,a7均為常數(shù)，取值為

(10)

a7=βΔt

(11)

研究表明，當(dāng)β≥0.5且γ≥0.25(0.5+β)2時(shí)，Newmark-β法將無條件穩(wěn)定，因此系數(shù)β與γ本文分別取為0.5和0.25。

橋梁渦激振動(dòng)為流固耦合問題，本文采用如圖3所示的分析流程。也即先通過Fluent求解扁平箱梁繞流的Navier-Stokes方程，得到作用在扁平箱梁上的升力，再基于Newmark-β求解加勁梁的運(yùn)動(dòng)方程，獲取加勁梁的振動(dòng)位移、速度和加速度，然后通過Fluent宏更新加勁梁姿態(tài)并重新進(jìn)行網(wǎng)格劃分，再進(jìn)入下一個(gè)時(shí)間步上計(jì)算域的求解，依次循環(huán)，如圖3所示。

圖2 加勁梁渦激振動(dòng)模型Fig.2 Vortex-induced vibration (VIV) model of stiffening girder

圖3 渦激振動(dòng)分析流程Fig.3 Flowchat of simulation of VIV

1.2 計(jì)算域及CFD相關(guān)參數(shù)設(shè)定

大帶東橋主橋加勁梁斷面CFD模擬的計(jì)算域，如圖4所示。CFD模擬的箱梁采用1∶80的縮尺比，并考慮了上下游欄桿和中央防撞欄。計(jì)算域的入口、上側(cè)和下側(cè)邊界到加勁梁剪切中心的距離均為10B，出口到斷面剪切中心的距離為25B，對(duì)應(yīng)的模型堵塞度為0.72%，顯著小于風(fēng)洞試驗(yàn)?zāi)Ｐ投氯鹊囊?。為保證模型運(yùn)動(dòng)時(shí)網(wǎng)格的質(zhì)量，本文采用“剛性運(yùn)動(dòng)區(qū)域+變形網(wǎng)格區(qū)域+靜止網(wǎng)格區(qū)域”的方法對(duì)計(jì)算區(qū)域進(jìn)行分區(qū)網(wǎng)格劃分。在扁平箱梁外圍劃分一個(gè)橢圓形區(qū)域，橢圓形區(qū)域的中心與扁平箱梁剪心重合，內(nèi)部采用結(jié)構(gòu)化網(wǎng)格，欄桿及分隔帶附近區(qū)域?yàn)榉墙Y(jié)構(gòu)化網(wǎng)格，整個(gè)橢圓形區(qū)域作為剛體隨扁平箱梁一起運(yùn)動(dòng)。因此，只有與橢圓區(qū)域相連的外部變形區(qū)域網(wǎng)格發(fā)生重構(gòu)，而緊鄰模型表面的邊界層網(wǎng)格采用固定不變的高質(zhì)量網(wǎng)格。通過用戶自定義函數(shù)給定橢圓區(qū)域內(nèi)剛性網(wǎng)格的運(yùn)動(dòng)方式，變形域內(nèi)的非結(jié)構(gòu)化網(wǎng)格通過“彈簧光順”和“區(qū)域重構(gòu)”進(jìn)行更新；另外，靜止網(wǎng)格區(qū)網(wǎng)格固定不動(dòng)，并采用結(jié)構(gòu)化網(wǎng)格劃分。

節(jié)段模型風(fēng)洞試驗(yàn)評(píng)價(jià)橋梁的渦激振動(dòng)通常基于均勻來流條件。為此，本文計(jì)算域入口邊界采用水平均勻流速度，出口邊界設(shè)置為壓力自由出流；扁平箱梁表面為無滑移壁面，上下邊界均設(shè)置為對(duì)稱邊界。采用入口速度對(duì)初始流場進(jìn)行全局初始化。數(shù)值計(jì)算采用循序漸進(jìn)的模擬方法，以平抑初始場計(jì)算導(dǎo)致的數(shù)值振蕩，并最終轉(zhuǎn)到SSTk-ω湍流模型計(jì)算，此時(shí)計(jì)算采用非定常二階隱式格式，采用速度-壓力解耦的SIMPLEC算法；動(dòng)量方程、湍動(dòng)能方程及耗散率方程均采用QUICK格式。數(shù)值模擬時(shí)先進(jìn)行扁平箱梁斷面的靜止繞流計(jì)算，待流場收斂后，再導(dǎo)入用戶自定義函數(shù)進(jìn)行扁平箱梁的高階渦激振動(dòng)模擬。

圖4 計(jì)算域分區(qū)Fig.4 Computational domain partition

1.3 網(wǎng)格及時(shí)間步長無關(guān)性檢查

合適的網(wǎng)格劃分對(duì)CFD精確模擬至關(guān)重要。為此，本文基于G1，G2和G3三套網(wǎng)格開展了網(wǎng)格無關(guān)性檢查。相關(guān)網(wǎng)格參數(shù)見表1所示。三套網(wǎng)格均滿足SST湍流模型對(duì)邊界層網(wǎng)格的要求(Y+<1)，其中最粗網(wǎng)格G1的物面最大Y+為0.96，位于模型前緣附近，如圖5所示。

圖5 G1網(wǎng)格主梁表面Y+分布Fig.5 Distribution of Y+ on girder surface of mesh G1

另外，以G2網(wǎng)格為例，通過GAMBIT之EquiAngle Skew評(píng)價(jià)的網(wǎng)格優(yōu)良質(zhì)量(0～0.25)率為96.25%，而網(wǎng)格質(zhì)量為好的百分率(0～0.5)達(dá)到100%。G2網(wǎng)格方案扁平箱梁周圍網(wǎng)格劃分見圖6；為提高節(jié)段模型前緣所劃分網(wǎng)格的正交性、保持網(wǎng)格邊長的平順變化并顯著減小網(wǎng)格數(shù)量，對(duì)前緣棱角采用半徑B/1550的圓化處理，這種處理對(duì)主梁氣動(dòng)特性的影響可忽略不計(jì)[16]。圖7展示了扁平箱梁前緣、中央防撞欄及上游欄桿周圍網(wǎng)格劃分細(xì)節(jié)，下游欄桿網(wǎng)格布置與圖7(c)相同。

圖6 G2網(wǎng)格在模型周圍的分布Fig.6 Arrangement of G2 mesh around stiffening girder

扁平箱梁斷面的氣動(dòng)升力系數(shù)CL、阻力系數(shù)CD和斯托羅哈數(shù)St分別定義為

CL=FL/(0.5ρU2B);CD=FD/(0.5ρU2B);

St=fsH/U

(12)

式中：U為計(jì)算域入口來流風(fēng)速；FL和FD分別為作用在加勁梁單位長度上的升力和阻力；fs為漩渦脫落頻率；H為扁平箱梁迎風(fēng)高度。

圖7 G2局部位置網(wǎng)格放大Fig.7 Close-up view of G2 mesh at local position

表1 不同網(wǎng)格固定模型繞流模擬結(jié)果Tab.1 Numerical results of flow around stationary stiffening girder on different grid arrangement

表2 自由振動(dòng)模型網(wǎng)格無關(guān)性檢查結(jié)果Tab.2 Numerical results of grid independence check for free oscillating stiffening girder

考慮時(shí)間步長的合適取值，本文采用G2網(wǎng)格，計(jì)算了0°攻角下四個(gè)不同時(shí)間步長0.001 s,0.002 s,0.005 s和0.010 s上的靜止繞流場，雷諾數(shù)為4×104。表3給出了不同時(shí)間步長的模擬結(jié)果，可見當(dāng)時(shí)間步長小于0.002 s時(shí)，氣動(dòng)力系數(shù)平均和脈動(dòng)值趨于穩(wěn)定。另外，對(duì)4°攻角下的扁平箱梁斷面渦激振動(dòng)繞流開展了0.002 s和0.001 s兩個(gè)不同時(shí)間步長下的數(shù)值模擬，模擬結(jié)果如表4所示，兩個(gè)時(shí)間步長上的模擬結(jié)果誤差很小，因此認(rèn)為采用0.002 s的時(shí)間步長是合適的。表3將模擬結(jié)果與Larsen和祝志文等的風(fēng)洞試驗(yàn)結(jié)果進(jìn)行了對(duì)比，可見本文數(shù)值模擬得到的漩渦脫落St數(shù)及阻力系數(shù)均值均與風(fēng)洞試驗(yàn)值吻合。

表3 不同時(shí)間步長靜止繞流數(shù)值結(jié)果Tab.3 Numerical results of stationary stiffening girder under different time step size

表4 自由振動(dòng)時(shí)間無關(guān)性研究數(shù)值結(jié)果Tab.4 Numerical results of time independence check of stiffening girder under free vibration

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

Larsen給出了大帶東橋節(jié)段模型風(fēng)洞試驗(yàn)獲得的實(shí)橋第八階豎向渦激振動(dòng)的完整鎖定曲線，見圖8。Larsen研究中節(jié)段模型試驗(yàn)結(jié)構(gòu)參數(shù)如下：模型寬度B和高度H分別為0.387 5 m和0.055 m，第八階模態(tài)質(zhì)量M=3.553 kg/m、模態(tài)頻率fn=3.9 Hz和模態(tài)阻尼比ζ=0.16%，來流風(fēng)攻角為4°。

2.1 渦激振動(dòng)鎖定曲線

定義折減風(fēng)速為Vr=U/(fnB)，本文對(duì)折減風(fēng)速在0.6～1.8內(nèi)的扁平箱梁高階渦激振動(dòng)進(jìn)行了數(shù)值模擬。從圖8無量綱渦振振幅RMS(root mean square)隨折減風(fēng)速的變化可見，大帶橋加勁梁第八階模態(tài)存在兩個(gè)渦振鎖定區(qū)間，高折減風(fēng)速渦振鎖定區(qū)間與試驗(yàn)有較好的一致性，表現(xiàn)出高階渦振幅值在鎖定區(qū)間內(nèi)先增大后減小的自限幅特征。低折減風(fēng)速渦振鎖定區(qū)間在Larsen研究中沒有體現(xiàn)，這可能與風(fēng)洞試驗(yàn)風(fēng)速縮尺比較大，而低折減風(fēng)速渦振鎖定區(qū)風(fēng)速太小有關(guān)。實(shí)際上對(duì)不少風(fēng)工程風(fēng)洞，試驗(yàn)段風(fēng)速開始穩(wěn)定的風(fēng)速較高，可能高于試驗(yàn)風(fēng)速縮尺對(duì)應(yīng)的實(shí)橋渦振風(fēng)速，因此將無法搜索出低折減風(fēng)速渦振鎖定區(qū)。這顯示了CFD模擬捕捉渦振研究上的優(yōu)勢。

圖8 無量綱渦振振幅RMS隨折減風(fēng)速的變化Fig.8 Non-dimensional RMS value of VIV versus reduced velocity

圖9(a)～圖9(c)分別為圖8中第二個(gè)高階渦振風(fēng)速鎖定區(qū)間的起始點(diǎn)、上升點(diǎn)和最高點(diǎn)對(duì)應(yīng)的加勁梁斷面高階渦激振動(dòng)響應(yīng)時(shí)程曲線，折減風(fēng)速分別為0.93，1.12和1.26，圖中重力加速度g=9.8 m/s2。從圖9中可以看出，高階渦激振動(dòng)時(shí)程曲線均是振幅先增大后趨于穩(wěn)定。圖10(a)～圖10(c)分別給出了三個(gè)折減風(fēng)速下，加勁梁斷面升力時(shí)程的歸一化功頻譜密度(power spectral density,PSD)曲線?？梢娙齻€(gè)卓越頻率皆與結(jié)構(gòu)的自振頻率3.9 Hz接近，即橋梁結(jié)構(gòu)頻率在較寬的風(fēng)速范圍內(nèi)“鎖定”了扁平箱梁的高階漩渦脫落頻率。

圖9 不同折減風(fēng)速下的渦激振動(dòng)響應(yīng)Fig.9 VIV response under different reduced wind velocities

圖10 不同折減風(fēng)速下的歸一化升力系數(shù)功率譜Fig.10 Normalized PSD of lift coefficient time histories under different reduced wind velocities

圖11和圖12分別為鎖定區(qū)間峰值點(diǎn)，即折減風(fēng)速為1.26時(shí)，高階渦激振動(dòng)達(dá)到穩(wěn)定狀態(tài)后一個(gè)周期內(nèi)升力系數(shù)和振動(dòng)位移時(shí)程曲線，以及扁平箱梁斷面漩渦脫落云圖。從圖可見，升力時(shí)程曲線與位移時(shí)程曲線并不完全同步，位移曲線滯后力系數(shù)曲線5°～10°的相位差。圖12(a)～圖12(d)描述了一個(gè)周期內(nèi)扁平箱梁斷面周圍漩渦的形成與發(fā)展過程(n為正整數(shù))，可見漩渦從扁平箱梁前緣和上游側(cè)欄桿形成，隨后沿箱梁頂面向下游漂移，并最終在下游側(cè)欄桿處脫落的過程。另外，高階渦振一個(gè)周期內(nèi)箱梁下表面未見明顯的漩渦形成和運(yùn)動(dòng)，可見上下游側(cè)橋面欄桿對(duì)扁平箱梁斷面漩渦的形成和脫落有重要影響，因此，大跨度橋梁需特別重視上下游側(cè)橋面欄桿的氣動(dòng)外形設(shè)計(jì)。

圖11 一個(gè)振動(dòng)周期內(nèi)升力及位移時(shí)程曲線Fig.11 Lift and displacement time histories during one cycle of VIV

圖12 渦振鎖定區(qū)一個(gè)周期內(nèi)的漩渦脫落云圖Fig.12 Vortex-shedding contours during one cycle of VIV

2.2 渦激振動(dòng)影響因素

表5 阻尼比對(duì)高階渦激振動(dòng)響應(yīng)的影響Tab.5 Effects of damping ratios on VIV response

表6 攻角對(duì)高階渦激振動(dòng)響應(yīng)的影響Tab.6 Effects of wind angles of attack on VIV response

3 結(jié) 論

本文開展了雷諾數(shù)為3.18×104～6.10×104的扁平箱梁高階模態(tài)渦激振動(dòng)響應(yīng)的CFD模擬，得到如下結(jié)論：

(1)采用本文CFD模擬和結(jié)構(gòu)響應(yīng)計(jì)算的耦合方法，能有效預(yù)測扁平箱梁的高階模態(tài)渦振鎖定區(qū)間和渦振幅值；在低折減風(fēng)速渦振鎖定區(qū)的識(shí)別上，CFD方法可能優(yōu)于風(fēng)洞試驗(yàn)。

(2)高階模態(tài)渦激振動(dòng)表現(xiàn)出與阻尼和來流風(fēng)攻角的相關(guān)性。當(dāng)阻尼比增大時(shí)，高階模態(tài)渦振逐漸減小甚至消失；另外，在研究的風(fēng)攻角范圍，而只有當(dāng)來流風(fēng)攻角大于2°時(shí)，才能觀察到扁平箱梁明顯的高階渦激振動(dòng)。

(3)上下游橋面欄桿對(duì)扁平箱梁上表面漩渦的產(chǎn)生和脫落有重要影響，因此，從預(yù)防渦激振動(dòng)發(fā)生的角度，大跨度橋梁需特別重視上下游橋面欄桿的氣動(dòng)外形設(shè)計(jì)。

由于計(jì)算資源的限制，本文僅開展了二維CFD模擬，進(jìn)一步研究可考慮扁平箱梁三維CFD模擬。當(dāng)未來計(jì)算資源足夠時(shí)，甚至可嘗試開展加勁梁真實(shí)的高階模態(tài)渦激振動(dòng)模擬。