劉修平， 劉煥舉， 趙 越， 李光玲

(1.天水師范學(xué)院 土木工程學(xué)院，甘肅 天水 741000;2.河北工程大學(xué) 土木工程學(xué)院，河北 邯鄲 056038;3.長(zhǎng)安大學(xué) 公路學(xué)院，西安 710064)

顫振是一種發(fā)散性氣動(dòng)力自激振動(dòng)，是橋梁風(fēng)致振動(dòng)當(dāng)中最危險(xiǎn)的振動(dòng)形式。自Scanlan引入顫振導(dǎo)數(shù)描述鈍體橋梁自激力以來(lái)，分離流顫振機(jī)理在橋梁顫振分析應(yīng)用中日趨完善[1-3]。然而隨著大跨纜索承重橋梁在跨越江海、深切峽谷的國(guó)家公路和鐵路交通控制性工程廣泛使用，結(jié)構(gòu)非線性、氣動(dòng)非線性等非線性因素對(duì)橋梁穩(wěn)定性的影響也日益突出，風(fēng)致顫振穩(wěn)定性成了大跨橋梁設(shè)計(jì)施工過(guò)程中需考慮的關(guān)鍵性問(wèn)題。由于Scanlan線性自激力模型不能很好的描述考慮各類(lèi)非線性因素在顫振中的影響，在橋梁抗風(fēng)研究日益精細(xì)的情況下，自激力時(shí)域模型研究成為了風(fēng)致顫振分析的熱點(diǎn)[4-6]。

自激力的時(shí)域表達(dá)式有兩種方法：其一是Scanlan提出的Wagner在航空中提出的經(jīng)典階躍函數(shù)表達(dá)法，其二是Buther & Lin提出的單位脈沖響應(yīng)函數(shù)融入Roger有理函數(shù)的表達(dá)方法。張志田等[7]分別研究了兩種自激力時(shí)域表達(dá)式的瞬態(tài)特性及極限特性，提出了通過(guò)限定有理函數(shù)表達(dá)式若干參數(shù)克服瞬態(tài)特性模擬失真的方法；Chen等[8]基于有理函數(shù)表達(dá)的自激氣模型，提出了橋梁大跨度橋梁結(jié)構(gòu)顫振分析的時(shí)域分析流程及特征值求解方法；郭增偉等[9-10]通過(guò)顫振導(dǎo)數(shù)的有理函數(shù)近似式得到自激力脈沖響應(yīng)函數(shù)在橋梁顫振分析中的表現(xiàn)方式和求解方案，提出了顫振分析中初始響應(yīng)失真及積分時(shí)間步長(zhǎng)設(shè)定的解決方法；伍波等[11-12]發(fā)現(xiàn)相比經(jīng)典階躍函數(shù)顫振分析方法，在考慮靜風(fēng)效應(yīng)的大跨度橋梁結(jié)構(gòu)顫振特性分析中，有理函數(shù)自激力表達(dá)式模型可以有效減少計(jì)算風(fēng)速下的迭代量，顯著提高計(jì)算效率。由于橋梁顫振非線性問(wèn)題分析和方程求解的復(fù)雜性，一直以來(lái)，橋梁顫振時(shí)域分析多是基于C/C++、Fortran開(kāi)發(fā)的專(zhuān)用程序進(jìn)行[13-17]。近些年，隨著大型通用有限元分析程序在橋梁工程中具有廣泛的應(yīng)用，基于通用有限元軟件ANSYS的顫振頻率分析方法日趨成熟，直接采用ANSYS進(jìn)行橋梁顫振時(shí)域分析的研究較少。華旭剛等[18]基于ANSYS中MATRIX27矩陣單元建立了橋梁結(jié)構(gòu)全模態(tài)顫振分析模型，提出以調(diào)整風(fēng)速和振動(dòng)頻率求解動(dòng)力系統(tǒng)各階復(fù)模態(tài)特性的顫振頻域分析方法，曾憲武等[19]基于MATRIX27矩陣單元對(duì)自激力以單元?dú)鈩?dòng)阻尼矩陣和單元?dú)鈩?dòng)剛度矩陣提出了基于ANSYS的大跨橋梁抖振分析方法；為避免在使用過(guò)程中需假定振動(dòng)頻率的步驟，謝煉等[20]利用了ANSYS內(nèi)置的重啟動(dòng)算法，迭代計(jì)算當(dāng)前時(shí)步真實(shí)振動(dòng)狀態(tài)，提出了橋梁顫振時(shí)域計(jì)算的重啟動(dòng)法。

