胡朋，陳屹林，蔡春聲，2，林偉，3，韓艷?

（1.長沙理工大學(xué)土木工程學(xué)院，湖南長沙 410114；2.東南大學(xué)交通學(xué)院，江蘇南京 211189；

3.廣州市城建規(guī)劃設(shè)計(jì)院有限公司，廣東廣州 510230）

大跨度鋼橋梁因跨度大、阻尼小、質(zhì)量輕、自振頻率低等特點(diǎn)常易引發(fā)低風(fēng)速下的渦激振動(dòng)（以下簡稱“渦振”）問題，這對(duì)大橋的施工安全及成橋狀態(tài)下的行車安全帶來了嚴(yán)重的安全隱患［1］.丹麥的大海帶橋東主橋?yàn)橹骺? 624 m 的帶風(fēng)嘴鋼箱梁懸索橋，在主梁架設(shè)及橋面鋪裝期間都發(fā)生過較大幅值的豎向渦振［2］，而且在之后的運(yùn)營期間又觀測到數(shù)十次的豎向渦振現(xiàn)象，其中最大幅值超過30 cm.日本東京灣聯(lián)絡(luò)橋?yàn)樽畲笾骺?40 m 的矩形鋼箱連續(xù)梁橋［3］，在約16 m∕s 的風(fēng)速下發(fā)生了豎向渦振，最大振幅更是超過了50 cm.在我國，近段時(shí)間也有幾座大橋接連發(fā)生了渦振現(xiàn)象，如2020 年4 月，主跨為850 m 采用雙邊主梁形式的武漢鸚鵡洲長江大橋發(fā)生了輕微渦振；2020 年5 月，主跨為888 m 采用帶風(fēng)嘴鋼箱主梁形式的廣東虎門大橋懸索橋也發(fā)生了較為明顯的渦振現(xiàn)象，并導(dǎo)致大橋交通關(guān)閉.因此，對(duì)于大橋渦振問題應(yīng)給予高度重視.

國內(nèi)外學(xué)者對(duì)于大橋主梁渦振問題開展了大量的研究工作，但對(duì)帶挑臂鋼箱梁橋渦振性能的研究還相對(duì)較少［4］.注意到，相對(duì)于流線型扁平鋼箱梁斷面，帶挑臂鋼箱梁斷面的腹板角度一般較大，在風(fēng)的作用下，容易在迎風(fēng)側(cè)腹板后緣產(chǎn)生較大的旋渦，同時(shí)在背風(fēng)側(cè)腹板與挑臂之間也容易形成較大的旋渦，極易引起大橋主梁的渦振問題.已有試驗(yàn)也表明［4-6］，帶挑臂箱梁斷面容易出現(xiàn)較明顯的渦振振幅.因此，圍繞帶挑臂鋼箱主梁斷面，對(duì)其渦振性能進(jìn)行深入研究非常有必要.

實(shí)際上，在大跨度橋梁抗風(fēng)設(shè)計(jì)時(shí)，很難完全消除橋梁的渦振現(xiàn)象，因此滿足結(jié)構(gòu)安全和使用舒適的小幅值渦振是可以接受的，但這要求我們能詳細(xì)了解橋梁渦振的非線性特性及細(xì)觀機(jī)理，并準(zhǔn)確預(yù)測渦振發(fā)生時(shí)的響應(yīng)值，為此需要建立準(zhǔn)確、可靠的非線性渦激力數(shù)學(xué)模型［7］.目前Scanlan 經(jīng)驗(yàn)非線性模型是橋梁渦振中應(yīng)用最為廣泛的，但孟曉亮［7］、Zhu 等［8］通過大比例節(jié)段模型風(fēng)洞試驗(yàn)發(fā)現(xiàn)，對(duì)于中央開槽箱梁和帶風(fēng)嘴的扁平鋼箱梁而言，采用Scanlan 非線性模型計(jì)算的渦振幅值和渦激力與試驗(yàn)值均存在較為顯著的差別，這也說明了Scanlan 非線性模型存在一定的局限性.針對(duì)帶風(fēng)嘴的扁平鋼箱梁斷面，Zhu 等［9］基于試驗(yàn)測得的渦激力時(shí)程及振動(dòng)響應(yīng)，深入探究了箱梁渦振的機(jī)理及非線性特性，并最終確立了此種斷面的渦激力數(shù)學(xué)模型.進(jìn)一步的，Zhu 等［10］通過分析各渦激力分量對(duì)帶風(fēng)嘴扁平鋼箱梁、中央開槽鋼箱梁以及半封閉鋼箱梁渦振響應(yīng)貢獻(xiàn)的大小，對(duì)箱梁的渦激力模型進(jìn)行了相應(yīng)的簡化，最終提出了一種統(tǒng)一的、形式較為簡單的渦激力模型，經(jīng)過與試驗(yàn)數(shù)據(jù)對(duì)比，表明該模型能較好地預(yù)測箱梁的渦振振幅.Chen 等［11］提出了一種自適應(yīng)非線性擬合方法，考察了各渦激力分量對(duì)帶風(fēng)嘴扁平鋼箱梁渦振貢獻(xiàn)的大小，最終建立了能夠較好預(yù)測箱梁節(jié)段模型渦振振幅的渦激力模型.由此可見，經(jīng)過學(xué)者們的努力，已經(jīng)建立了幾種扁平箱梁斷面形式的渦激力數(shù)學(xué)模型.隨著帶挑臂鋼箱梁斷面形式越來越多地應(yīng)用于大跨度橋梁中，為了能準(zhǔn)確地預(yù)測該類主梁斷面的渦振振幅和進(jìn)一步了解其渦激力的非線性特性，有必要建立可以準(zhǔn)確描述其渦激力的數(shù)學(xué)模型，但目前針對(duì)帶挑臂鋼箱梁斷面渦激力數(shù)學(xué)模型的研究尚未見報(bào)道.

