卿前志

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

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

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

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

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

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

(1)

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

Lse=

(2)

Mse=

(3)

當(dāng)采用脈沖響應(yīng)函數(shù)對(duì)斷面自激力進(jìn)行表達(dá)時(shí)，同樣，忽略斷面?zhèn)认蛘駝?dòng)的影響，將豎向運(yùn)動(dòng)和扭轉(zhuǎn)運(yùn)動(dòng)對(duì)斷面自激力的貢獻(xiàn)分開(kāi)，則斷面單位長(zhǎng)度氣動(dòng)自激力可表示為：

Lae=Lsea(t)+Lseh(t)

(4)

Mae=Msea(t)+Mseh(t)

(5)

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

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

(6)

(7)

(8)

(9)

(10)

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

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

(11)

代入式(4)中，有：

Lseh(t)=

(12)

對(duì)上式中的積分項(xiàng)進(jìn)行一次分步積分，可以得到Lseh(t)的最終表達(dá)式如下：

(13)

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

令

(14)

則

Zi(t+Δt)

(15)

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

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

(16)

(17)

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

圖1 ANSYS顫振時(shí)程分析流程

2 數(shù)值算例

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

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

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

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

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

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

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

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

3 結(jié) 論

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

[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(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中實(shí)現(xiàn)顫振時(shí)程分析的方法[J]．中國(guó)公路學(xué)報(bào)，2002, 15(4): 32-34．

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