為避免MATRIX27矩陣單元輸入振動(dòng)頻率的誤差干擾，同時(shí)簡(jiǎn)化計(jì)算時(shí)步內(nèi)的自激力反復(fù)迭代確定的重復(fù)過(guò)程。本文基于傳統(tǒng)有理函數(shù)自激力表達(dá)式推導(dǎo)處理，使得自激力表達(dá)式中的“記憶效應(yīng)”項(xiàng)僅與上一時(shí)刻的橋面運(yùn)動(dòng)狀態(tài)相關(guān)聯(lián)，繼而重新形成有理函數(shù)自激力有表達(dá)式的質(zhì)量項(xiàng)、剛度項(xiàng)、阻尼項(xiàng)和時(shí)間歷程項(xiàng)，建立基于MATRIX27矩陣單元的自激力時(shí)域模擬方法，并通過(guò)經(jīng)典平板顫振分析成功驗(yàn)證了所提方法的準(zhǔn)確性。

1 橋梁斷面自激力

1.1 自激力有理函數(shù)表達(dá)

(1)

p=K(ξ+i)，ξ為橋梁結(jié)構(gòu)阻尼比一般較小，可以忽略；其他參數(shù)含義如表1所示。

表1 8個(gè)顫振導(dǎo)數(shù)的參與振動(dòng)方向Tab.1 Directions of eight aerodynamic derivatives

Scanlan自激力列式和振動(dòng)圓頻率關(guān)聯(lián)，不具備直接進(jìn)行時(shí)域分析計(jì)算的條件。Lin提出可以完全在時(shí)域中描述風(fēng)致振動(dòng)中產(chǎn)生的耦合氣動(dòng)自激力的關(guān)系式，將之表達(dá)為

(4)

根據(jù)式(3)，有理函數(shù)推廣用于近似橋梁斷面自激力的氣動(dòng)矩陣Q，形式為

(5)

式中：C1，…Cl+3均與頻率無(wú)關(guān)的氣動(dòng)力系數(shù)；C1為自激力的氣動(dòng)剛度項(xiàng)；C2為氣動(dòng)阻尼項(xiàng)；C3為氣動(dòng)質(zhì)量項(xiàng)，其值一般較小，可以忽略；第4項(xiàng)以后都為自激力的記憶效應(yīng)，與結(jié)構(gòu)前一時(shí)刻和當(dāng)前時(shí)刻的運(yùn)動(dòng)有關(guān)；λl為大于0的非定常項(xiàng)衰減常數(shù)。

1.2 氣動(dòng)力系數(shù)識(shí)別

顫振導(dǎo)數(shù)的擬合對(duì)象為高度非線性的方程組，為避免大量的復(fù)數(shù)運(yùn)算過(guò)程，令上述方程的實(shí)部和虛部相等，以升力的計(jì)算為例，則可以得到如下的關(guān)系式

(6)

該組方程式對(duì)應(yīng)7個(gè)未知量，擬合過(guò)程需要參數(shù)共享，同時(shí)也要避免因初始值的影響陷入局部最優(yōu)的問(wèn)題。對(duì)于二維彎扭耦合顫振導(dǎo)數(shù)擬合共需識(shí)別28個(gè)氣動(dòng)力系數(shù)，因此，優(yōu)化迭代算法的選擇至關(guān)重要。

以理想平板顫振導(dǎo)數(shù)擬合為例，在最大迭代次數(shù)2 000、收斂誤差為1×10-10的設(shè)定下，分別采用常用的麥夸特法、擬牛頓法、遺傳算法、模擬退火算法、擬牛頓+通用全局優(yōu)化算法求解上述方程組(2)式，求解結(jié)果和過(guò)程參數(shù)對(duì)比如表2所示。

表2 不同算法求解結(jié)果比較Tab.2 Comparison of the results of different algorithms

由表2可知，在顫振導(dǎo)數(shù)擬合的實(shí)際計(jì)算中，求解此類(lèi)方程組，麥夸特法、擬牛頓法和擬牛頓+通用全局優(yōu)化算法都能快速的得到收斂解，而遺傳算法和模擬退火算法的計(jì)算效率較低，達(dá)到限定迭代步數(shù)后仍未收斂，同時(shí)對(duì)非定常項(xiàng)衰減常數(shù)λ計(jì)算值出入較大，甚至出現(xiàn)負(fù)值，顯然，考慮卷積積分Φsel的存在，λ應(yīng)是大于0的實(shí)數(shù)。因此，在氣動(dòng)力系數(shù)識(shí)別過(guò)程中應(yīng)優(yōu)先采用擬牛頓+通用全局優(yōu)化算法。

