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

(1.湖南大學(xué) 風(fēng)工程試驗(yàn)研究中心，湖南長沙410082;2.上海市政工程設(shè)計(jì)研究總院(集團(tuán))有限公司，上海200092)

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

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

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

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

結(jié)構(gòu)體系的運(yùn)動(dòng)方程為:

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

忽略斷面?zhèn)认蛘駝?dòng)的影響，則作用在斷面單位長度上氣動(dòng)自激力可表示如下:

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

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

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

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

在橋梁風(fēng)致振動(dòng)問題研究中，可用如下更加靈活的形式進(jìn)行表達(dá)［6］:

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

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

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

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

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

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

式中:Lseα(s)，Mseα(s)分別表示扭轉(zhuǎn)運(yùn)動(dòng)引起的自激升力與升力矩時(shí)程;Lseh(s)，Mseh(s)分別表示豎向運(yùn)動(dòng)引起的自激升力與升力矩時(shí)程。

以式(10)最后一個(gè)等式右邊的第1項(xiàng)，即扭轉(zhuǎn)運(yùn)動(dòng)引起的自激升力為例，假設(shè)s時(shí)刻(s=Ut/B)的值已知，通過變換可得式(12):

式中:

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

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

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

根據(jù)式(15)、式(16)，可通過前一時(shí)刻的速度、加速度及當(dāng)前時(shí)刻的位移求得當(dāng)前時(shí)刻的速度和加速度大小，逐步提高風(fēng)速，便可得到結(jié)構(gòu)的顫振臨界風(fēng)速。

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

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

利用ANSYS軟件進(jìn)行橋梁顫振時(shí)域分析的基本步驟如下:

2 數(shù)值算例

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

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

首先對(duì)理論平板氣動(dòng)導(dǎo)數(shù)進(jìn)行階躍函數(shù)參數(shù)擬合，表1列出了階躍函數(shù)參數(shù)擬合值。

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

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

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

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

圖3、圖4分別為風(fēng)速137 m/s及139 m/s簡支梁跨中點(diǎn)位移響應(yīng)曲線。

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

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

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

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

圖5 相軌跡Fig.5 Phase diagrams

3 結(jié)論

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

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

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

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

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

［2］項(xiàng)海帆.現(xiàn)代橋梁抗風(fēng)理論與實(shí)踐［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］張志田，陳政清，李春光.橋梁氣動(dòng)自激力時(shí)域表達(dá)式的瞬態(tài)與極限特性［J］.工程力學(xué)，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中實(shí)現(xiàn)顫振時(shí)程分析的方法［J］.中國公路學(xué)報(bào)，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.