典型帶挑臂鋼箱主梁渦激力模型研究

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

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

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

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

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

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

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

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

1.1 工程簡介

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

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

1.2 試驗概況

本橋的渦振試驗是在長沙理工大學風工程與風環(huán)境研究中心的高速試驗段中進行（如圖2 所示），采用Cobra Probe 三維脈動風速儀測試來流風速，采用激光位移計測試主梁節(jié)段模型的位移時程.彈性懸掛的節(jié)段模型縮尺比為1∶50，風速比為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 放置于風洞中的試驗模型Fig.2 Testing model in the wind tunnel

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

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

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

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

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

圖5 全局及細部網(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é)果驗證

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

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

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

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

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

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

2.1 主梁渦激力的提取

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2 渦激力數(shù)學模型的建立及驗證

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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