2 基于MATRIX27的自激力時(shí)域模擬

2.1 MATRIX27矩陣單元的組集

(7)

(8)

(9)

MATRIX27矩陣單元無(wú)幾何外形特性，其運(yùn)動(dòng)學(xué)響應(yīng)用剛度、阻尼或者質(zhì)量系數(shù)來(lái)指定，該矩陣單元連接2個(gè)節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)有6個(gè)自由度。在通過(guò)式(8)、(9)集成式(7)中的各矩陣后，將矩陣內(nèi)的元素以實(shí)常數(shù)的方式賦予MATRIX27單元，并將其添加至結(jié)構(gòu)單元的節(jié)點(diǎn)處，如圖3所示，時(shí)間歷程項(xiàng)作為外荷載形式輸入，從而形成時(shí)域自激力橋梁顫振分析有限元模型。

2.2 非線性項(xiàng)處理方法

時(shí)域自激力計(jì)算的難點(diǎn)在于非線性項(xiàng)的荷載受當(dāng)前的時(shí)步ti和上一時(shí)步ti-1節(jié)點(diǎn)運(yùn)動(dòng)狀態(tài)共同影響，而商業(yè)有限元軟件ANSYS的瞬態(tài)分析是基于隱式算法，無(wú)法推算出當(dāng)前時(shí)步的荷載。雖然可以利用重啟動(dòng)實(shí)現(xiàn)了當(dāng)前時(shí)步荷載的計(jì)算，但重啟動(dòng)每個(gè)計(jì)算時(shí)步需迭代確定，在通用有限元程序中計(jì)算步數(shù)相對(duì)較多，耗時(shí)較長(zhǎng)，對(duì)于大跨橋梁風(fēng)致抖振分析而言效率并不高。因此，真正意義地實(shí)現(xiàn)基于MATRIX27矩陣單元的ANSYS時(shí)域顫振計(jì)算的前提就是將有理函數(shù)自激力表達(dá)式重新改寫(xiě)，具體如下：

(1) 將非線性項(xiàng)重新改寫(xiě)，以單位長(zhǎng)度上扭轉(zhuǎn)運(yùn)動(dòng)對(duì)扭轉(zhuǎn)的影響為例。

(10)

(11)

引入，

那么，

(12)

利用Φsel(ti)和Φsel(ti-1)的關(guān)系可以得到

(13)

(2) 形成自激力有理函數(shù)表達(dá)的質(zhì)量項(xiàng)、剛度項(xiàng)、阻尼項(xiàng)和時(shí)間歷程項(xiàng)。

在非線性項(xiàng)改寫(xiě)的基礎(chǔ)上，那么，式(10)和式(11)可以表達(dá)為

(14)

(15)

令

(16)

2.3 時(shí)程分析實(shí)現(xiàn)過(guò)程

橋梁顫振時(shí)域分析過(guò)程中，首先要顫振分析有限元模型施加一個(gè)初始激勵(lì)，由于自激振動(dòng)系統(tǒng)的特征量，如頻率、衰減率等與初始激勵(lì)無(wú)關(guān)，在顫振分析時(shí)可對(duì)結(jié)構(gòu)施加初始激勵(lì)，激勵(lì)形式可為集中力或者位移等，在ANSYS中的具體實(shí)現(xiàn)過(guò)程為：

(1) 假定計(jì)算風(fēng)速。

(2) 施加任意初始荷載激勵(lì)。

(3) 計(jì)算第ti時(shí)刻的位移、速度和加速度向量。

(4) 計(jì)算ti+1時(shí)刻的荷載向量。

(5) 計(jì)算ti+1時(shí)刻的位移、速度和加速度向量。

(6) 依次重復(fù)步驟(4)、(5)，直到設(shè)置時(shí)間步完成。

(7) 計(jì)算此風(fēng)速下結(jié)構(gòu)振動(dòng)時(shí)程響應(yīng)的對(duì)數(shù)衰減率r，若r>0增加風(fēng)速返回步驟(2)，若r<0減小風(fēng)速返回步驟(2)，若r=0，則為顫振臨界風(fēng)速，程序結(jié)束。