針對(duì)上述問題，本文以某典型帶挑臂鋼箱主梁的斜拉橋?yàn)楣こ瘫尘?，首先，基于雙向流固耦合的數(shù)值模擬方法對(duì)該橋主梁的渦振全過程進(jìn)行了數(shù)值模擬計(jì)算，并與相應(yīng)的風(fēng)洞試驗(yàn)結(jié)果進(jìn)行對(duì)比驗(yàn)證.針對(duì)渦振振幅峰值對(duì)應(yīng)的來流風(fēng)速，計(jì)算了該帶挑臂鋼箱主梁的渦激力時(shí)程，并建立了該主梁形式的渦激力數(shù)學(xué)模型.進(jìn)一步的，基于該渦激力數(shù)學(xué)模型，推導(dǎo)了主梁穩(wěn)態(tài)渦振振幅的表達(dá)式.

1 帶挑臂鋼箱主梁渦振的數(shù)值模擬

1.1 工程簡介

某斜拉橋跨徑布置為85 m+90 m+160 m=335 m，其中主跨為160 m，大橋主梁斷面采用帶挑臂的鋼箱梁形式，梁寬37 m，中心線處梁高3.282 m，寬高比為11.3∶1，如圖1所示.由圖1可知，該橋主梁腹板角度較大，在風(fēng)的作用下，容易引起氣流分離與旋渦脫落，從而引起大橋的渦振問題.鑒于渦振對(duì)成橋態(tài)與施工態(tài)的危害性，有必要對(duì)該種斷面形式的主梁進(jìn)行深入的渦振性能分析.

圖1 主梁截面示意（單位：cm）Fig.1 Schematic diagram of the main girder（unit：cm）

1.2 試驗(yàn)概況

本橋的渦振試驗(yàn)是在長沙理工大學(xué)風(fēng)工程與風(fēng)環(huán)境研究中心的高速試驗(yàn)段中進(jìn)行（如圖2 所示），采用Cobra Probe 三維脈動(dòng)風(fēng)速儀測試來流風(fēng)速，采用激光位移計(jì)測試主梁節(jié)段模型的位移時(shí)程.彈性懸掛的節(jié)段模型縮尺比為1∶50，風(fēng)速比為1∶3.318，其他主要參數(shù)如表1所示.

表1 帶挑臂型鋼箱主梁節(jié)段模型參數(shù)取值Tab.1 Parameter values for the section model of the steel box girder with projecting slab

圖2 放置于風(fēng)洞中的試驗(yàn)?zāi)Ｐ虵ig.2 Testing model in the wind tunnel

數(shù)值模擬時(shí)，為了與風(fēng)洞試驗(yàn)的節(jié)段模型保持一致，數(shù)值模型中包含了主梁、欄桿、檢修車軌道等構(gòu)件，并將其簡化為二維模型，如圖3所示.采用Fluent前處理軟件Gambit對(duì)主梁斷面進(jìn)行網(wǎng)格劃分.計(jì)算區(qū)域尺寸太小會(huì)導(dǎo)致流動(dòng)難以充分發(fā)展，影響邊界條件的設(shè)定，而太大的計(jì)算區(qū)域會(huì)浪費(fèi)計(jì)算資源，增加計(jì)算時(shí)間.綜合考慮計(jì)算精度與計(jì)算效率，并使計(jì)算區(qū)域滿足阻塞率小于5%的要求，本文選擇28D×17B的矩形計(jì)算域，其中D為主梁總高度，B為主梁總寬度.計(jì)算區(qū)域分為剛性運(yùn)動(dòng)區(qū)域、動(dòng)網(wǎng)格區(qū)域和靜止網(wǎng)格區(qū)域［12］.參考已有研究經(jīng)驗(yàn)［11］，剛性運(yùn)動(dòng)區(qū)域尺寸為1.25D× 1.15B，主梁至剛體運(yùn)動(dòng)區(qū)域左邊緣、右邊緣、上邊緣和下邊緣的距離分別為0.05B、0.1B、0.125D和0.125D，動(dòng)網(wǎng)格區(qū)域尺寸為10D×3B；為使尾流流場充分發(fā)展，尾流域尺寸為10D×9.5B，具體區(qū)域劃分如圖4 所示.對(duì)于邊界條件設(shè)置，計(jì)算域入口設(shè)置為速度入口邊界條件，因?yàn)楣?jié)段模型風(fēng)洞試驗(yàn)結(jié)果表明該橋在+3°風(fēng)攻角下發(fā)生了較明顯的豎彎渦振現(xiàn)象.針對(duì)該風(fēng)攻角，將下邊界和左邊界均設(shè)置為速度入口邊界條件，相應(yīng)的右邊界和上邊界均設(shè)置為壓力出口邊界條件，而主梁斷面及附屬結(jié)構(gòu)模型表面則設(shè)置為無滑移的壁面邊界條件，具體如圖4 所示.研究中，采用對(duì)鈍體繞流比較適合的非定常SSTk-ω湍流模型進(jìn)行計(jì)算［13］，由于試驗(yàn)中來流湍流強(qiáng)度在0.5%左右，因此將入口處的湍流強(qiáng)度值設(shè)置為0.5%，同時(shí)根據(jù)軟件常用設(shè)置，將入口處的湍流黏性比設(shè)置為5.0［13］.

