基于Küssner函數的不同氣動導納模型對大跨橋梁抖振響應的影響

張志田， 陳添樂， 吳長青

(湖南大學 風工程試驗研究中心，長沙 410082)

橋梁抖振分析有時域和頻域兩種算法，頻域分析具有簡單快速的優點，但只適用于線性問題，只能得到結構響應的統計特性。時域分析能考慮各種非線性影響，并能得到抖振位移時程，但不便于氣動導納的靈活應用。抖振時域分析分為以下三個環節：橋梁脈動風場模擬、抖振力的表達、抖振位移計算。脈動風場模擬一般采用諧波合成法，國內外許多學者已對此做過詳盡研究[1-4]。抖振力的關鍵問題是氣動導納函數的確定。自20世紀60年代Davenport將Sears函數和Liepmann的機翼抖振理論應用于橋梁抖振分析至今[5-6]，橋梁顫抖振計算通常采用Davenport準定常模型、氣動導納、以及Scanlan建議自激力表達式[7]。對于氣動導納，文獻中通常要么取值為1(即不考慮抖振力的1非定常特性)，要么取Sears函數或者試驗導納[8-11]。但是三種氣動導納計算結果到底有多大的差別是一個值得深入研究的問題。在目前抖振時域計算中，一般是采用等效風譜法進行抖振力的計算[12]，即把氣動導納乘以風譜看成一個等效風譜，再用得到的等效風譜進行脈動風場模擬。等效風譜法是通過修正風譜來間接考慮抖振力的非定常性，當導納為試驗值時，需要分別對升力、阻力、扭矩的等效風速時程進行模擬。故本文放棄采用等效風譜法，而是從氣動力演變特性出發采用Küssner函數進行抖振力模擬。并在此基礎上分析不同氣動導納模型引起的抖振響應差別。

1 非定常抖振力的卷積表達

以水平速度U飛行的斷面，穿過幅值為w0的豎向階躍陣風時，其氣動力隨時間演變的公式為[13]

(1)

(2)

式中：a0,ai,di為常數。

文獻[14-15]中指出式(2)中的a0項應該等于1，本文抖振力表達式最終是要和Scanlan建議抖振力模型進行功率譜等效，故不詳細討論其所代表的物理意義。

在線性疊加原理滿足的前提下，任意分布的豎向脈動風w(t)引起的非定常氣動升力可表示為

(3)

式中：ψ′(σ)為Küssner函數對時間t的導數。

式(3)中只考慮了豎向脈動風的影響，在橋梁風工程中，同時考慮水平向和豎向脈動風所引起的非定常氣動力為

(4)

(5)

(6)

式中：C′L,C′D,C′M分別為升力系數、阻力系數和升力矩系數對攻角的導數； 函數ψFx(F=L,D,M;x=u,w)表示x方向的單位數值脈動風引起的每延米F氣動力的瞬態演變，都為Küssner類型的函數；ψ′Fx(F=L,D,M;x=u,w)為Küssner函數對時間t的導數；ψFx表達式如下

(7)

2 Küssner函數的擬合

斷面氣動導納函數已知的情況下，引入氣動導納修正的Scanlan建議抖振力模型如下[16-17]

(8)

(9)

(10)

式中：CL,CD,CM分別為升力系數、阻力系數和升力矩系數；χLu,χLw,χDu,χDw,χMu,χMw為6個氣動導納。

根據式(8)～(10)，求得抖振力功率譜如下

(11)

(12)

(13)

式中：SL,SD,SM分別為升力、阻力、升力矩的功率譜；Suu,Sww分別為水平和豎向的脈動風功率譜；Suw,Swu為脈動風互功率譜；式中*號表示求復共軛。

同樣，對本文中采用Küssner函數表示的抖振力表達式(4)～(6)求得抖振力功率譜如下

(14)

(15)

(16)

式中的-號表示求傅里葉變換。

根據兩種抖振力表達式功率譜相等的原則，忽略互功率譜的影響，分別比較式(11)與式(14)，式(12)與式(15)，式(13)與式(16)，可以得到如下關系

(17)

(18)

(19)

(20)

(21)

(22)

(23)

式中：ω位脈動風的圓頻率。

將式(23)代入式(17)～(22)中得到各Küssner函數參數與氣動導納的關系后，可根據已知的氣動導納確定Küssner函數參數。這是一個最優化問題，可以通過以下六個函數的極小值問題求得所需結果

(24)

(25)

(26)

(27)

(28)

(29)

式中： 求和符號上標n為氣動導納試驗點(不同的折算頻率K)的數目。

