卿前志

(上海市政工程設(shè)計研究總院(集團)有限公司， 上海 200092)

隨著橋梁結(jié)構(gòu)跨徑的日益增加，結(jié)構(gòu)的顫振穩(wěn)定性能成為設(shè)計過程中重要的控制因素之一。對結(jié)構(gòu)顫振穩(wěn)定性能研究的主要目的是確定結(jié)構(gòu)的顫振臨界狀態(tài)，揭示顫振發(fā)生的機理，從而保證結(jié)構(gòu)的抗風(fēng)穩(wěn)定性。目前評價橋梁的顫振性能主要有三種方法，即經(jīng)典理論法、風(fēng)洞試驗測試法和試驗加理論法[1, 2]。

至今為止，結(jié)構(gòu)顫振分析求解方法包含頻域和時域兩大類方法。頻域范圍內(nèi)對顫振問題進行求解的基本思路為：首先將氣動力用Scanlan自激力表達式[3]進行表達，然后通過求解結(jié)構(gòu)運動特征方程得到結(jié)構(gòu)顫振臨界風(fēng)速和顫振頻率。頻域分析屬于線性分析方法，與之相比，時域分析能方便地考慮各類非線性因素的影響，并反映結(jié)構(gòu)在顫振后的振幅演變規(guī)律，有利于基于過程性能的橋梁氣動性能研究，因而近年來發(fā)展迅速。在時域范圍對橋梁顫振穩(wěn)定性進行求解首先需要得到時域化的氣動自激力表達式，再在動力有限元中進行求解[2]。對于明顯存在流動分離的鈍體橋梁斷面而言，以顫振導(dǎo)數(shù)表示的自激力是一種時頻混合表達式，通常不能直接用于時域顫振分析，需要轉(zhuǎn)換成等效的時域表達。目前，顫振時域分析氣動力的表達方法有兩種。其一是通過脈沖響應(yīng)函數(shù)結(jié)合有理函數(shù)進行表達[4, 5]；其二是采用階躍函數(shù)結(jié)合橋梁斷面的氣動導(dǎo)數(shù)進行自激力模擬[6]。本文采用第一種形式，首先得到橋梁斷面的脈沖響應(yīng)函數(shù)所表達的氣動自激力，再對其進行動力有限元求解。

由于ANSYS軟件具有良好的二次開發(fā)功能，可以應(yīng)用于橋梁頻域顫振分析中[7]，但相關(guān)文獻對ANSYS應(yīng)用于橋梁時域顫振分析的報道較少。文獻[8]中雖然提出了一種在ANSYS中實現(xiàn)顫振時程分析的方法，但是，其必須同時對頻率和風(fēng)速進行搜索，因此實質(zhì)上仍然是一種頻域方法。本文提供了在ANSYS中實現(xiàn)顫振時程分析的另一方法，在瞬態(tài)動力學(xué)分析過程中嵌入氣動自激力的遞推算法來考慮其運動狀態(tài)歷史記憶特性，得到結(jié)構(gòu)的顫振臨界風(fēng)速。本文的方法也適合其他風(fēng)振問題，如顫抖振時域分析中的氣動自激力模擬。

1 氣動自激力脈沖響應(yīng)函數(shù)表達

結(jié)構(gòu)體系的運動方程為：

(1)

式中[M]、[C]、[K]分別為結(jié)構(gòu)質(zhì)量矩陣、阻尼矩陣和剛度矩陣，{Fa}為結(jié)構(gòu)上受到的荷載列向量，當(dāng)用于顫振分析時，其可以用氣動自激力的形式表達。忽略斷面?zhèn)认蛘駝拥挠绊?，斷面單位長度氣動自激力可采用如下形式：

Lse=

(2)

Mse=

(3)

當(dāng)采用脈沖響應(yīng)函數(shù)對斷面自激力進行表達時，同樣，忽略斷面?zhèn)认蛘駝拥挠绊懀瑢⒇Q向運動和扭轉(zhuǎn)運動對斷面自激力的貢獻分開，則斷面單位長度氣動自激力可表示為：

Lae=Lsea(t)+Lseh(t)