在步驟(7)中的顫振臨界風(fēng)速確定時(shí)，可采用折半法查找，首先確定顫振臨界風(fēng)速大致區(qū)間[Ul,Uh]，計(jì)算風(fēng)速取Um=(Ul+Uh)/2，再進(jìn)行迭代計(jì)算。同時(shí)，由于ANSYS中的APDL語(yǔ)言編程的限制，節(jié)點(diǎn)振動(dòng)時(shí)程響應(yīng)對(duì)數(shù)衰減率的計(jì)算可由MATLAB中的峰值讀數(shù)法函數(shù)直接進(jìn)行，通過(guò)增加最小間隔條件，可以過(guò)濾掉單周期內(nèi)的干擾極值，更準(zhǔn)確的判斷振動(dòng)是否衰減，在程序中迭代框架具體見(jiàn)圖1。

圖1 顫振時(shí)域分析實(shí)現(xiàn)框架Fig.1 The implementation framework of time-domain flutter analysis

3 算例驗(yàn)證

3.1 理想平板斷面的簡(jiǎn)支梁

為驗(yàn)證所提時(shí)域自激力計(jì)算的準(zhǔn)確性，以理想平板的顫振分析為例，并將計(jì)算結(jié)果和理論解對(duì)比分析。該簡(jiǎn)支梁長(zhǎng)度為300 m，橋面寬度40 m，材料彈性模量取2.1×1012Pa，橫向抗彎慣性取8.57 m4，豎向抗彎慣矩取1.0×104m4，扭轉(zhuǎn)慣矩取0.386 m4，質(zhì)量取2×104kg/m，質(zhì)量慣矩取4.5×106kg·m，空氣密度取1.225 kg/m3，簡(jiǎn)支梁前10階的自振特性如表3所示，理想平板顫振導(dǎo)數(shù)如圖2所示。

(a) H*

(b) A*圖2 理想平板的顫振導(dǎo)數(shù)Fig.2 Aerodynamic derivatives of ideal plate

表3 簡(jiǎn)支梁前10階自振特性Tab.3 The first ten order frequencies of simply supported beam bridge

在MATRIX27矩陣單元的顫振分析模型建立過(guò)程中，首先將簡(jiǎn)支梁劃分為30個(gè)單元，將質(zhì)量等效到節(jié)點(diǎn)上，梁端節(jié)點(diǎn)約束扭轉(zhuǎn)自由度，引入29個(gè)氣動(dòng)剛度矩陣、29個(gè)氣動(dòng)阻尼矩陣和29個(gè)氣動(dòng)質(zhì)量矩陣，顫振分析有限元模型如圖3所示，模型共包括單元210個(gè)，節(jié)點(diǎn)122個(gè)，基于擬牛頓+通用全局優(yōu)化算法對(duì)理想平板顫振導(dǎo)數(shù)擬合后得到的氣動(dòng)力系數(shù)如表4所示。

圖3 簡(jiǎn)支梁顫振分析有限元模型Fig.3 Flutter analysis FE model of simply supported beam bridge

表4 理想平板氣動(dòng)力系數(shù)Tab.4 Aerodynamic coefficients of ideal plate

從表4中可見(jiàn)：φl(shuí)α和φmh的氣動(dòng)質(zhì)量項(xiàng)足夠小，可以忽略不計(jì)，而φl(shuí)h和φmα的擬合值并非如此，因此，在進(jìn)行氣動(dòng)導(dǎo)數(shù)擬合時(shí)，不考慮氣動(dòng)質(zhì)量項(xiàng)可能使得擬合結(jié)果出現(xiàn)偏差，可能影響顫振臨界風(fēng)速的計(jì)算結(jié)果。

簡(jiǎn)支梁在顫振發(fā)生及前后時(shí)的跨中豎向位移和扭轉(zhuǎn)位移的響應(yīng)時(shí)程如圖4所示，由圖可見(jiàn)：隨著風(fēng)速的增加，簡(jiǎn)支梁跨中位置振動(dòng)時(shí)程曲線由衰減到逐步發(fā)散，在顫振臨界風(fēng)速處，簡(jiǎn)支梁維持等幅振動(dòng)。

(a) 豎向位移時(shí)程

(b) 扭轉(zhuǎn)位移時(shí)程圖4 顫振前后簡(jiǎn)支梁跨中位移響應(yīng)Fig.4 Displacement responses of mid-span of simply supported beam bridge before and after flutter

表5為本文方法與其他方法計(jì)算結(jié)果的比較，由表可知：本文所提方法的計(jì)算結(jié)果與理論解吻合良好，誤差不足1%。