根據式(24)～(29)就能識別得到與水平以及豎向脈動風相關的升力、阻力和升力矩6個Küssner函數的待定參數，需要說明的是參數擬合過程中可能存在多個折算頻率點，即存在很多擬合點，為此作者使用遺傳算法進行參數擬合[18]，在設定初始種群、適應度函數以及控制參數后，遺傳算法會對每一次進化后的種群進行計算，使式(24)～式(29)的值最小的一些個體稱為優異個體，下一次進化時優異個體會得以保留，其它個體在會淘汰或者變異，將進化后的種群代入式(24)～式(29)重新計算。隨著種群的進化，最終優異個體會越來越多，直至產生滿足設定條件的個體。

3 抖振時域分析

求出升力、阻力和升力矩的Küssner函數表達式后，可采用式(4)～式(6)對結構進行動力有限元抖振分析。但由于式(4)～式(6)中存在與脈動風速相關的時間積分，直接使用其求解抖振力的計算效率非常低。參考文獻[14-15]的時域自激力遞推表達式可得到抖振力的遞推表達如下：

首先，將抖振力分解成脈動風u、v兩部分的貢獻之和

(30)

(31)

(32)

式中：Fbu(t),Fbw(t)(F=L,D,M)分別為順風向和豎風向脈動風速引起的抖振力(升力、阻力和升力矩)。以式(30)中順風向脈動風速引起的抖振升力Lbu(t)為例進行如下推導，對積分項中進行變量代換后得

(33)

式中：

(34)

由式(33)可知，對于任意時刻的抖振力，第一項只與當前時刻的風速有關，而第二項則為時間積分項，與整個時間歷程風速有關，故根據t時刻的抖振力表達式遞推得到t+Δt時刻的抖振力，只需要求式(34)的遞推表達即可，該式的遞推表達如下

(35)

同理可得到其它幾個抖振力的遞推表達式

(36)

(37)

(38)

(39)

(40)

4 數值算例

以某特大懸索橋的初步設計方案為例，計算其在0°風攻角下的抖振位移響應。該橋為主跨920 m的單跨扁平鋼箱梁懸索橋，主塔采用混凝土澆筑，共設75對吊桿，橋址位于我國西南部山區峽谷，受脈動風影響較大。該橋立面圖和箱梁斷面，如圖1與圖2所示。

圖1 懸索橋立面圖(cm)Fig.1 The elevation of the suspension bridge(cm)

圖2 橋梁箱梁斷面(cm)Fig.2 Section of the box girder(cm)

4.1 三分力系數和氣動導納試驗值

采用幾何縮尺比為1∶50的主梁節段模型進行風洞試驗。模型的長度為1.54 m，寬度為0.64 m。模型本身具有足夠大的剛度，測力試驗以及氣動導納試驗中無出現明顯變形。三分力試驗采用高頻測力天平測得試驗風速下不同風攻角的靜力三分力，通過計算得到該斷面的三分力系數以及三分力系數對攻角的導數。氣動導納測試時采用水平格柵產生被動紊流。三分力試驗以及氣動導納試驗的風洞布置分別見圖3(a)與3(b)。通過風洞試驗測得0°攻角時的三分力系數為CL=-0.190 7，CD=0.029，CM=0.008 5，三分力系數對攻角的導數為C′L=3.095 1，C′D=0.280 2，C′M=1.057 1。

圖3 風洞試驗Fig.3 Wind tunnel tests

氣動導數的識別方法主要有等效導納法、交叉譜法以及自功率譜-交叉譜總體最小二乘法。等效導納法是基于抖振力功率譜識別導納，需要假設順風向導納與豎風向導納相等[19]；交叉譜法通過求氣動力和脈動風速的交叉譜來識別氣動導納[20]，自功率譜-交叉譜總體最小二乘法是結合上述兩種方法，采用最小二乘法識別氣動導納[21]。此外，對于存在分享流的鈍體橋梁斷面，氣動導納并非僅僅由氣動外形以及頻率決定，同時也受來流的湍流特性的影響，即不同的風場特性會產生不同的氣動導納[22]。本文關注的側重點不是氣動導納的識別方法以及不同識別方法所帶來的差異，因此這里采用最簡單也最穩定的等效導納法識別法。等效導納法假設χL=χLu=χLw,χD=χDu=χDw,χM=χMu=χMw[19], 將其代入式(11)～式(13)中，忽略互功率譜的影響可得

(41)

(42)

(43)

圖4 等效氣動導納試驗值Fig.4 The test value of equivalent aerodynamic admittance

4.2 Küssner函數參數擬合

在Küssner函數的擬合過程中，本例中指數項取三項，需要擬合的參數為7個，采用遺傳算法進行參數擬合時種群規模取為5 000，每一代的精英數目為300，交叉概率取0.8，變異率取0.1，停止條件為容許誤差小于0.000 1。表1和表2分別列出了Sears函數和試驗氣動導納的各Küssner函數參數擬合值。

表1 Sears函數的Küssner函數參數擬合值Tab.1 Fitted value of Sears function of Küssner function

