王向陽，張帥輝

（武漢理工大學(xué) 交通與物流工程學(xué)院，武漢 430063）

現(xiàn)代橋梁由于跨徑的增大以致結(jié)構(gòu)剛度不斷下降，而結(jié)構(gòu)的輕柔致使其對外部風(fēng)環(huán)境十分敏感，目前風(fēng)致橋梁振動(dòng)問題已經(jīng)成為大跨橋梁施工及運(yùn)營階段抗風(fēng)設(shè)計(jì)必須著重考慮的控制因素[1]。橋梁渦激振動(dòng)是風(fēng)與結(jié)構(gòu)相互作用產(chǎn)生的一種振動(dòng)形式，是由流經(jīng)橋梁上下表面的流體所產(chǎn)生的交替脫落旋渦引起的限幅振動(dòng)現(xiàn)象。當(dāng)旋渦脫落頻率接近結(jié)構(gòu)的某階自振頻率時(shí)就會(huì)發(fā)生共振，此時(shí)橋面振幅驟增嚴(yán)重影響行車安全以及降低結(jié)構(gòu)疲勞壽命[2]。近幾十年來，由于橋梁渦振現(xiàn)象頻發(fā)，國內(nèi)外學(xué)者對此展開了大量研究，從而在渦振的機(jī)理、影響因素以及抑振措施等方面取得一定的進(jìn)展[3]。其主要研究方法逐漸從完全依賴傳統(tǒng)的風(fēng)洞試驗(yàn)過渡到數(shù)值模擬與試驗(yàn)相結(jié)合，這一定程度上彌補(bǔ)了風(fēng)洞試驗(yàn)耗資大、周期長等不足，同時(shí)也說明數(shù)值模擬技術(shù)的可行性。國外學(xué)者Frandsen 等[4]基于有限元方法和動(dòng)網(wǎng)格技術(shù)對大帶東橋進(jìn)行渦激振動(dòng)的全耦合模擬，其計(jì)算結(jié)果同其他研究數(shù)據(jù)相比較為吻合；Uenal 等[5]通過對一個(gè)近尾跡流場的計(jì)算分析發(fā)現(xiàn)4 種基于RANS（Reynolds Average Navier-Stokes）的湍流模型中SST（Shear Stress Transport）湍流模型模擬結(jié)果最準(zhǔn)確；Daniels 等[6]采用大渦模型對矩形柱進(jìn)行強(qiáng)迫振動(dòng)模擬，得到了比較準(zhǔn)確的計(jì)算結(jié)果；國內(nèi)學(xué)者徐楓 等[7]將Newmark-β法寫入U(xiǎn)DF（Users Defined Functions）并結(jié)合FLUENT 中的動(dòng)網(wǎng)格技術(shù)成功觀測到柱體的鎖定、“拍”以及位移失諧的現(xiàn)象；陳文禮等[8]基于雷諾時(shí)均應(yīng)力模型模擬得到了圓柱渦致振動(dòng)的風(fēng)速鎖定區(qū)域，且與試驗(yàn)結(jié)果吻合較好；李永樂等[9]基于CFD（Computational Fluid Dynamics）和CSD（Computational Structural Dynamics）耦合的方法較好地模擬了方柱風(fēng)致振動(dòng)中的流固耦合效應(yīng)。

以上研究的斷面形狀較為規(guī)則，且多為風(fēng)洞試驗(yàn)的驗(yàn)證性研究。本文以π 形主梁斷面為研究對象，基于松耦合方法通過Newmark-β數(shù)值解法結(jié)合FLUENT中的二次開發(fā)接口來獲取斷面的渦激振動(dòng)響應(yīng)，采用“剛性運(yùn)動(dòng)區(qū)域+動(dòng)網(wǎng)格區(qū)域+靜止網(wǎng)格區(qū)域”的方法分塊劃分網(wǎng)格[10]，通過網(wǎng)格光順和重構(gòu)改善網(wǎng)格運(yùn)動(dòng)過程中的質(zhì)量問題。在驗(yàn)證該方法正確性的基礎(chǔ)上，建立一座擬建斜拉橋的有限元?jiǎng)恿δＰ?，根?jù)獲取的動(dòng)力特性參數(shù)對其進(jìn)行風(fēng)致渦振模擬，根據(jù)該模擬結(jié)果可初步判斷該橋梁的渦振性能，同時(shí)也可為后續(xù)的風(fēng)洞試驗(yàn)奠定一定基礎(chǔ)。

