基于ANSYS的橋梁顫振時域分析及程序?qū)崿F(xiàn)

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

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

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

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

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

1 橋梁斷面自激力

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

(1)

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

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

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

(4)

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

(5)

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

1.2 氣動力系數(shù)識別

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

(6)

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

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

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

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

2 基于MATRIX27的自激力時域模擬

2.1 MATRIX27矩陣單元的組集

(7)

(8)

(9)

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

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

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

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

(10)

(11)

引入，

那么，

(12)

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

(13)

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

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

(14)

(15)

令

(16)

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

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

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

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

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

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

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

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

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

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

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

3 算例驗(yàn)證

3.1 理想平板斷面的簡支梁

為驗(yàn)證所提時域自激力計(jì)算的準(zhǔn)確性，以理想平板的顫振分析為例，并將計(jì)算結(jié)果和理論解對比分析。該簡支梁長度為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，簡支梁前10階的自振特性如表3所示，理想平板顫振導(dǎo)數(shù)如圖2所示。

(a) H*

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

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

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

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

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

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

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

(a) 豎向位移時程

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

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

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

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

以某大跨度斜拉橋顫振分析為例，設(shè)計(jì)采用三塔斜拉橋方案，橋梁跨徑具體布置為110 m+129 m+258 m+258 m+129 m+110 m，全長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é)果(擬合值為無符號的曲線)，表7為氣動力系數(shù)計(jì)算結(jié)果。

表7 橋梁斷面平板氣動力系數(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) 豎向位移時程

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

0°風(fēng)攻角下，試驗(yàn)顫振臨界風(fēng)速>17 m/s，換算實(shí)橋顫振臨界風(fēng)速>110.5 m/s，通過本文方法計(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，說明該橋的顫振穩(wěn)定性良好。圖7為顫振臨界風(fēng)速條件下橋梁主梁斷面豎向、扭轉(zhuǎn)位移響應(yīng)時程曲線，可以直觀地判斷結(jié)構(gòu)正處于臨界運(yùn)動狀態(tài)。

4 結(jié) 論

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

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

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

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