圖3 二維主梁斷面Fig.3 Two-dimensional main girder section

圖4 計(jì)算域分塊及邊界條件Fig.4 Computational domain divisions and boundary conditions

在保證數(shù)值模擬計(jì)算精度的前提下，應(yīng)盡可能縮減網(wǎng)格數(shù)量，從而縮短計(jì)算時(shí)間.對(duì)于剛性運(yùn)動(dòng)區(qū)域，為縮減網(wǎng)格數(shù)量并保證網(wǎng)格初始質(zhì)量，本研究中在該區(qū)域采用四邊形網(wǎng)格進(jìn)行劃分.因剛性運(yùn)動(dòng)區(qū)域引發(fā)的位移會(huì)傳遞至動(dòng)網(wǎng)格區(qū)域，故動(dòng)網(wǎng)格區(qū)域的網(wǎng)格需不斷地重構(gòu)更新.為防止網(wǎng)格在重構(gòu)更新時(shí)發(fā)生畸變以致計(jì)算中斷，本研究在動(dòng)網(wǎng)格區(qū)域采用適應(yīng)性更好的三角形網(wǎng)格進(jìn)行劃分，且同時(shí)采用彈簧光順和網(wǎng)格重構(gòu)兩種方法對(duì)動(dòng)網(wǎng)格部分進(jìn)行重構(gòu)更新.對(duì)于靜止網(wǎng)格區(qū)域，為有效地減少網(wǎng)格數(shù)量，對(duì)該區(qū)域全部采用四邊形網(wǎng)格進(jìn)行劃分.對(duì)主梁斷面及附屬設(shè)施設(shè)置邊界層網(wǎng)格，以更好地模擬主梁斷面周圍的流場情況，從而得到更準(zhǔn)確的數(shù)值模擬結(jié)果.同時(shí)，為驗(yàn)證計(jì)算區(qū)域的網(wǎng)格無關(guān)性，先后共劃分生成3套網(wǎng)格，網(wǎng)格總數(shù)分別為18萬、26萬和36 萬.通過試算發(fā)現(xiàn)，26 萬和36 萬網(wǎng)格的計(jì)算結(jié)果非常接近，而18 萬網(wǎng)格的計(jì)算結(jié)果與之相差較大，如表2所示，這說明當(dāng)網(wǎng)格數(shù)達(dá)到26萬時(shí)，再繼續(xù)增加已經(jīng)不會(huì)顯著影響計(jì)算結(jié)果.兼顧計(jì)算效率和計(jì)算精度，最終選用26 萬的網(wǎng)格作為最終計(jì)算網(wǎng)格，此時(shí)第一層網(wǎng)格高度為0.016 mm，在主梁及附屬設(shè)施周圍設(shè)置15 層邊界層網(wǎng)格，以提高計(jì)算精度，具體網(wǎng)格劃分情況如圖5 所示.圖6 為經(jīng)過2 000 步穩(wěn)態(tài)計(jì)算后主梁壁面處無量綱Y+值，由圖可知，本計(jì)算模型的Y+值總體小于1.0，滿足要求.

表2 不同網(wǎng)格計(jì)算結(jié)果對(duì)比Tab.2 Computed results with different meshing schemes

圖5 全局及細(xì)部網(wǎng)格劃分圖Fig.5 Global and local mesh schemes

圖6 沿主梁表面Y+值分布圖Fig.6 Distribution of Y+value along the surface of main girder

1.3 數(shù)值模擬結(jié)果驗(yàn)證