表2 試驗導納的Küssner函數參數擬合值Tab.2 Fitted value of test admittance of Küssner function

將擬合得到的Küssner函數代入式(17)～式(22)中，可以反算出氣動導納擬合值，圖5給出了Sears函數擬合值與Sears函數的比較，圖6給出了試驗氣動導納擬合值與試驗值的比較。從圖5和圖6可以看出，不論是Sears函數還是試驗氣動導納，反算得到的氣動導納與目標值都具有很高的吻合度。因此，本文的抖振力模型與引入氣動導納修正的Scanlan建議模型兩者的抖振力功率譜是等效的。

圖5 Sears函數擬合值與Sears函數的比較Fig.5 Comparison of fitted Sears with Sears function

圖6 試驗導納擬合值與試驗導納的比較Fig.6 Comparison of fitted equivalent aerodynamic admittance with test values

4.3 脈動風場模擬

本文的風場模擬中，順風向和豎風向風譜分別采用Kaimal譜、Panofsky譜[23-24]，風場模擬采用諧波合成法，具體過程已有學者做過詳細研究，在此不再贅述。需要說明的是本文主要是研究不同氣動導納對橋梁抖振響應的影響，因此在合成功率譜密度矩陣時，無需討論衰減因子λ取值對脈動風速時程的影響，本文中取值為λ=7，且風速時程取為600 s。表3給出了風場模擬的主要參數。圖7、圖8分別為跨中處模擬的順風向和豎風向脈動風速時程，圖9、圖10分別為對應的脈動風速功率譜和目標功率譜的比較。

表3 風場模擬的主要參數Tab.3 Main parameters of the wind field simulation

圖7 跨中處順風向脈動風速時程Fig.7 Time history of fluctuating wind velocity of the mid-span section

圖8 跨中處豎風向脈動風速時程Fig.8 Time history of vertical wind velocity of the mid-span section

圖9 跨中處順風向脈動風功率譜密度Fig.9 The power spectral density of fluctuating wind of the mid-span section

4.4 抖振力計算

“4.3”中已經說明本文主要目的是研究不同氣動導納對橋梁抖振響應的影響，因此在動力有限元計算時，荷載只考慮了抖振力。

圖10 跨中處豎風向脈動風功率譜密度Fig.10 The power spectral density of vertical wind of the mid-span section

對氣動導納等于1即不考慮氣動導納的情況，直接采用準定常模型計算得到抖振力時程；對于考慮氣動導納修正的情況，采用本文推導的抖振力遞推表達式計算。將Sears函數的Küssner函數參數擬合值與模擬脈動風速時程帶入本文抖振力遞推表達式中，得到考慮Sears函數修正的抖振力時程；同樣，將試驗氣動導納的Küssner函數參數擬合值與模擬脈動風速時程帶入遞推表達式中，得到考慮試驗氣動導納修正的抖振力時程。至此，已獲得考慮不同氣動導納影響抖振計算所需的三組抖振力。跨中處不同氣動導納計算得到的抖振力時程，如圖11所示。

4.5 本文küssner表達式抖振位移計算結果與分析

圖12為跨中處不同氣動導納函數計算得到的抖振位移時程，根據圖10可知，不同氣動導納計算得到的抖振位移響應時程在變化趨勢上大體一致，但是響應大小有明顯差別。對于本例中3種抖振位移響應，考慮氣動導納修正的抖振位移遠小于不考慮導納修正的位移；導納為1時計算得到的抖振位移最大，試驗氣動導納修正計算得到的抖振位移最小，Sears函數修正的結果位于導納為1和試驗導納計算結果之間。

圖13更清楚地反應了不同氣動導納引起的抖振位移的差別。①對比導納為Sears函數和導納為1時的情況，豎向位移均方根前者比后者小了32.4%；扭轉位移均方根前者比后者小了45.2%；側向位移均方根前者比后者小了50.1%。②對比導納為試驗值和導納為1時的情況，豎向位移均方根前者比后者小了35.2%；扭轉位移均方根前者比后者小了66.3%；側向位移均方根前者比后者小了65.8%。③本例中三種抖振位移響應均方根，考慮氣動導納修正的結果均小于氣動導納為1的情況，且氣動導納為試驗值的結果最小。④Sears函數對不同的抖振位移響應均方根影響程度不一樣，Sears函數對側向位移響應的均方根影響最大，對豎向位移響應的均方根影響最小。⑤試驗導納對側向和扭轉位移響應的均方根影響較大，對豎向位移響應的均方根影響最小。

圖11 跨中處的抖振力時程Fig.11 Time history of buffeting force of the mid-span section

圖12 跨中處的抖振位移時程(küssner表達)Fig.12 Time history of buffeting response of the mid-span section(by küssner expression)

圖13 跨中處的抖振響應均方根Fig.13 RMS of buffeting response of the mid-span section