(4)

Mae=Msea(t)+Mseh(t)

(5)

式中Ifx(f=M，L；x=α，h)為相應(yīng)的脈沖響應(yīng)函數(shù)，其直觀意義為單位脈沖位移引起的氣動自激力。

根據(jù)脈沖響應(yīng)函數(shù)表達的氣動自激力與Scanlan氣動自激力表達式兩者頻譜特性一致，可得脈沖響應(yīng)函數(shù)與顫振導(dǎo)數(shù)之間的關(guān)系式：

(6)

(7)

(8)

(9)

(10)

式中，ALh1、ALh2、ALh3、ALhi、dLhi(dLhi≥0；i=4,…,m)都是待擬合的參數(shù)。第一項和第二項分別表示由位移項引起的氣動力和速度項引起的氣動力，第三項表示由加速度項引起的氣動力，該項通常較小而被忽略，第四項用于描述滯后于速度項的氣動力非定常部分，m的大小決定了這種近似的精度和附加方程的數(shù)量，所有這些參數(shù)可以通過試驗提取的顫振導(dǎo)數(shù)進行非線性擬合而確定[9]。

當(dāng)?shù)玫接欣砗瘮?shù)的相關(guān)系數(shù)后，對式(10)做傅里葉變換，可得到脈沖響應(yīng)函數(shù)的表達式：

(11)

代入式(4)中，有：

Lseh(t)=

(12)

對上式中的積分項進行一次分步積分，可以得到Lseh(t)的最終表達式如下：

(13)

有了式(13)的表達式后，在進行計算時，必須將卷積積分表達式轉(zhuǎn)化為遞推關(guān)系表達式。從式(13)易知，只需要解決Lseh2(t)的計算問題，該式可通過以下方法完成遞推：

令

(14)

則

Zi(t+Δt)

(15)

同理可以得出Lseα(s)、Mseα(s)、Mseh(s)相應(yīng)的表達式。

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

(16)

(17)

式中α=1/4(1+γ)2，δ=1/2+γ，γ為振幅衰減系數(shù)，即數(shù)值阻尼因子，ANSYS中默認為0.005，{un}及{un+1}為結(jié)構(gòu)前一時刻運動狀態(tài)量和結(jié)構(gòu)當(dāng)前時刻運動狀態(tài)量，Δt為計算時間步長，根據(jù)式(16)、(17)即可通過前一時刻的速度、加速度及當(dāng)前時刻的位移求得當(dāng)前時刻的速度和加速度大小，逐步提高風(fēng)速，便可得到結(jié)構(gòu)的顫振臨界風(fēng)速。圖1給出了在ANSYS有限元分析軟件中實現(xiàn)顫振時程分析的計算流程。

圖1 ANSYS顫振時程分析流程

2 數(shù)值算例