基于計(jì)算流體動(dòng)力學(xué)軟件Fluent，通過編寫UDF計(jì)算程序?qū)ι鲜龆S主梁斷面模型進(jìn)行雙向流固耦合計(jì)算，并采用四階Runge-Kutta 法求解主梁的運(yùn)動(dòng)狀態(tài)，以進(jìn)一步提高計(jì)算精度［14］.同時(shí)利用大型工作站并行計(jì)算以減少運(yùn)行時(shí)間.求解時(shí)采用速度-壓力耦合算法，以二階迎風(fēng)格式對(duì)控制方程中的各參數(shù)進(jìn)行離散，計(jì)算時(shí)先采用SIMPLE 算法進(jìn)行穩(wěn)態(tài)求解，待穩(wěn)態(tài)計(jì)算結(jié)果穩(wěn)定后，再采用SIMPLEC算法進(jìn)行瞬態(tài)計(jì)算求解.瞬態(tài)求解時(shí)，時(shí)間步長經(jīng)過多次試算，最終設(shè)置為0.000 5 s.

研究中采用與節(jié)段模型風(fēng)洞試驗(yàn)相同的試驗(yàn)參數(shù)，具體如表1 所示.考慮風(fēng)洞試驗(yàn)中主梁在+3°風(fēng)攻角下的豎彎渦振最為明顯，因此對(duì)該風(fēng)攻角下二維主梁斷面的渦振過程進(jìn)行了雙向流固耦合數(shù)值模擬，并將數(shù)值模擬結(jié)果與節(jié)段模型試驗(yàn)結(jié)果進(jìn)行對(duì)比，如圖7所示.由圖7可知，數(shù)值模擬計(jì)算得到的帶挑臂鋼箱主梁斷面在不同風(fēng)速下的渦振幅值變化趨勢與通過節(jié)段模型試驗(yàn)得到的基本一致，且二者測得的渦振鎖定區(qū)間基本相同，大致位于5.5～7.5 m∕s之間，其中最大渦振振幅發(fā)生的風(fēng)速均在7.2 m∕s 左右.但同時(shí)也注意到，兩者在相同風(fēng)速下的渦振幅值存在一定的差別，這可能是由風(fēng)洞實(shí)驗(yàn)中模型制作不夠精準(zhǔn)，節(jié)段模型有一定的展向長度或者數(shù)值模擬本身精度不足等原因造成的.但總體而言，數(shù)值模擬計(jì)算得到的渦振振幅與渦振區(qū)間結(jié)果與試驗(yàn)測試結(jié)果基本一致，其中渦振最大振幅相差約6.2%，渦振區(qū)間長度相差約3.6%，可知數(shù)值模擬結(jié)果精度滿足要求.

圖7 主梁渦振振幅隨來流風(fēng)速的變化Fig.7 Amplitudes of the vortex-induced vibration of the main girder changing with the oncoming wind speed

為進(jìn)一步驗(yàn)證數(shù)值模擬結(jié)果，針對(duì)7.2 m∕s 來流風(fēng)速下所發(fā)生的最大渦振振幅，對(duì)比了該風(fēng)速下數(shù)值模擬的主梁位移時(shí)程與試驗(yàn)時(shí)的主梁位移時(shí)程，如圖8（a）所示，從圖8（a）中可以看到，二者幅值和振動(dòng)頻率均存在細(xì)微的差別，試驗(yàn)測得的最大振幅值略小于數(shù)值模擬值，二者對(duì)應(yīng)的均方差分別為1.770 mm 和1.905 mm，相差約7.63%.同時(shí)，對(duì)數(shù)值模擬時(shí)程結(jié)果進(jìn)行頻譜分析得到的渦振位移頻率為9.498 Hz，而二維節(jié)段模型試驗(yàn)的渦振頻率為9.745 Hz，如圖8（b）所示，兩者相差約2.53%，這可能與試驗(yàn)中振動(dòng)系統(tǒng)的質(zhì)量、剛度等參數(shù)難以準(zhǔn)確測量有關(guān).但總體上，數(shù)值模擬的誤差較小，能夠較為準(zhǔn)確地模擬帶挑臂鋼箱主梁的渦振過程.

圖8 主梁渦振最大幅值處的響應(yīng)對(duì)比圖Fig.8 Comparisons of responses of the girder at the maximum amplitude during the vortex-induced vibration

2 帶挑臂鋼箱主梁渦激力模型

2.1 主梁渦激力的提取

數(shù)值模擬計(jì)算過程中，帶挑臂鋼箱主梁渦振發(fā)生時(shí)作用于每延米主梁上的氣動(dòng)力如圖9 所示，圖9中FH、FV和M表示體軸坐標(biāo)系中的氣動(dòng)力.以主梁斷面豎向受力為例，作用于其上的升力FV包含了兩部分，即脈動(dòng)部分fV和平均部分fV0，而脈動(dòng)部分fV則由渦激力fVI和非風(fēng)致附加力fn兩部分組成，有以下表達(dá)式成立：

圖9 帶挑臂鋼箱主梁斷面渦振時(shí)所受氣動(dòng)力示意Fig.9 The aerodynamic forces acting on the steel box girder with projecting slab during the vortex-induced vibration