1 松耦合方法

松耦合法在不同學(xué)科中均有廣泛應(yīng)用，對于氣動(dòng)彈性問題，王少波等[11]通過CFD+CSD結(jié)合的方式給出基于松耦合法求解該問題的一般思路。針對本文中的渦激振動(dòng)問題，首先在分析結(jié)構(gòu)的豎向振動(dòng)時(shí)將其看作是單自由度的彈簧振子系統(tǒng)，然后利用FLUENT 中的二次開發(fā)功能，將求解結(jié)構(gòu)豎向動(dòng)力方程的Newmark-β方法以及FLUENT中特有的宏命令寫入自編UDF 中與FLUENT 進(jìn)行對接。開始計(jì)算時(shí)先求解每個(gè)時(shí)間步內(nèi)系統(tǒng)的流體控制方程，并通過宏Compute_Force_And_Moment 提取結(jié)構(gòu)上的升力，然后通過上述UDF中的Newmark-β數(shù)值方法求解在升力作用下結(jié)構(gòu)的豎向運(yùn)動(dòng)方程從而得到位移、速度等物理量；再將在上一步得到的物理量通過宏DEFINE_CG_MOTIONG賦予給結(jié)構(gòu)及周圍的動(dòng)網(wǎng)格由此更新結(jié)構(gòu)的運(yùn)動(dòng)狀態(tài)，待該時(shí)間步內(nèi)計(jì)算收斂后進(jìn)入下一時(shí)間步，由此不斷循環(huán)迭代直至所監(jiān)控的物理量趨于穩(wěn)定。具體求解思路如圖1所示。

圖1 基于松耦合法求解豎向渦激振動(dòng)示意圖

2 渦激振動(dòng)的模擬驗(yàn)證

為了檢驗(yàn)本文基于松耦合方法求解豎向渦激振動(dòng)數(shù)值模擬方法的正確性，選取文獻(xiàn)[12]中的矩形斷面進(jìn)行渦激響應(yīng)計(jì)算。其中矩形斷面高為0.075 m，寬為0.3 m，為與風(fēng)洞試驗(yàn)數(shù)據(jù)對比只考慮豎向自由度，因此將其簡化為只有豎向自由度的彈簧振子系統(tǒng)[13]。渦激振動(dòng)模擬中的設(shè)計(jì)參數(shù)均取自文獻(xiàn)[12]，其中單位長度質(zhì)量為1.362 kg，豎彎頻率fn為3.905 Hz，豎彎阻尼比為0.177%。

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

為了使進(jìn)口處流動(dòng)能夠發(fā)展完全以及出口流完全流出，從而不對流動(dòng)產(chǎn)生影響，要保證數(shù)值模擬的阻塞率小于3%[14]，因此將計(jì)算域大小取為18B×12B的矩形，速度進(jìn)口邊界和上下對稱邊界距離矩形斷面形心6B，壓力出口邊界距離斷面形心12B，其中B為模型特征長度（B=0.3 m）。計(jì)算域以及邊界條件的設(shè)置如圖2所示，其中矩形斷面為無滑移邊界。

圖2 計(jì)算域及邊界條件設(shè)置

在劃分網(wǎng)格時(shí)考慮到梁斷面在運(yùn)動(dòng)過程中所導(dǎo)致的網(wǎng)格畸變、負(fù)體積等問題，采用分塊的方式將網(wǎng)格分為剛性運(yùn)動(dòng)、動(dòng)網(wǎng)格以及靜止網(wǎng)格3 大區(qū)域。其中剛性運(yùn)動(dòng)區(qū)域距離梁斷面最近，設(shè)置較密的網(wǎng)格，該部分網(wǎng)格跟隨梁斷面一起運(yùn)動(dòng)，且為了能夠捕捉到氣流流經(jīng)主梁斷面時(shí)發(fā)生的“分離”和“再附”現(xiàn)象，對主梁斷面設(shè)置邊界層網(wǎng)格；在動(dòng)網(wǎng)格區(qū)域設(shè)置尺寸較大三角形網(wǎng)格來解決剛性域運(yùn)動(dòng)帶來的網(wǎng)格質(zhì)量問題；靜止網(wǎng)格區(qū)域距離梁斷面最遠(yuǎn)，對其采用漸變率較大的網(wǎng)格以提高計(jì)算效率，整個(gè)計(jì)算域的網(wǎng)格總數(shù)為46 893，其中第一層網(wǎng)格高度設(shè)為0.001B，近壁面網(wǎng)格劃分如圖3所示。