為驗證本文計算模型和求解方法的正確性，采用文獻[7]中具有理想平板斷面的簡支板梁進行ANSYS顫振時程分析驗證。模型參數(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=20000 kg/m，質(zhì)量慣矩Im=4.5×106kg·m2/m，空氣的密度ρ=1.248 kg/m3。建模時質(zhì)量矩陣采用集中質(zhì)量矩陣形式。橋面主梁采用BEAM4單元，質(zhì)量慣矩采用MASS21單元模擬，整個模型具有30個橋面梁單元和29個扭轉(zhuǎn)質(zhì)量單元。為便于和文獻[7]對比，結(jié)構(gòu)阻尼設(shè)為零。

對理論平板氣動導(dǎo)數(shù)進行有理函數(shù)參數(shù)擬合，表1列出了擬合結(jié)果。圖2為擬合曲線和理論曲線的對照。從圖2可以看出，由擬合有理函數(shù)參數(shù)反算顫振導(dǎo)數(shù)和理論顫振導(dǎo)數(shù)之間差別較小，擬合結(jié)果具有較高的精度。

表1 理想平板有理函數(shù)擬合參數(shù)

圖2 理想平板顫振導(dǎo)數(shù)擬合值與理論值的比較

圖3 跨中位移時程(U=135.5 m/s)

對有理函數(shù)參數(shù)擬合后，給結(jié)構(gòu)施加初始激勵進行動力響應(yīng)分析，逐步提高風(fēng)速，便可得到結(jié)構(gòu)的顫振臨界風(fēng)速。圖3、圖4、圖5分別為風(fēng)速135.5 m/s、136.6 m/s及138 m/s簡支梁跨中點位移響應(yīng)曲線。由圖3可知風(fēng)速為135.5 m/s時跨中點豎向及扭轉(zhuǎn)位移是衰減的。由圖4可知風(fēng)速為136.6 m/s時跨中點豎向及扭轉(zhuǎn)位移具有等幅特性。從圖5可知當(dāng)風(fēng)速為138 m/s時，結(jié)構(gòu)振動位移出現(xiàn)明顯的發(fā)散現(xiàn)象。得到結(jié)構(gòu)的顫振臨界風(fēng)速136.6 m/s，對顫振臨界狀態(tài)位移響應(yīng)時程做頻譜分析得到結(jié)構(gòu)的顫振頻率為f=0.3906 Hz，與文獻[7]中提供的顫振臨界風(fēng)速精確解U=136.3 m/s及顫振頻率f=0.3914 Hz非常接近。

圖4 跨中位移時程(U=136.6 m/s)

圖5 跨中位移時程(U=138 m/s)

3 結(jié) 論

橋梁時域顫振分析在非線性特性模擬、后顫振形態(tài)以及結(jié)構(gòu)風(fēng)振全過程再現(xiàn)等方面具有頻域方法不可代替的優(yōu)勢，因此時域顫振分析在將來的大跨橋梁抗風(fēng)研究中將會得到越來越廣泛的應(yīng)用。本文的算例表明，基于擬合的氣動自激力，在ANSYS中能方便地實現(xiàn)顫振自激力的施加，計算得到的顫振臨界風(fēng)速和相關(guān)文獻中的報道結(jié)果吻合良好。本文的研究表明，利用ANSYS平臺的二次開發(fā)功能進行橋梁氣動穩(wěn)定性能的時域分析是可行的，這一方法可為橋梁顫振分析或顫抖振理論分析提高效率，降低自主開發(fā)大型軟件過程中的算法錯誤或疏漏所帶來的風(fēng)險。最后，需要指出的是，采用本文時域方法對顫振臨界風(fēng)速求解的關(guān)鍵在于脈沖響應(yīng)函數(shù)參數(shù)的擬合，文中采用最小二乘法。若斷面顫振導(dǎo)數(shù)隨折算風(fēng)速變化非常劇烈，盡管能得到折算風(fēng)速范圍內(nèi)全局最優(yōu)參數(shù)，仍可能引起自激力的模擬失真，可嘗試其他參數(shù)擬合方法，繼而利用ANSYS數(shù)值計算平臺進行時域顫振分析。

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

[2] 項海帆. 現(xiàn)代橋梁抗風(fēng)理論與實踐[M]. 北京：人民交通出版社, 2005.

[3] Scanlan R H, Tomko J J. Airfoil and bridge deck flutter derivatives[J]. Journal of Engineering Mechanics, ASCE, 1971,97(6): 1717-1737.

[4] Li Q C, Lin Y K. New stochastic theory for bridge stability in turbulent flow[J]. Journal of Engineering Mechanics, ASCE, 1995, 121(1): 102-116.

[5] Lin Y K, Yang J N. Multimode bridge response to wind excitations [J]. Journal of Engineering Mechanics, ASCE,1983, 109(2): 586-603.

[6] Scanlan R H. Motion-related body force functions in two-dimensional low-speed flow[J]. Journal of Fluids and Structures, 2000, 14: 49-63.

[7] 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(1): 61-82.

[8] 華旭剛，陳政清，祝志文．一種在ANSYS中實現(xiàn)顫振時程分析的方法[J]．中國公路學(xué)報，2002, 15(4): 32-34．

[9] 張志田. 大跨度橋梁非線性抖振及其對抗風(fēng)穩(wěn)定性影響的研究[D]. 上海：同濟大學(xué)，2004.