式中：非風(fēng)致附加力fn是指振動(dòng)的模型與周圍空氣相互作用所產(chǎn)生的氣動(dòng)力.如當(dāng)節(jié)段模型向下振動(dòng)時(shí)，會(huì)擠壓模型下方的空氣一起運(yùn)動(dòng).與此同時(shí)，周圍的空氣又會(huì)及時(shí)補(bǔ)充到模型上方.如此一來，振動(dòng)的模型就會(huì)對(duì)周圍空氣產(chǎn)生受迫振動(dòng)，而周圍空氣反過來會(huì)對(duì)模型產(chǎn)生作用力，為了與來流所產(chǎn)生的自激力相區(qū)別，這種力就稱為非風(fēng)致附加力，它一般可表示為非風(fēng)致附加質(zhì)量和非風(fēng)致附加阻尼的表達(dá)式.

式中：mn和cn分別表示非風(fēng)致附加質(zhì)量和阻尼，可通過零風(fēng)速下的主梁速度、加速度以及所受升力求得［7］.

將式（2）代入式（1）可得主梁渦激力表達(dá)式為：

式（3）即為在CFD 數(shù)值模擬中，求解主梁渦激力的表達(dá)式.

前文已述，在+3°風(fēng)攻角下，當(dāng)來流風(fēng)速為7.2 m∕s時(shí)，主梁的渦振振幅達(dá)到最大值.同時(shí)，也證明了本文數(shù)值模擬的準(zhǔn)確性.因此，針對(duì)該工況，以下就討論渦激力的求解過程.由式（3）可知，在求解渦激力之前，首先應(yīng)求解非風(fēng)致附加質(zhì)量mn與非風(fēng)致附加阻尼cn參數(shù).考慮在零風(fēng)速來流下，主梁的渦激力為零，即此時(shí)主梁所受的脈動(dòng)升力與非風(fēng)致附加力之和應(yīng)為零，如式（3）所示.為此，研究中首先開展了零風(fēng)速下帶挑臂鋼箱主梁的自由衰減振動(dòng)的數(shù)值模擬，計(jì)算中將結(jié)構(gòu)阻尼比設(shè)置為0.3%，初始位移設(shè)置為0.002 5 m，初始速度和初始加速度則都為0.計(jì)算得到的零風(fēng)速來流下該帶挑臂鋼箱梁的自由衰減振動(dòng)的豎向位移時(shí)程如圖10 所示，而計(jì)算得到的主梁脈動(dòng)升力時(shí)程如圖11 所示.同時(shí)，通過圖10 的位移時(shí)程，采用中心差分法可不難得到其速度和加速度時(shí)程.根據(jù)式（3），利用非線性最小二乘法擬合得到的主梁脈動(dòng)升力時(shí)程也在圖11 中給出，由圖11 可知，兩者吻合較好，此時(shí)擬合得到的非風(fēng)致附加質(zhì)量mn和非風(fēng)致附加阻尼cn的值分別為0.522 7 kg∕m 和6.656 5 N·s·m-2.

圖10 主梁斷面零風(fēng)速自由衰減振動(dòng)位移時(shí)程Fig.10 Time-history of displacement of the main girder during free attenuation vibration at zero wind speed

圖11 主梁斷面零風(fēng)速自由衰減振動(dòng)脈動(dòng)升力時(shí)程Fig.11 Time-history of fluctuating lift force of the main girder during free attenuation vibration at zero wind speed

根據(jù)+3°風(fēng)攻角、7.2 m∕s 來流風(fēng)速下計(jì)算的帶挑臂鋼箱主梁的穩(wěn)態(tài)渦振時(shí)程響應(yīng)，可得到其速度時(shí)程如圖12 所示，而加速度時(shí)程則可通過中心差分法計(jì)算得到.圖13 則表示了該帶挑臂鋼箱梁斷面發(fā)生渦振時(shí)的脈動(dòng)升力時(shí)程.根據(jù)式（3），并利用上述求得的非風(fēng)致附加質(zhì)量mn和非風(fēng)致附加阻尼cn參數(shù)，可最終求得該帶挑臂鋼箱主梁斷面的渦激力時(shí)程如圖14 所示.為了分析其渦激力的頻率組成，其頻譜分析如圖15 所示，由圖15 可知，可以看到它存在多個(gè)頻率，除一階主頻外，二階倍頻譜峰值為一階主頻的5.1%，這說明帶挑臂鋼箱主梁的渦激力表現(xiàn)出了一定的非線性特征.

圖12 主梁斷面渦振速度時(shí)程Fig.12 Time-history of velocity of the main girder during the vortex-induced vibration

圖13 主梁斷面渦振脈動(dòng)升力時(shí)程Fig.13 Time-history of fluctuating lift force of the main girder during the vortex-induced vibration

圖14 主梁斷面渦振時(shí)渦激力時(shí)程Fig.14 Time-history of vortex-induced force of the main girder during the vortex-induced vibration

圖15 主梁渦激力時(shí)程頻譜分析圖Fig.15 Spectrum analysis on the time-history of the vortex-induced force of the main girder