若不重新改寫(xiě)非線性Φsel，直接由Newmark-β進(jìn)行時(shí)程求解上述自激力模型時(shí)，則會(huì)出現(xiàn)文獻(xiàn)[9]中的結(jié)果，結(jié)構(gòu)的顫振臨界風(fēng)速值受時(shí)間步影響較大，這是由于計(jì)算時(shí)步的真實(shí)荷載受當(dāng)前和上一時(shí)刻橋梁運(yùn)動(dòng)狀態(tài)共同的影響，直接遞推得到的自激力在相位差和大小上與實(shí)際值存在一定的差距。為研究時(shí)間步長(zhǎng)對(duì)本文所提方法的影響程度，分別計(jì)算不同的時(shí)間步長(zhǎng)下的顫振臨界風(fēng)速值，計(jì)算結(jié)果如表6所示。由表6可知：對(duì)非線性Φsel處理后，時(shí)間步長(zhǎng)的選取對(duì)顫振臨界風(fēng)速的計(jì)算結(jié)果依賴性不明顯。

3.2 某大跨斜拉橋顫振臨界風(fēng)速計(jì)算

以某大跨度斜拉橋顫振分析為例，設(shè)計(jì)采用三塔斜拉橋方案，橋梁跨徑具體布置為110 m+129 m+258 m+258 m+129 m+110 m，全長(zhǎng)994 m，橋梁寬度為38.8 m，箱梁高為4.5 m，主梁標(biāo)準(zhǔn)斷面如圖5所示。

圖5 橋梁標(biāo)準(zhǔn)斷面(mm)Fig.5 Standard section of a bridge (mm)

圖6為該主梁斷面節(jié)段模型0°風(fēng)攻角下的(幾何縮尺比為1∶50)顫振導(dǎo)數(shù)試驗(yàn)及擬合結(jié)果(擬合值為無(wú)符號(hào)的曲線)，表7為氣動(dòng)力系數(shù)計(jì)算結(jié)果。

表7 橋梁斷面平板氣動(dòng)力系數(shù)Tab.7 Aerodynamic coefficients of bridge section

圖6 0°風(fēng)攻角下的顫振導(dǎo)數(shù)和擬合結(jié)果Fig.6 Flutter derivatives and fitting results under 0° wind attack angle

(a) 豎向位移時(shí)程

(b) 扭轉(zhuǎn)位移時(shí)程

0°風(fēng)攻角下，試驗(yàn)顫振臨界風(fēng)速>17 m/s，換算實(shí)橋顫振臨界風(fēng)速>110.5 m/s，通過(guò)本文方法計(jì)算得到顫振臨界風(fēng)速Ucr=198.6 m/s，顫振臨界頻率fcr=0.930 8 Hz，遠(yuǎn)高于成橋狀態(tài)顫振檢驗(yàn)風(fēng)速79.88 m/s，說(shuō)明該橋的顫振穩(wěn)定性良好。圖7為顫振臨界風(fēng)速條件下橋梁主梁斷面豎向、扭轉(zhuǎn)位移響應(yīng)時(shí)程曲線，可以直觀地判斷結(jié)構(gòu)正處于臨界運(yùn)動(dòng)狀態(tài)。

4 結(jié) 論

(1) 基于自激力有理函數(shù)表達(dá)式，推導(dǎo)了時(shí)域自激力遞推計(jì)算過(guò)程，提出了一種基于ANSYS的橋梁顫振時(shí)域分析模型和求解方法，并采用MATLAB和ANSYS混合編程技術(shù)實(shí)現(xiàn)了顫振臨界風(fēng)速的程序求解。

(2) 在相同約束條件的限定下，相比其他幾種優(yōu)化算法，擬牛頓+通用全局優(yōu)化算法在氣動(dòng)力系數(shù)識(shí)別過(guò)程中迭代步數(shù)更少，收斂更快，擬合效果更好。

(3) 以經(jīng)典平板顫振為例，采用本文所提方法計(jì)算得到的顫振臨界風(fēng)速和顫振頻率與理論吻合度很高，并且本文所提方法對(duì)計(jì)算時(shí)間步長(zhǎng)的選取不具備依賴性，某大跨橋梁顫振臨界風(fēng)速計(jì)算說(shuō)明所提方法在工程應(yīng)用中具備實(shí)用性和可操作性。

(4) 伴隨著橋梁抗風(fēng)精細(xì)化研究的進(jìn)程，結(jié)合通用有限元程序在處理非線性問(wèn)題上優(yōu)勢(shì)，本文所提方法為復(fù)雜的非線性顫振問(wèn)題研究分析提供了參考。