基于ANSYS平臺的橋梁時域顫振階躍函數(shù)算法

卿前志，張志田，肖 瑋，朱明坤

(1.湖南大學 風工程試驗研究中心，湖南長沙410082;2.上海市政工程設計研究總院(集團)有限公司，上海200092)

結構的顫振穩(wěn)定性能是超大跨度橋梁設計中至關重要的控制因素之一。目前橋梁顫振性能的評價方法主要有3種，即經(jīng)典理論法、直接風洞試驗測試法和試驗加理論法［1－2］。經(jīng)典理論法是以Theodorsen機翼顫振理論為基礎的，主要針對平板機翼的古典耦合顫振問題，后經(jīng) Bleich、Kl?ppel和 Thiele、Selberg和Van der Put等的努力，這一理論被應用于懸索橋顫振的近似計算;直接試驗法一般通過節(jié)段模型或全橋氣彈模型直接在風洞中測試橋梁的顫振性能;試驗加理論法則通過強迫振動或自由振動裝置測試斷面的顫振導數(shù)［3］，再將顫振導數(shù)應用到結構顫振理論中分析結構的氣動穩(wěn)定性能。

迄今為止，顫振理論分析包含頻域和時域兩大類方法。頻域法又可分為多模態(tài)法和全模態(tài)法，其通過將顫振運動方程轉化為特征值問題，求解方程的特征值確定顫振臨界風速。時域法則首先需獲得時域化的氣動自激力，再在動力有限元分析中進行求解［2］。頻域法屬于線性分析方法，與之相比，時域分析能方便地考慮各類非線性因素的影響，并能反映結構在顫振后的振幅演變特性，有利于基于過程性能的橋梁氣動性能研究而日益受到重視。對于存在流動分離的鈍體橋梁斷面，以顫振導數(shù)表示的自激力是一種時頻混合形式，因而通常不能直接用于橋梁時域顫振分析，需要將其轉化為時域化的表達式。目前，時域分析氣動自激力表達的方法主要有兩種:①通過脈沖響應函數(shù)結合有理函數(shù)進行表達［4］;②沿用機翼理論的方法，采用階躍函數(shù)結合橋梁斷面氣動導數(shù)進行表達［5－6］。采用有理函數(shù)法對斷面自激力進行表達存在的主要缺陷有兩個方面:①與斷面的準定常風荷載特性(平均風荷載特性)不相容;②其極限特性不具備明確的物理意義，只適合于結構響應均值為零的情況［7］。鑒于此，筆者采用階躍函數(shù)模擬自激力進行時域顫振分析。

自20世紀70年代以來，隨著計算機技術的飛速發(fā)展，國內(nèi)外相繼出現(xiàn)了許多大型通用有限元分析軟件。ANSYS軟件因其強大的前后處理和計算分析功能而在橋梁工程中有著廣泛的應用。大跨度橋梁的氣彈問題通常不能在一般商業(yè)軟件中進行分析，但由于ANSYS軟件具有良好的二次開發(fā)接口功能而被應用于橋梁結構的頻域顫振分析中［8］。但是，很少有利用ANSYS二次開發(fā)功能進行橋梁時域顫振分析的報道。文獻［9］雖然提出了一種在ANSYS中實現(xiàn)顫振時程分析的方法，但是，其必須同時對頻率和風速進行搜索，因此實質(zhì)上仍然是一種頻域方法。為了充分利用ANSYS強大的各類非線性分析功能進行大跨度橋梁的非線性氣動穩(wěn)定性研究，筆者提供了在ANSYS中實現(xiàn)顫振時程分析的另一方法，即在瞬態(tài)動力學分析過程中嵌入氣動自激力的遞推算法來考慮其運動狀態(tài)歷史記憶特性，得到結構的顫振臨界風速。本文的方法也適合其他風振問題如顫抖振時域分析中的氣動自激力模擬，為利用大型通用軟件進行大跨橋梁結構抗風分析提供參考。