2.2 渦激力數(shù)學(xué)模型的建立及驗(yàn)證

2.2.1 渦激力數(shù)學(xué)模型的建立

由于本文不考慮紊流的影響，因此鋼箱梁的渦激力可由風(fēng)速、時(shí)間、鋼箱梁的位移、速度及加速度的函數(shù)加以描述.另一方面，由前文數(shù)值模擬計(jì)算結(jié)果可知，本文研究的帶挑臂鋼箱梁主要以一階主頻發(fā)生豎向渦激振動(dòng)，此時(shí)鋼箱梁的加速度變量可通過位移變量表示.由此，帶挑臂鋼箱梁豎向渦激力可表示為：

從Scanlan 的經(jīng)驗(yàn)線性、非線性渦激力模型以及后續(xù)的針對(duì)扁平鋼箱梁渦激力模型可以看出，渦激力模型通常含有位移y項(xiàng)、速度y˙項(xiàng)、位移y和速度y˙的交叉項(xiàng)以及渦脫力項(xiàng)等.故渦激力在廣義上可表示如下：

式中：ρ表示空氣密度；U表示來流風(fēng)速；D表示主梁梁高；y表示主梁渦振位移；˙表示主梁渦振速度；Pij（i，j不同時(shí)為0）、Vs、ωs和φ表示待定的參數(shù).

由式（5）可知，模型中的一階氣動(dòng)力除渦脫力項(xiàng)Vssin（ωst+φ）外，還包含線性氣動(dòng)阻尼項(xiàng)P10˙∕U和線性氣動(dòng)剛度項(xiàng)P01y∕D，其他的項(xiàng)次均為非線性高階項(xiàng)，這些非線性高階項(xiàng)代表了主梁渦激力的高階倍頻成分.由圖15 可知，渦激力時(shí)程基本上只包含前三階成分，因此研究中渦激力模型式（5）中i與j之和不需超過3.同時(shí)，為了保證渦激力模型的準(zhǔn)確性，渦激力模型中應(yīng)含若干二階和三階項(xiàng)，因此需要合理選取渦激力中的高階項(xiàng)次.

為合理取舍帶挑臂鋼箱梁的渦激力模型項(xiàng)次，本文以數(shù)學(xué)模型重構(gòu)的渦激力與實(shí)際渦激力之間的殘差大小為判別標(biāo)準(zhǔn)，殘差表達(dá)式如式（6）所示.研究中首先使重構(gòu)的渦激力模型中僅包含渦脫力項(xiàng)，然后逐步增加項(xiàng)次并計(jì)算其與實(shí)際渦激力的殘差大小，殘差越小表明通過數(shù)學(xué)模型重構(gòu)的渦激力與實(shí)際渦激力越接近.

式中：R表示殘差值；N表示時(shí)間步數(shù)量；為數(shù)學(xué)模型重構(gòu)的渦激力；fVI表示實(shí)際渦激力.

基于式（5）的渦激力表達(dá)式，重構(gòu)的過程及相應(yīng)的殘差大小如表3 所示，由表3 中工況2 和工況3 可知，一階項(xiàng)中，P10相比P01而言，前者對(duì)渦激力的影響較大；由工況4～工況9 可知，二階項(xiàng)中，P11對(duì)渦激力的影響最大；由工況10～工況13 可知，三階項(xiàng)中，P12和P30對(duì)渦激力的影響均較大.對(duì)比工況14 與工況15 可知，當(dāng)P21與P12搭配時(shí)，總體殘差達(dá)到最小.由以上分析，可建立出帶挑臂鋼箱主梁渦激力數(shù)學(xué)模型表達(dá)式如式（7）所示，即對(duì)應(yīng)于最小殘差的工況14，此時(shí)其殘差僅為0.018 1，這說明重構(gòu)的渦激力與實(shí)際渦激力非常接近.對(duì)于該工況，相應(yīng)的各系數(shù)值如表4所示.

表4 主梁渦激力數(shù)學(xué)模型式（7）中的參數(shù)值Tab.4 Parameter values in the mathematical model of the vortex-induced force for the main girder by Eq.（7）

文獻(xiàn)［10］針對(duì)帶風(fēng)嘴扁平鋼箱梁、中央開槽鋼箱梁以及半封閉鋼箱梁提出了一個(gè)統(tǒng)一的簡化渦激力模型，如下式所示：