圖3 近壁面網(wǎng)格劃分

2.2 數(shù)值結(jié)果驗(yàn)證分析

采用k-ωSST湍流模型進(jìn)行豎向振動(dòng)瞬態(tài)繞流計(jì)算，通過UDF 控制梁體，經(jīng)過充分繞流后釋放梁體，分別對折算風(fēng)速為4、6、8、8.5、9、9.5、10、12、14、16 共計(jì)10 個(gè)風(fēng)速工況進(jìn)行豎向單自由度的渦激振動(dòng)計(jì)算，由于篇幅所限部分折算風(fēng)速下的位移時(shí)程計(jì)算結(jié)果如圖4所示。

圖4 不同折算風(fēng)速下矩形斷面位移時(shí)程曲線

由圖4可知矩形斷面在不同折算風(fēng)速下的振動(dòng)幅值差別很大，在Vr=6時(shí)振幅較小，在Vr=9.5時(shí)振幅達(dá)到最大，當(dāng)風(fēng)速增加到Vr=12 時(shí)振幅又減小，可見渦激振動(dòng)時(shí)的風(fēng)速在6～12之間。

圖5所示為無量綱振幅隨折減風(fēng)速變化的折線圖，從圖中可以看出本文模擬結(jié)果中除渦振結(jié)束時(shí)的振幅有一定差異外，起振風(fēng)速、最大振幅以及風(fēng)速“鎖定區(qū)間”同文獻(xiàn)試驗(yàn)數(shù)據(jù)吻合度較高；圖6所示為漩渦脫落頻率與結(jié)構(gòu)固有頻率比（fs/fn）隨折減風(fēng)速變化的折線圖，從圖中可以看風(fēng)速在6～12 之間時(shí)的fs/fn值在1.0附近，此時(shí)漩渦脫落頻被結(jié)構(gòu)頻率“捕獲”。以上結(jié)果驗(yàn)證了本文求解豎向渦激振動(dòng)問題數(shù)值模擬方法的可靠性，該方法可用于后續(xù)的分析。

圖5 無量綱振幅隨折減風(fēng)速的變化曲線

圖6 頻率比隨折減風(fēng)速的變化曲線

3 漢江大橋渦激振動(dòng)模擬

3.1 有限元模型建立及動(dòng)力特性計(jì)算

某擬建漢江大橋橋型為雙塔雙索面組合梁斜拉橋，主橋立面采用對稱布置，主梁采用半漂浮體系。全橋的跨徑布置為（50+112）m+370 m+（112+50）m，主橋長度為694 m，主跨跨度為370 m，邊跨設(shè)置1個(gè)輔助墩和1個(gè)過渡墩，邊跨跨度為162 m，邊中跨比為0.437。橋面全幅寬28.3 m，橋面雙向橫坡的坡率為2%，π形主梁標(biāo)準(zhǔn)斷面如圖7所示。

圖7 π形主梁標(biāo)準(zhǔn)橫斷面

建立該斜拉橋有限元?jiǎng)恿δＰ偷墓ぷ鹘柚浖嗀NSYS來完成。主梁采用單脊梁模式模擬，其中主梁、橋塔和橋墩采用BEAM188 梁單元，斜拉索采用LINK8 桿單元，二期與檢修車道等附屬設(shè)施采用MASS21質(zhì)量點(diǎn)單元，全橋有限元模型如圖8所示。

圖8 漢江大橋有限元?jiǎng)恿δＰ?/p>