圖14為跨中處不同氣動導納函數計算得到的位移功率譜密度，本算例中豎向位移功率譜第一個峰值橫坐標為0.171 8，與結構一階正對稱豎彎自振頻率(0.174 5)相對應，扭轉位移功率譜第一個峰值橫坐標為0.403 3，與結構一階正對稱扭轉自振頻率(0.401 9)相對應，側向位移功率譜第一個峰值橫坐標為0.083，與結構一階正對稱側彎自振頻率(0.082 9)相對應。位移功率譜密度最大值均出現在第一個峰值處。

圖14 跨中處抖振位移功率譜密度Fig.14 The power spectral density of buffeting response of the mid-span section

從圖14可以發現： ①抖振位移響應主要是受若干階低頻模態的控制，尤其是受基頻的影響最大。②在主要頻率范圍內，導納函數為1時位移響應功率譜密度最大，Sears其次，試驗導納的位移功率譜密度最小。這與抖振位移時程的結論是一致的。

4.6 本文küssner表達式法與等效風譜法的比較

通常采用等效風譜法進行考慮氣動導納修正的抖振時程計算，即把氣動導納函數乘以風譜看成等效風譜，然后用等效風譜進行脈動風場模擬。本小節中首先得到考慮sears函數和試驗導納函數修正的等效風譜，再采用和“4.3”節中相同的參數進行脈動風場模擬。對于試驗氣動導納修正的情況，由于升力、扭矩和阻力的導納函數不同，需分別模擬升力、扭矩、阻力的順風向和豎風向共6組脈動風速時程。

圖15為等效風譜法計算得到的跨中抖振位移時程，與圖14中的Küssner方法所得結果相比可以發現：①定性結果一致，即考慮sears函數修正的抖振位移小于不考慮修正結果，考慮試驗導納修正的抖振位移小于sears函數修正結果。②兩者計算得到的隨機抖振位移的變化范圍基本一致。③本文中采用Küssner函數時采用同一組脈動風速時程，不同導納計算的抖振位移隨時間變化具有同步性；而采用等效風譜法時，不同的導納函數和抖振力均需要單獨模擬脈動風時程，因而不同導納計算的抖振位移隨時間變化無同步性。

圖16為分別采用Küssner表達式和等效風譜法計算得到考慮sears函數修正的跨中抖振位移時程，可見雖然兩者抖振位移隨時間變化并不一致，但是兩者的變化范圍基本相同。

圖15 跨中處的抖振位移時程(等效風譜法)Fig.15 Time history of buffeting response of the mid-span section(by the equivalent wind spectrum)

圖16 跨中處Küssner表達式和等效風譜法計算得到的考慮sears函數修正的抖振位移時程Fig.16 Time history of buffeting response of the mid-span section calculated by Küssner expression and equivalent wind spectrum method which sears function is modified

圖13比較了本文Küssner表達式和等效風譜法計算得到的抖振響應均方根。①導納為sears函數時，豎向位移兩者計算相差13.1%，扭轉位移相差9.1%，側向位移相差 11.8%。②導納為試驗值時，豎向位移兩者計算相差16.3%，扭轉位移相差2.7%，側向位移相差 4.3%。兩種方法在計算豎向位移均方根時相差較大，作者認為主要原因有兩方面：其一是有限長度的風場模擬時存在個體差異；其二是兩種方法形成的抖振力相關譜不等效。

5 結 論

扁平鋼箱梁是大跨度橋梁常用的加勁梁型式。由于其良好的氣動外形與平板具有一定的相似性，對這類橋梁進行抖振分析時常采用源于古典機翼理論的Sears氣動導納函數代替斷面的實際氣動導納函數，或者偏于安全不考慮氣動導納函數。橋梁抖振分析的傳統方法有頻域與時域兩類方法。頻域法可方便地考慮氣動導納，但忽略了結構的非線性；時域法能全面考慮非線性，但不便于氣動導納的應用。本文基于Küssner函數時域模型，在時域內研究了這三種不同的導納選取方案對千米級橋梁抖振分析結果的影響，并得到以下結論：

(1) 本文提出Küssner階躍函數法可將氣動導納靈活地從頻域轉換到時域后進行非線性動力有限元分析。

(2) 與等效風譜法相比，本文的Küssner函數時域抖振模型只需模擬一組脈動風速時程，提高了計算效率。此外，應用于實測風速時程時，避免了根據時程反算風譜、風譜等效后再模擬時程的繁瑣過程。

(3) 基于Sears函數氣動導納與試驗氣動導納得到的豎向響應基本一致，但扭轉與側向響應相差明顯。以扭轉響應為例，前者在后者的基礎上高出約63%。顯然，即使是扁平鋼箱梁，采用源于機翼理論的Sears函數進行抖振分析不合適。