式（8）渦激力模型對(duì)應(yīng)的工況如表3 中的工況16，由表3 可知，此時(shí)殘差為0.188 0；而式（7）渦激力模型對(duì)應(yīng)的是工況14，此時(shí)殘差僅為0.018 1.為更直觀地反映式（7）與式（8）模型的差異，兩種模型的渦激力時(shí)程與目標(biāo)值的對(duì)比如圖16（a）所示，其中目標(biāo)值即為CFD 計(jì)算的原始渦激力時(shí)程.由圖16（a）可知，雖然3 條時(shí)程曲線總體上吻合，但在-0.5～1.5 N∕m 范圍內(nèi)，相同時(shí)間步時(shí)，式（8）的計(jì)算值與目標(biāo)值有一定的差距.此外，在-2.5～-2.0 N∕m 范圍內(nèi)，式（8）計(jì)算的渦激力絕對(duì)值也小于目標(biāo)值（如圖16（b）所示）.而對(duì)于式（7）計(jì)算的渦激力，由圖16（a）與圖16（b）可知，相同時(shí)間步時(shí)，式（7）的計(jì)算值與目標(biāo)值幾乎一致，幅值也非常吻合.因此，式（7）的渦激力模型能比較準(zhǔn)確地模擬帶挑臂鋼箱主梁的渦激力值.

表3 不同重構(gòu)渦激力模型殘差大小Tab.3 Residual values of different reconstituted vortex-induced force models

需要說明的是，Simiu 與Scanlan 于1986 年提出了一種線性氣動(dòng)力模型，并被廣泛應(yīng)用，如式（9）所示.該氣動(dòng)力模型對(duì)應(yīng)的工況如表3 中的工況4 所示，可知該工況的殘差為0.097 2，遠(yuǎn)大于式（7）所對(duì)應(yīng)的工況14 的殘差值.為考察式（9）所代表的經(jīng)典線性氣動(dòng)力模型的適用性，圖17 給出了式（9）線性渦激力時(shí)程與式（7）非線性渦激力時(shí)程的局部對(duì)比.由圖17 可知，線性渦激力模型的絕對(duì)值相比非線性渦激力模型的絕對(duì)值要明顯偏小，其幅值偏小約8.8%，而由圖16 可知，式（7）所代表的非線性渦激力與目標(biāo)值基本一致.由此可知，采用線性渦激力模型預(yù)測的渦激力時(shí)程有一定的誤差.由于式（9）的渦激力模型相比式（7）的渦激力模型缺少了高階項(xiàng)，這說明了盡管渦激力的二階及高階成分占比很?。ㄈ鐖D15所示），但在一定程度上影響了渦激力時(shí)程的建模精度.

圖16 不同渦激力模型對(duì)比Fig.16 Comparisons between different vortex-induced force models

圖17 線性渦激力模型與非線性渦激力模型對(duì)比Fig.17 Comparisons between the linear and nonlinear vortex-induced force models

2.2.2 非線性與經(jīng)典線性渦激力數(shù)學(xué)模型位移反算值的對(duì)比

為進(jìn)一步驗(yàn)證式（7）非線性渦激力數(shù)學(xué)模型的準(zhǔn)確性，基于四階Runge-Kutta 法利用MATLAB 自編程序通過式（7）確定的非線性渦激力數(shù)學(xué)模型反算出帶挑臂鋼箱主梁的渦振響應(yīng)，并與CFD 數(shù)值模擬計(jì)算的渦振位移時(shí)程進(jìn)行對(duì)比，如圖18（a）所示.從圖中可以看出，在振幅達(dá)到穩(wěn)定后，通過式（7）渦激力模型反算的渦振位移值與CFD 數(shù)值模擬計(jì)算的目標(biāo)值吻合良好.但另一方面，在渦振振幅達(dá)到穩(wěn)態(tài)前，渦振位移的反算值與目標(biāo)值有一定區(qū)別，但總體上差異較小.圖18（b）給出了由式（9）線性渦激力模型計(jì)算的位移值，由圖可知，在渦振振幅達(dá)到穩(wěn)態(tài)后，由線性渦激力模型得到的反算值相比目標(biāo)值要偏小7%左右.同時(shí)，在渦振振幅達(dá)到穩(wěn)態(tài)前，其渦振位移也明顯大于目標(biāo)值，由此可知，由經(jīng)典線性渦激力模型反算的主梁渦振位移值的精度不如非線性渦激力模型.同理，對(duì)于其他風(fēng)速情況，也有類似趨勢，在此不再贅述.

圖18 不同渦激力模型主梁位移反算值與目標(biāo)值對(duì)比Fig.18 Comparisons of displacements of the main girder inversely solved by different vortex-induced force models with the target values

此外，研究中還發(fā)現(xiàn)，不同風(fēng)速時(shí)兩種渦激力模型反算的渦振位移頻率值與目標(biāo)值較為吻合，說明兩種模型對(duì)振動(dòng)系統(tǒng)的氣動(dòng)剛度模擬較準(zhǔn)確.通過上述分析可知，非線性渦激力模型比線性模型的精度更高，反映出式（7）渦激力模型的準(zhǔn)確性.

3 渦振穩(wěn)態(tài)振幅計(jì)算公式

基于以上建立的帶挑臂鋼箱主梁渦激力的數(shù)學(xué)模型，以下討論其主梁的渦振穩(wěn)態(tài)振幅.主梁渦振時(shí)，振動(dòng)系統(tǒng)的動(dòng)力學(xué)控制方程可表示為：

式中：m為振動(dòng)系統(tǒng)的質(zhì)量；ζ為振動(dòng)系統(tǒng)的阻尼比；ω為系統(tǒng)的固有圓頻率.