表1中給出了橋梁豎彎、扭轉(zhuǎn)的1階自振頻率以及對應(yīng)的等效質(zhì)量等參數(shù)，與表1所對應(yīng)的振型圖如圖9所示。

圖9 漢江大橋振型圖

表1 動(dòng)力特性計(jì)算結(jié)果

3.2 單自由度豎向渦激振動(dòng)模擬

對于渦激振動(dòng)，文獻(xiàn)[15]中表明當(dāng)扭彎基頻比大于1.4時(shí)難以發(fā)生彎扭耦合振動(dòng)的情況，這時(shí)兩種渦振形式相互獨(dú)立。根據(jù)表1可知文中漢江大橋的扭彎基頻比為1.52，兩者基頻相差較大，且根據(jù)渦激振動(dòng)的特點(diǎn)可知豎彎渦振會(huì)先于扭轉(zhuǎn)渦振發(fā)生，因此本文針對豎向自由度渦振進(jìn)行模擬。

3.2.1 計(jì)算域及網(wǎng)格劃分

模擬時(shí)對主梁進(jìn)行一定簡化，暫不考慮護(hù)欄等附屬設(shè)施，斷面幾何參數(shù)如表2所示，計(jì)算域的大小及邊界條件的設(shè)置同圖2中設(shè)置相同，網(wǎng)格劃分方式同2.1小節(jié)中一樣將網(wǎng)格分為剛性運(yùn)動(dòng)、動(dòng)網(wǎng)格以及靜止網(wǎng)格3大區(qū)域，具體網(wǎng)格劃分如圖10所示，第一層網(wǎng)格高度為0.000 2B（B為梁寬）。數(shù)值模擬相關(guān)參數(shù)設(shè)置如表2所示，其中豎彎阻尼比取為0.79%。

表2 模擬相關(guān)參數(shù)設(shè)置

圖10 π梁斷面網(wǎng)格劃分

3.2.2 豎向渦激振動(dòng)模擬結(jié)果

先對主梁斷面進(jìn)行靜止繞流模擬，數(shù)值計(jì)算時(shí)采用k-ωSST雙方程模型，通過Simplec算法求解耦合的速度場和壓力場，時(shí)間步長取為0.001 s。通過計(jì)算來流風(fēng)速為2 m/s～5 m/s范圍內(nèi)的St（斯托羅哈數(shù)）得到其均值為0.06，根據(jù)公式U=fnD/St估計(jì)發(fā)生最大渦振振幅時(shí)的來流風(fēng)速在3.6 m/s左右，因此分別對來流風(fēng)速為1.9 m/s、2.5 m/s、3.1 m/s、3.3 m/s、3.5 m/s、3.7 m/s、4.0 m/s、4.3 m/s、4.9 m/s的9個(gè)工況在0°風(fēng)攻角下進(jìn)行渦激響應(yīng)計(jì)算。限于篇幅，本文只給出了渦振起振振幅、較大振幅、最大振幅和振幅開始減小時(shí)所對應(yīng)風(fēng)速下的位移時(shí)程曲線以及FFT（Fast Fourier Transform）變化頻譜曲線，如圖11至圖14所示。

圖11 來流風(fēng)速為2.5 m/s時(shí)的計(jì)算結(jié)果

圖12 來流風(fēng)速為3.3 m/s時(shí)的計(jì)算結(jié)果

圖13 來流風(fēng)速為3.7 m/s時(shí)的計(jì)算結(jié)果

如圖11至圖14所示，在不同風(fēng)速下π梁斷面渦激振動(dòng)幅值變化很大，在來流風(fēng)速為2.5 m/s時(shí)振幅很小，對應(yīng)梁斷面的旋渦脫落頻率（卓越頻率）為2.09 Hz，此時(shí)還未被結(jié)構(gòu)自振頻率捕獲；當(dāng)來流風(fēng)速增大為3.3 m/s 時(shí)，梁斷面振幅迅速增大，渦脫頻率接近結(jié)構(gòu)自振頻率，此時(shí)發(fā)生“鎖定現(xiàn)象”；當(dāng)來流風(fēng)速進(jìn)一步增大到3.7 m/s時(shí)渦激振幅達(dá)到最大值，此時(shí)渦脫頻率與結(jié)構(gòu)豎頻基本一致；當(dāng)風(fēng)速增大至4.0 m/s 時(shí)梁斷面的旋渦脫落頻率繼續(xù)增大，振幅回落到較小值，“鎖定現(xiàn)象”消失。

