祝志文， 石亞光， 顏 爽

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

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

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

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

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

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

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

1.1 控制方程

計算風工程中，橋梁繞流問題的二維黏性不可壓雷諾時均Navier-Stokes方程表示為

(1)

(2)

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

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

(3)

利用Newmark-β法求解式(3)的運動微分方程，對下一時刻t+Δt，振動控制方程可以寫為

FL(t+Δt)/M

(4)

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

(5)

(6)

其中,

(8)

(9)

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

(10)

a7=βΔt

(11)

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

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

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

圖3 渦激振動分析流程Fig.3 Flowchat of simulation of VIV

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

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

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

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

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

合適的網(wǎng)格劃分對CFD精確模擬至關(guān)重要。為此，本文基于G1，G2和G3三套網(wǎng)格開展了網(wǎng)格無關(guān)性檢查。相關(guān)網(wǎng)格參數(shù)見表1所示。三套網(wǎng)格均滿足SST湍流模型對邊界層網(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評價的網(wǎng)格優(yōu)良質(zhì)量(0～0.25)率為96.25%，而網(wǎng)格質(zhì)量為好的百分率(0～0.5)達到100%。G2網(wǎng)格方案扁平箱梁周圍網(wǎng)格劃分見圖6；為提高節(jié)段模型前緣所劃分網(wǎng)格的正交性、保持網(wǎng)格邊長的平順變化并顯著減小網(wǎng)格數(shù)量，對前緣棱角采用半徑B/1550的圓化處理，這種處理對主梁氣動特性的影響可忽略不計[16]。圖7展示了扁平箱梁前緣、中央防撞欄及上游欄桿周圍網(wǎng)格劃分細節(jié)，下游欄桿網(wǎng)格布置與圖7(c)相同。

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

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

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

St=fsH/U

(12)

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

圖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 自由振動模型網(wǎng)格無關(guān)性檢查結(jié)果Tab.2 Numerical results of grid independence check for free oscillating stiffening girder

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

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

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

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

Larsen給出了大帶東橋節(jié)段模型風洞試驗獲得的實橋第八階豎向渦激振動的完整鎖定曲線，見圖8。Larsen研究中節(jié)段模型試驗結(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%，來流風攻角為4°。

2.1 渦激振動鎖定曲線

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

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

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

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

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

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

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

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

2.2 渦激振動影響因素

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

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

3 結(jié) 論

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

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

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

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

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