在渦振鎖定區(qū)域，主梁的位移可近似為簡諧振動(dòng)，此時(shí)可假設(shè)其振動(dòng)位移為：

式中：A是振幅；ω1為鎖定區(qū)的振動(dòng)頻率；ψ是相位角.

根據(jù)式（11），得其渦振時(shí)的振動(dòng)速度為：

進(jìn)一步假設(shè)振幅A與相位角ψ是時(shí)間t的函數(shù)［15-16］，對(duì)式（11）求導(dǎo)并代入式（12）可得：

對(duì)式（12）求導(dǎo)，并聯(lián)立式（11）與式（10）可得：

聯(lián)立式（14）與式（13）可求解得到：

引入慢變函數(shù)的概念［17］，即假設(shè)在［0，2π］的一個(gè)整周期內(nèi)，函數(shù)A與ψ的變化都很小，此時(shí)函數(shù)值可用一個(gè)整周期內(nèi)的平均值來代替.令p=（ω1tψ），即可得：

將式（11）、式（12）、式（15）及式（7）代入式（17），并考慮式（15）或式（7）中的簡諧渦脫力項(xiàng)占比較小，在推導(dǎo)中予以忽略，由此可得：

求解式（18）的微分方程，可得：

式中：A0為初始位移；C1與C2為式（18）中的系數(shù)，其表達(dá)式如下：

對(duì)于渦振穩(wěn)態(tài)振幅，可令式（19）中的t=+∞，即可得到渦振穩(wěn)態(tài)振幅的表達(dá)式為：

根據(jù)上述分析，若將各參數(shù)代入式（22），可計(jì)算得到主梁模型渦振時(shí)的穩(wěn)態(tài)振幅A（+∞）=0.003 01 m；而由圖8（a）可以看出，CFD 計(jì)算出的渦振穩(wěn)態(tài)位移為0.002 69 m，兩者相差約11.9%.分析原因，這可能是由于當(dāng)前主梁模型的渦振振幅較小，從而導(dǎo)致相對(duì)誤差較大；另外，在推導(dǎo)過程中，忽略了簡諧渦脫力項(xiàng)，這也會(huì)部分影響計(jì)算結(jié)果的精度.但總體而言，通過式（22）計(jì)算的渦振穩(wěn)態(tài)振幅值與CFD 計(jì)算值大致相當(dāng).更為重要的是，從式（22）中可看出，渦振穩(wěn)態(tài)振幅隨著一階阻尼和二階剛度交叉項(xiàng)P12（為負(fù)數(shù)）絕對(duì)值的增大而減小，同時(shí)隨著一階阻尼項(xiàng)P10的增大而增大.

4 結(jié)論

本文對(duì)+3°風(fēng)攻角下帶挑臂鋼箱主梁斷面進(jìn)行了雙向流固耦合數(shù)值模擬計(jì)算，基于計(jì)算結(jié)果開展了該主梁渦激力數(shù)學(xué)模型的研究，可得出以下主要結(jié)論：

1）基于雙向流固耦合數(shù)值模擬方法對(duì)帶挑臂鋼箱主梁的渦振全過程進(jìn)行了數(shù)值模擬計(jì)算，并且數(shù)值模擬的渦振鎖定區(qū)間及渦振振幅值和節(jié)段模型風(fēng)洞試驗(yàn)結(jié)果吻合良好，表明文中的計(jì)算模型、網(wǎng)格設(shè)置及流固耦合方法等滿足精度要求.

2）基于已驗(yàn)證的數(shù)值模擬計(jì)算結(jié)果，發(fā)現(xiàn)帶挑臂鋼箱主梁的渦激力時(shí)程頻譜圖中，除最顯著的一階主頻外，還存在二階和三階倍頻，這種現(xiàn)象表明該主梁渦激力具有一定的非線性特性.在此基礎(chǔ)上，根據(jù)重構(gòu)的渦激力值與目標(biāo)值的殘差大小，建立了帶挑臂鋼箱主梁的非線性渦激力數(shù)學(xué)模型.

3）基于上述建立的渦激力數(shù)學(xué)模型，反算出的渦振位移時(shí)程與CFD 計(jì)算的位移時(shí)程總體吻合良好，進(jìn)一步驗(yàn)證了所建立渦激力模型的準(zhǔn)確性；同時(shí)基于慢變函數(shù)理論，推導(dǎo)了該主梁斷面渦振穩(wěn)態(tài)振幅的表達(dá)式，并發(fā)現(xiàn)該主梁渦振穩(wěn)態(tài)振幅隨著一階阻尼和二階剛度交叉項(xiàng)P12絕對(duì)值的增大而減小，而隨著一階阻尼項(xiàng)P10的增大而增大.

本文目前只研究了一種形式的帶挑臂鋼箱主梁斷面的渦激力模型，未來將對(duì)更多類似斷面進(jìn)行研究，以驗(yàn)證該渦激力模型的準(zhǔn)確性與適用性.