圖14 來流風(fēng)速為4.0 m/s時(shí)的計(jì)算結(jié)果

圖15至圖16分別為梁斷面在1/45 縮尺下振幅隨風(fēng)速變化的折線圖和旋渦脫落頻率隨來流風(fēng)速變化的折線圖。

圖15 振幅隨風(fēng)速的變化曲線

從圖15中可以看出，π梁斷面在2.5 m/s～4 m/s的風(fēng)速范圍內(nèi)振幅較大，發(fā)生渦激振動(dòng)現(xiàn)象，在來流為3.7 m/s時(shí)振動(dòng)幅值最大為2.57 mm；

從圖16可以看出，風(fēng)速在3.0 m/s～4.0 m/s的區(qū)間內(nèi)時(shí)旋渦脫落頻率變化很小且接近主梁豎向固有頻率，即發(fā)生了“鎖定現(xiàn)象”。

圖16 旋渦脫落頻率隨來流風(fēng)速的變化曲線

圖17為來流風(fēng)速為3.7 m/s 時(shí)π 梁斷面某一周期內(nèi)的瞬時(shí)渦脫變化圖。

圖17 π梁斷面瞬時(shí)渦脫演化圖

從圖17中可以看出，在梁的不同位置處均有不同尺度的漩渦出現(xiàn)且主要產(chǎn)生自梁前緣的上下側(cè)以及梁后緣下側(cè)。對于梁前緣下側(cè)的旋渦，從圖17(a)中可以看出，在初始nT時(shí)刻剛生成的旋渦準(zhǔn)備向下游移動(dòng)；在nT+1/3時(shí)刻移動(dòng)至梁1/2處；在nT+2/3時(shí)刻移至梁后緣并且此時(shí)梁前緣下側(cè)的旋渦已形成較大規(guī)模；在nT+T時(shí)刻已經(jīng)脫離梁斷面向右下方移動(dòng)，此時(shí)梁前緣下側(cè)的旋渦經(jīng)過發(fā)展已具有足夠規(guī)模又開始準(zhǔn)備移動(dòng)，至此完成一個(gè)脫落周期。對于梁前緣上側(cè)以及后緣下側(cè)的旋渦而言，前者在經(jīng)歷分離、再附后與后者形成交替出現(xiàn)在梁尾流區(qū)域的卡門渦脫。由以上分析可知，氣流在梁迎風(fēng)面發(fā)生分離后分別在梁上下側(cè)的某些固定位置形成旋渦，旋渦不斷后移至梁背，此時(shí)梁背上側(cè)與下側(cè)的旋渦相互合并后在梁背風(fēng)面形成上下不斷脫落的旋渦，從而導(dǎo)致主梁發(fā)生豎向的渦激振動(dòng)。

4 結(jié)語

本文借助流體分析軟件FLUENT，基于松耦合方法通過嵌入Newmark-β法的UDF程序，求解某擬建漢江斜拉橋π 型主梁斷面的豎向渦激振動(dòng)問題，得到如下結(jié)論：

（1）通過對寬高比為4 的矩形斷面的渦激振動(dòng)模擬，證明了采用該模擬方法求解渦激振動(dòng)問題的可行性和正確性。

（2）通過對9 個(gè)風(fēng)速工況下的π 型主梁進(jìn)行渦激振動(dòng)模擬確定了“鎖定區(qū)間”，發(fā)現(xiàn)來流在2.5 m/s～4.0 m/s風(fēng)速范圍內(nèi)時(shí)主梁振動(dòng)幅值較大，發(fā)生“鎖定現(xiàn)象”。

（3）根據(jù)對π 型主梁斷面某一完整周期內(nèi)旋渦脫落過程的分析，發(fā)現(xiàn)在梁背風(fēng)面不斷產(chǎn)生上下脫落的旋渦導(dǎo)致了π型主梁斷面的豎向渦激振動(dòng)。