1 氣動自激力階躍函數(shù)表達

結構體系的運動方程為:

式中:［M］，［C］，［K］分別為結構質(zhì)量矩陣、阻尼矩陣和剛度矩陣;{Fa}為結構上受到的荷載列向量。

忽略斷面?zhèn)认蛘駝拥挠绊?，則作用在斷面單位長度上氣動自激力可表示如下:

式中:ρ為空氣密度;U為來流風速;K=Bω/U為折算頻率;B為斷面寬度;(i=1～4)為顫振導數(shù)，一般通過風洞試驗的方法進行識別得到;˙h，˙α分別表示對實時間t的導數(shù)。

從式(2)、式(3)中可以看出，氣動自激力是折算頻率K及風速U的函數(shù)。

為了得到時域化自激力表達式，可采用階躍函數(shù)對斷面氣動自激力進行表達。假設斷面在某一時刻相對來流的攻角姿態(tài)突然產(chǎn)生一個階躍位移α0，那么，形成的氣動升力時程可表示如下:

式中:C'L為斷面升力系數(shù)對攻角的導數(shù)，C'L=dCL/dα(CL為升力系數(shù)，是風攻角α的函數(shù));φ(s)稱為階躍升力函數(shù)(s=Ut/B為無量綱時間變量)。

在橋梁風致振動問題研究中，可用如下更加靈活的形式進行表達［6］:

階躍函數(shù)所描述的階躍攻角變化引起的瞬態(tài)氣動力變化形式確定后，對于任意形式的小幅振動時程，所引起的氣動升力可根據(jù)線性疊加的原理表達如下:

式中:α'(σ)=dα/ds。

當采用階躍函數(shù)來表示橋梁斷面階躍位移引起的氣動力時，通常將豎向運動以及扭轉運動對氣動力的貢獻分開處理［10－11］。于是，單位長度梁體所受的氣動升力與扭矩分別有:

式(7)、式(8)即為采用階躍函數(shù)表示的橋梁斷面氣動自激力表達式，考慮了扭轉位移與豎向速度兩種運動狀態(tài)的影響，也有文獻考慮了側向運動對自激力的貢獻，即將式(7)與式(8)中積分項數(shù)擴充至3項，分別增加了 φLp與 φMp的貢獻［6］。

對比式(2)、式(3)與式(7)、式(8)橋梁斷面顫振力表達式后，根據(jù)兩種表達式頻譜特性一致特點，可通過積分變換的方法求得顫振導數(shù)與階躍函數(shù)中待定參數(shù)的關系，其詳細過程可見文獻［7，12］。

求出橋梁斷面升力以及升力矩對應各運動狀態(tài)的階躍函數(shù)表達式后，采用動力有限元該法進行三維顫振分析過程中，可采用式(7)與式(8)進行自激力求解。由于兩式中有關于運動狀態(tài)過程的時間積分，因此須化解成為遞推的計算式，否則自激力的計算將成為一個十分耗時的過程［7］。將自激升力與升力矩分解成與豎向以及扭轉運動相關的兩部分:

式中:Lseα(s)，Mseα(s)分別表示扭轉運動引起的自激升力與升力矩時程;Lseh(s)，Mseh(s)分別表示豎向運動引起的自激升力與升力矩時程。

以式(10)最后一個等式右邊的第1項，即扭轉運動引起的自激升力為例，假設s時刻(s=Ut/B)的值已知，通過變換可得式(12):

式中:

從式(12)可知，對于任意時刻的自激力，其求解式中的第1項只與當前運動狀態(tài)(扭轉位移或者豎向速度)有關，而對于后面若干項積分式，則與從結構運動開始到當前的整個運動時間歷程有關。因此，要根據(jù) s時刻的 Lseα(s)遞推出 s+Δs時刻的Lseα(s+Δs)，只需要解決式(13)的遞推計算即可，該式可通過以下方法完成遞推:

同理，可以得出 Lseh(s)，Msea(s)，Mseh(s)的遞推表達式。

得到斷面的氣動自激力表達式后，即可方便地在ANSYS中實現(xiàn)顫振時程分析。值得一提的是ANSYS中節(jié)點速度和加速度的求取，參照ANSYS幫助文件［13］，對于采用Newmark法的時程分析，節(jié)點加速度和速度可表示為如下差分格式:

根據(jù)式(15)、式(16)，可通過前一時刻的速度、加速度及當前時刻的位移求得當前時刻的速度和加速度大小，逐步提高風速，便可得到結構的顫振臨界風速。

圖1給出了在ANSYS有限元分析軟件中實現(xiàn)顫振時程分析的計算流程。

圖1 ANSYS顫振時程分析流程Fig.1 Flow chart of time history flutter analysis

利用ANSYS軟件進行橋梁顫振時域分析的基本步驟如下:

2 數(shù)值算例

2.1具有理想平板斷面的簡支梁

為驗證本文模型和求解方法的正確性，首先采用文獻［8］中具有理想平板斷面的簡支板梁進行驗證。模型參數(shù)為:簡支板梁長度L=300 m，寬40 m，約束兩端扭轉自由度。平板斷面的豎向和橫向剛度分別為 EIz=2.1× 106MPa·m4，EIy=1.8×107MPa·m4，扭轉剛度GIx=4.1× 105MPa·m4。每延米的質(zhì)量m=20 000 kg/m，質(zhì)量慣矩 Im=4.5 × 106kg·m2/m，空氣密度ρ=1.248 kg/m3。建模時質(zhì)量矩陣采用集中質(zhì)量矩陣形式。主梁選用BEAM4單元模擬，質(zhì)量慣矩采用MASS 21單元模擬，整個模型具有30個橋面梁單元和29個扭轉質(zhì)量單元。為便于和文獻［8］計算值進行對比，結構阻尼設為0。

首先對理論平板氣動導數(shù)進行階躍函數(shù)參數(shù)擬合，表1列出了階躍函數(shù)參數(shù)擬合值。

表1 平板氣動導數(shù)階躍函數(shù)參數(shù)擬合值Table 1 Thin foil fitted parameters of indicial functions

圖2為通過階躍函數(shù)擬合參數(shù)反算得到的顫振導數(shù)值和理論顫振導數(shù)值的對比。

圖2 理想平板顫振導數(shù)擬合值與理論值的比較Fig.2 Comparison between fitted flutter derivatives and theoretic results of thin foil

從圖2可以看出，反算的顫振導數(shù)值和理論值之間差別很小，計算結果具有較好的精度。得到階躍函數(shù)參數(shù)后，給結構施加初始激勵做動力響應分析，逐步提高風速，便可得到結構的顫振臨界風速。

圖3、圖4分別為風速137 m/s及139 m/s簡支梁跨中點位移響應曲線。

圖3 跨中位移時程(U=137 m/s)Fig.3 Time history displacements of mid-span section(U=137 m/s)

圖4 跨中位移時程(U=139 m/s)Fig.4 Time history displacements of mid-span section(U=139 m/s)

由圖3可知風速為137 m/s時跨中點豎向及扭轉運動具有等幅特性。從圖4可知當風速為139 m/s時，結構振動位移出現(xiàn)明顯的發(fā)散現(xiàn)象。得到結構的顫振臨界風速137 m/s，對顫振臨界狀態(tài)位移響應時程做頻譜分析得到結構的顫振頻率為f=0.386 4 Hz，與文獻［8］中提供的顫振臨界風速精確解U=136.3 m/s及顫振頻率f=0.391 4 Hz非常接近。

圖5為給出了風速137 m/s及139 m/s時相軌跡。對比圖5(a)、(b)的相軌跡，可知風速達到139 m/s時表征結構運動狀態(tài)的極限環(huán)存在明顯發(fā)散特性。

圖5 相軌跡Fig.5 Phase diagrams

3 結論

橋梁時域顫振分析在非線性特性模擬、后顫振形態(tài)以及結構風振全過程再現(xiàn)等方面具有頻域方法不可代替的優(yōu)勢，因此時域顫振分析在將來的大跨橋梁抗風研究中將會得到越來越廣泛的應用。

筆者討論了階躍函數(shù)對結構斷面自激力進行擬合的方法，在此基礎上，利用有限元分析軟件ANSYS的二次開發(fā)語言(APDL)對結構進行時域顫振瞬態(tài)動力有限元分析，實現(xiàn)了具有記憶特性的階躍函數(shù)自激力遞推算法與高效、可靠的大型商用軟件的結合。

算例表明，基于擬合的氣動自激力，在ANSYS三維有限元分析軟件中能方便地實現(xiàn)顫振自激力的施加，計算得到的顫振臨界風速和相關文獻中的報道結果吻合良好。

研究表明，利用ANSYS平臺的二次開發(fā)功能進行橋梁氣動穩(wěn)定性能的時域分析是可行的，這一方法可為橋梁顫振分析或顫抖振理論分析提高效率，降低自主開發(fā)大型軟件過程中的算法錯誤或疏漏所帶來的風險。

［1］陳政清.橋梁風工程［M］.北京:人民交通出版社，2005.

［2］項海帆.現(xiàn)代橋梁抗風理論與實踐［M］.北京:人民交通出版社，2005.

［3］Scanlan R H，Tomko J J.Airfoil and bridge deck flutter derivatives［J］.Journal of Engineering Mechanics，ASCE，1971(97):1717－1737.

［4］Scanlan R H.Motion－related body force functions in two－dimensional low－speed flow［J］.Journal of Fluids and Structures，2000，14:49－63.

［5］Theodorsen T.NACA Report 496:General Theory of Aerodynamic Instability and the Mechanism of Flutter［R］.VA，USA:Advisory Committee for Aeronautics，1934.

［6］Caracoglia L，Jones N P.Time domain vs.frequency domain characterization of aeroelastic forces for bridge deck sections［J］.Journal of Wind Engineering and Industrial Aerodynamics，2003，91:371－402.

［7］張志田，陳政清，李春光.橋梁氣動自激力時域表達式的瞬態(tài)與極限特性［J］.工程力學，2011，28(2):75－85.Zhang Zhitian，Chen Zhengqing，Li Chunguang.Limiting and transient characteristics of time－domain expressions for bridge selfexcited aerodynamic force［J］.Engineering Mechanics，2011，28(2):75－85.

［8］Hua X G，Chen Z Q，Ni Y Q，et al.Flutter analysis of long－span bridges using ANSYS［J］.Wind and Structures，2007，10:61 －82.

［9］華旭剛，陳政清，祝志文.一種在ANSYS中實現(xiàn)顫振時程分析的方法［J］.中國公路學報，2008，15(4):32 －34.Hua Xugang，Chen Zhengqing，Zhu Zhiwen.An approach of timehistory analysis of flutter in ANSYS［J］.China Journal of Highway and Transport，2008，15(4):32 －34.

［10］Costa C，Borri C.Application of indicial functions in bridge deck aero elasticity［J］.Journal of Wind Engineering and Industrial Aerodynamics，2006，94(9):859 －881.

［11］Salvatori L，Borri C.Frequency and time－domain methods for the numerical modeling of full－bridge aero elasticity［J］.Computers and Structures，2007，85:675 －687.

［12］Zhang Z T，Chen Z Q，Cai Y Y，et al.Indicial functions for bridge aero－elastic forces and time－domain flutter analysis［J］.Journal of Bridge Engineering，10.1061/(ASCE)BE.1943 －5592.0000176(Sep.6.2010).

［13］ANSYS Inc..ANSYS Theory Reference Documentation for ANSYS 11.0［M］.PA:ANSYS Inc.，2007.