謝蘭博， 廖海黎

(1.中鐵大橋勘測設(shè)計院集團(tuán)有限公司，武漢 430050； 2.西南交通大學(xué) 風(fēng)工程試驗研究中心，成都 610036)

橫風(fēng)向馳振是一種發(fā)生在細(xì)長結(jié)構(gòu)上較常見的氣動發(fā)散性振動。原因是在一定風(fēng)速下氣動力對結(jié)構(gòu)的氣動阻尼為負(fù)并大于結(jié)構(gòu)阻尼，對結(jié)構(gòu)做正功，使得結(jié)構(gòu)振動振幅逐漸變大，而隨著振幅的逐漸增大，氣動負(fù)阻尼降低，有可能和結(jié)構(gòu)阻尼達(dá)到平衡，使得結(jié)構(gòu)振動振幅達(dá)到穩(wěn)定。馳振的振幅較大，往往達(dá)到數(shù)倍以至于十幾倍橫風(fēng)向尺寸。而且振動頻率往往較小，一般遠(yuǎn)小于相同截面的漩渦脫落頻率[1]。

氣動發(fā)散振動的振幅計算一般認(rèn)為開始于Parkinson，他將升力系數(shù)表達(dá)為關(guān)于風(fēng)攻角的多項式形式，代入振動方程，利用漸進(jìn)法求出了方程的近似解析解，發(fā)現(xiàn)和試驗吻合的非常好[2-3]。隨后，Novak等[4-6]探討了紊流對馳振的影響，Lanevile[7]對這方面有非常詳細(xì)的闡述。Blevins[8]分析了多自由度的情況。關(guān)于馳振的振型展開，Novak[9]和Sullivan[10]都做過詳細(xì)的工作，而且又以Sullivan最具有代表性。

馳振相關(guān)的計算理論基本上都來源于Parkinson，他的基本思路如下。

當(dāng)來流風(fēng)和振動方向垂直時，將升力系數(shù)表達(dá)為風(fēng)攻角的奇數(shù)次多項式形式，簡記為

Fv=f(α)=c1α+c3α3+c5α5+…

(1)

那么振動方程可以寫為

(2)

式中：ξ為結(jié)構(gòu)阻尼比；ω為結(jié)構(gòu)自振圓頻率；ρ為空氣密度；B為特征尺寸；m為結(jié)構(gòu)質(zhì)量；vr為等效風(fēng)速；v為來流風(fēng)速。

令μ=0.5(ρv2BL)/m，則振動方程為

(3)

(4)

令

(5)

則式(4)可以進(jìn)行參數(shù)無量綱化

Y″+Y=ηf(Y′)=

(6)

上式等號右邊線性項為0時，即是馳振振動發(fā)散的臨界狀態(tài)，經(jīng)過化簡即可得到Den-Hartog公式。當(dāng)線性項大于0時，馳振振動發(fā)散。一般情況下ηA1遠(yuǎn)小于1，式(6)滿足弱非線性條件，可以求得近似解析解。根據(jù)非線性振動理論，有如下等式近似成立

(7)

式中：a是馳振振幅。實際工程中，我們更關(guān)注的是穩(wěn)定后的振幅大小，令等式兩邊為0，可得

(8)

求解式(8)即可得到穩(wěn)定的振動振幅，事實上求得風(fēng)速位移曲線也就得到了臨界風(fēng)速。

馳振的理論計算主要包括預(yù)測臨界風(fēng)速和計算馳振振幅，關(guān)于前者現(xiàn)在通用的依然是Den-Hartog公式，但是Den-Hartog公式嚴(yán)格來說僅僅適合0°攻角，即振動方向和來流風(fēng)垂直，而對于非零攻角只是一種近似。國內(nèi)大都直接利用Den-Hartog公式計算振動方向和來流方向不垂直的馳振臨界風(fēng)速，關(guān)于兩者的差別，本文后續(xù)將做詳細(xì)的討論。關(guān)于馳振的振幅計算，國內(nèi)外研究也主要集中在0°攻角的情況，對于非零攻角情況大都是一些試驗研究[11]，還未有理論計算出現(xiàn)。本文首先推導(dǎo)了振動方向和來流風(fēng)不垂直的馳振計算方法，然后利用H形截面模型開展實驗，發(fā)現(xiàn)理論計算和試驗結(jié)果吻合的比較好，證明了本文計算方法的可靠性。

1 馳振臨界風(fēng)速的計算

如圖1所示，結(jié)構(gòu)有x和y兩個振動主軸，馳振一般放生在風(fēng)攻角為0°附近或者90°附近，當(dāng)風(fēng)攻角較小，在0°附近時，結(jié)構(gòu)馳振振動為y方向；當(dāng)風(fēng)攻角較大，在90°附近時，結(jié)構(gòu)馳振振動為x方向，且振動方向都和模型主軸垂直，這兩種情況研究方法一致，因此本文選用風(fēng)攻角在0°附近的情況，即振動為y方向。結(jié)構(gòu)的風(fēng)攻角為α0，振動方向來流風(fēng)速為v，振動方向向上為正，位移為y，則根據(jù)幾何關(guān)系可得等效風(fēng)攻角α滿足

圖1 有風(fēng)攻角的馳振計算Fig.1 Calculation of galloping vibration with wind attack angle

(9)

等效風(fēng)速vr滿足

(10)

根據(jù)式(10)可得

(11)

結(jié)構(gòu)振動方程可以寫成

(12)

式中：Cv是體軸下升力系數(shù)；ξ是阻尼比；ω是振動圓頻率；m為模型質(zhì)量；ρ是空氣密度，；本文中取1.225 kg/m3；B是特征尺寸；L為模型長度。

將式(9)和式(11)代入式(12)，可得

(13)

(14)

則可以得到體軸下的升力系數(shù)為

(15)

(16)

觀察式(16)，當(dāng)?shù)仁絻蛇呑枘嵯嗟葧r是馳振發(fā)散的臨界點，即可得臨界風(fēng)速為

(17)

上述推導(dǎo)都是基于體軸坐標(biāo)系，而風(fēng)軸坐標(biāo)系下的推導(dǎo)思路和上述一致，即升力為

Fv=

(18)

將其按照一次多項式展開，并去掉常數(shù)項，可得

(19)

根據(jù)式(19)，可以得到馳振臨界風(fēng)速為

v0=

(20)

觀察式(20)，當(dāng)風(fēng)攻角為0時，就可以退化為Den-Hartog公式。從以上推導(dǎo)可以看出，Den-Hartog公式嚴(yán)格說僅適用于振動方向和來流風(fēng)垂直的情況，對于0°攻角是一種近似。

2 馳振振幅的計算

根據(jù)式(13)，可以得到當(dāng)風(fēng)攻角為α0時振動方程為

(21)

觀察等效風(fēng)攻角和振動速度的關(guān)系，牽涉反正切計算，因此本文將升力系數(shù)表達(dá)為風(fēng)攻角正切值的多項式形式，以方便公式推導(dǎo)。和Parkinson一致，本文也取到七次多項式，即

(22)

結(jié)合式(21)，可得升力表達(dá)式為

(23)

聯(lián)立式(22)和式(23)，將升力系數(shù)展開，去掉偶次冪，保留到7次冪，可得升力表達(dá)式為

(24)

其中：

c7A0=c6s(c1-c3+c5-c7)+c4s(c3-2c5+3c7)+

c2s(c5-3c7)+sc7，

c7A1=c8(2c1-2c3+2c5-2c7)+c6(-3c1+7c3-

11c5+15c7)+c4(-5c3+16c5-33c7)+

c2(29c7-7c5)-9c7，

c7A3=c6(-c1+9c3-25c5+49c7)+c4(-10c3+

60c5-182c7)+c2(217c7-35c5)-84c7，

c7A5=c4(-c3+20c5-105c7)+c2(231c7-

21c5)-126c7，

c7A7=c2(35c7-c5)-36c7，

c=cosα0,s=sinα0

(25)

由于常數(shù)項只是產(chǎn)生一個不變動的位移，所以可以刪去，令

(26)

則振動方程可以寫為

(27)

假定結(jié)構(gòu)馳振時振動為簡諧振動，振幅保持不變，且振動頻率和0風(fēng)速時保持一致，即

Y=asinτ,Y′=acosτ

(28)

則可以得到當(dāng)馳振發(fā)生時，在振動一個周期阻尼力和氣動力做功為

(29)

令振幅的平方為x，即：a2=x，則可得

T(x)=mh2ω2η×

(30)

馳振達(dá)到穩(wěn)定的振幅時一個周期內(nèi)外力做功是0，所以令式(30)兩端為0，就可以求出振動穩(wěn)定時的振幅。整理可得振幅滿足如下方程

B1x+B3x2+B5x3+B7x4=0，

(31)

令

(32)

則消去方程的二次項，式(31)可以轉(zhuǎn)化為如下形式

z3-αz-β=0

(33)

參考恒等式

(m+n)3-3mn(m+n)-(n3+m3)=0

(34)

可得

α=3mn，

β=n3+m3

(35)

所以m+n即是式(33)的一個根，根據(jù)韋達(dá)定理可以求出另外兩個根。

根據(jù)以上推導(dǎo)，可得式(31)有三個根如下所示

(36)

其中

(37)

因為p和q是實數(shù)，所以x1肯定是實數(shù)。觀察x3和x2的形式，可以看出如果三次開方下是虛數(shù)的話，則x3和x2的虛數(shù)部分正好可以抵消掉，也就是p2+q3<0時，三個根全部都是實數(shù)，而且各不相等；p2+q3>0時，方程有一個實數(shù)根和一對共軛復(fù)數(shù)根；p2+q3=0時方程有一個實根和一對相等的實根，其中如果p=0，q=0，方程有三重根。

如果僅在數(shù)學(xué)上考慮式(30)，求其振幅本質(zhì)上是研究其極循環(huán)。根據(jù)近似解析解可以證明其最多有三個極循環(huán)，其實在實數(shù)域上他最多有三個，而到復(fù)數(shù)域上，他一定是有三個極循環(huán)。此處采用七次多項式是因為他最多可能有三個極循環(huán)，而只有兩個穩(wěn)定，這和試驗比較吻合；如果采用五次多項式，則在實數(shù)域上他最多有兩個極循環(huán)，其中還有一個不穩(wěn)定；如果采用三次多項式，則極循環(huán)個數(shù)就會降到一個。因為根據(jù)實測豎向馳振最多有兩個穩(wěn)定的極循環(huán)，因此采用七次多項式擬合是比較合理的。

馳振達(dá)到穩(wěn)定的振幅時一個周期內(nèi)外力做功必然是0，但是外力做功是0并不一定能達(dá)到穩(wěn)定的馳振狀態(tài)。只有當(dāng)振幅稍微增大外力做功為負(fù)，振幅稍微減小而外力做功為正時，這個振幅才是穩(wěn)定的振幅，否則這個振幅只能說理論上存在而已。這個條件可以表達(dá)為

(38)

只要能滿足式(38)，就可以認(rèn)為xi是穩(wěn)定的振幅，否則就不是穩(wěn)定振幅。

以上計算方法可以求得馳振穩(wěn)定的振幅，但是不能求解馳振振幅逐漸增大的過程，下面給出求解發(fā)散過程振動的求解方法。

在振動過程中系統(tǒng)的總能量為

(39)

能量對時間的導(dǎo)數(shù)為外力的功率，因此可以得到

(40)

式中：Fξ代表阻尼力。假定馳振振幅在一個周期內(nèi)變化不是很大，則將式(40)在一個周期上取平均值，即可得到

(41)

結(jié)合式(39)對式(41)進(jìn)行計算并化簡a2=x，則可得

(42)

觀察式(42)，當(dāng)馳振穩(wěn)定時，振幅保持不變，上式左右都為0，即可得到穩(wěn)定的振幅。式(42)是常微分方程，求解較復(fù)雜，對其分離變量可得

dτ=

(43)

求解式(43)需要先求解出式(31)的三個解，然后分解因式，即可以求出時間關(guān)于振幅的函數(shù)。

3 工程實例

在我國，大跨度拱橋大量采用H型截面剛性吊桿，如佛山東平大橋的H型吊桿長達(dá)40.8 m。相對于圓形截面的平行鋼絲或者鋼絞線這一類較柔的吊桿，H型截面吊桿的空氣動力學(xué)性能更差，更容易發(fā)生各類風(fēng)致振動。因此本文選用H型截面模型作為研究對象開展實驗。

本文實驗在西南交通大學(xué)單回流串聯(lián)雙試驗段工業(yè)風(fēng)洞(XNJD-1)第二試驗段中進(jìn)行，該試驗段斷面為2.4 m(寬)×2 m(高)的長方形，最大來流風(fēng)速為45 m/s，最小來流風(fēng)速為0.5 m/s。

實驗?zāi)Ｐ徒孛鏋镠形，高寬比為0.845∶1，模型長度為1 m，豎向為y方向，其尺寸如圖2所示。

圖2 模型截面尺寸(mm)Fig.2 Model section size (mm)

升力和阻力系數(shù)定義如下

(44)

式中：U=14.6 m/s，B=0.1 m，L=0.3 m。所測風(fēng)軸下阻力和升力系數(shù)如圖3所示。

圖3 升力和阻力系數(shù)Fig.3 Lift and drag coefficient

圖3中散點是試驗所測結(jié)果，而曲線是利用25次多項式擬合所得結(jié)果，數(shù)據(jù)相對較光滑，本文插值計算采用25次多項式擬合的結(jié)果。

定義馳振力系數(shù)如下所示

Fdsin2α0+Fd

(45)

式(45)中SDe代表的是利用Den-Hartog公式所定義的馳振力系數(shù)，Sα是根據(jù)本文式(20)所定義的馳振力系數(shù)，由于本文推導(dǎo)的臨界風(fēng)速公式振動方向是確定的，而振動結(jié)構(gòu)一般都是有兩個主軸，所以式(20)只能代表一個主軸方向，另一個方向的公式需要做一個旋轉(zhuǎn)變換。根據(jù)式(45)和本文所測風(fēng)軸的三分力系數(shù)，可得兩種定義方法的馳振力系數(shù)對比，如圖4所示。

圖4 馳振力系數(shù)對比Fig.4 Comparison of galloping force coefficient

如圖4所示，本文所推導(dǎo)結(jié)果和Den-Hartog公式結(jié)果在0°攻角和90°攻角附近非常接近，在70°～80°攻角范圍內(nèi)兩者誤差稍大，而在20°～65°攻角范圍內(nèi)兩者則有很大的誤差，由于Den-Hartog公式假定振動方向和來流風(fēng)垂直，因此，在有較大風(fēng)攻角的臨界風(fēng)速判定建議采用本文所推導(dǎo)公式。觀察圖4還可以發(fā)現(xiàn)，0°攻角附近比90°攻角附近更容易發(fā)生馳振，因此本文重點研究0°攻角附近的馳振振幅。

風(fēng)洞試驗對于風(fēng)攻角設(shè)置方面大致有兩種思路，第一種是改變模型的傾斜角度，而彈簧依然和來流風(fēng)垂直；另一種是改變模型和彈簧的傾斜角度，使得能夠?qū)崿F(xiàn)和事實符合的風(fēng)攻角。第一種情況容易實現(xiàn)，而且在風(fēng)攻角比較小的情況下，所得結(jié)果和實際相差很小，但是在風(fēng)攻角較大時，所得結(jié)果將會和實際產(chǎn)生較大誤差。第二種情況實現(xiàn)較麻煩，但是和實際一致，本文實驗采用第二種方法。

如圖5所示，風(fēng)攻角設(shè)置和實際一致，實驗風(fēng)攻角為0°、6°、8°和20°四個工況。試驗?zāi)Ｐ唾|(zhì)量為6.5 kg，y方向振動頻率為2.82 Hz，扭轉(zhuǎn)頻率為5.4 Hz，y方向振動阻尼比約為0.05%，扭轉(zhuǎn)阻尼比約為2.5%。所測y方向位移曲線如圖6所示。

圖5 風(fēng)攻角設(shè)置示意圖(mm)Fig.5 Schematic plot of wind attack angle (mm)

從圖6可知，在0°攻角工況下位移增長最快，攻角越大，增長越慢，而20°攻角情況下，馳振消失。

(a) 0°攻角風(fēng)速位移曲線

(b) 6°攻角風(fēng)速位移曲線

(c) 8°攻角風(fēng)速位移曲線

(d) 20°攻角風(fēng)速位移曲線圖6 風(fēng)速位移關(guān)系Fig.6 The relationship between wind speed and displacement

已知結(jié)構(gòu)三分力系數(shù)，利用多項式擬合進(jìn)行計算無論如何擬合數(shù)據(jù)都會和原數(shù)據(jù)有一定的差異，所以本節(jié)先介紹利用三次樣條插值進(jìn)行升力系數(shù)擬合而進(jìn)行數(shù)值計算的求解結(jié)果。圖7是模型分別在0°、6°和8°攻角下的理論計算和實驗的對比。因為20°攻角未發(fā)現(xiàn)馳振現(xiàn)象，所以在此不做討論。

(a) 0°攻角試驗和理論馳振振幅對比

(b) 6°攻角試驗和理論馳振振幅對比

(c) 8°攻角試驗和理論馳振振幅對比圖7 試驗和理論馳振振幅對比Fig.7 Experimental and theoretical comparison of the amplitude

由圖7可知，8°攻角位移曲線和試驗吻合非常好，其次是6°攻角，而0°攻角理論計算和試驗所測位移有一定的偏離。但是三者在較高的風(fēng)速時，理論計算和實驗都趨向于一致，說明在高風(fēng)速下準(zhǔn)定常理論成立的較好。這說明了本文所推導(dǎo)有風(fēng)攻角計算馳振的公式的正確性。在低風(fēng)速時，0°攻角和6°攻角計算和理論偏移較大，甚至0°攻角計算的臨界風(fēng)速和實驗值相差到了40%～50%，產(chǎn)生以上誤差的原因可能是：

(1)準(zhǔn)定常理論的局限性。準(zhǔn)定常理論本身就會有一定的誤差，氣動力總會存在滯后現(xiàn)象。

(2)渦激力的影響。如果在馳振的風(fēng)速下漩渦脫落頻率和結(jié)構(gòu)振動頻率接近，則會導(dǎo)致試驗振幅和計算的結(jié)果偏離很遠(yuǎn)。本實驗馳振臨界風(fēng)速非常低，只有1.5 m/s左右，為了弄清楚渦激力對馳振的影響，需要計算出結(jié)構(gòu)的斯托羅哈數(shù)。本文借助Fluent軟件利用CFD進(jìn)行數(shù)值模擬，求取結(jié)構(gòu)的斯托羅哈數(shù)。網(wǎng)格劃分如圖8所示。

邊界條件設(shè)置為：上游速度入口，使用Velocity-inlet邊界條件，湍流模型采用SSTk-ω，湍流強度取0.2%，湍流黏性比取10%，下游使用pressure-outlet邊界條件，出口相對壓強平均值取0，上下側(cè)邊界條件是Symmetry邊界條件，H形截面的表面采用No-Slip Wall邊界條件。

圖8 網(wǎng)格劃分Fig.8 Mesh partition

求解設(shè)置為：靜力計算壓力-速度耦合算法使用PISO算法，離散格式控制方程采用QUICK格式進(jìn)行求解，計算時間步長0.005 s。采用Interface邊界利用滑移網(wǎng)格方便旋轉(zhuǎn)模型，求得結(jié)果如表1所示。

觀察表1可以看出，對于0°攻角，St≈0.13，當(dāng)振動頻率為2.82 Hz時，對應(yīng)渦激共振風(fēng)速為2.16 m/s，而這個風(fēng)速恰好是0°攻角馳振剛發(fā)生不久，同理計算出6°攻角及8°攻角對應(yīng)渦激共振風(fēng)速分別為1.96 m/s和1.86 m/s，都是對應(yīng)馳振發(fā)散剛剛發(fā)生，所以有可能渦激力影響了理論計算和實驗結(jié)果之間的差別。

利用樣條插值計算馳振位移雖然比多項式擬合更精確，但是它不能寫出解析形式，因此其應(yīng)用有一定的局限性。圖9是利用七次多項式進(jìn)行升力系數(shù)擬合結(jié)果。

表1 斯托羅哈數(shù)計算結(jié)果Tab.1 Results of calculation of Strouhal number

利用最小二乘法，可得升力表達(dá)式為

Cv=c1tanα+c3tan3α+c5tan5α+c7tan7α

(46)

式中：c1=-7.6，c3=131.8，c5=-635.9，c7=1 016.9，代入到式(25)，求出各攻角之下各風(fēng)速之下的振幅，并用式(38)判斷是不是穩(wěn)定的振幅。利用多項式求得馳振振幅如圖10所示。

圖9 多項式擬合升力系數(shù)結(jié)果Fig.9 Polynomial fitting lift coefficient results

(a) 0°攻角馳振振幅試驗和理論對比

(b) 6°攻角馳振振幅試驗和理論對比

(c) 8°攻角馳振振幅試驗和理論對比圖10 馳振振幅試驗和理論對比Fig.10 Experimental and theoretical comparison of the amplitude

觀察圖10，可以看出多項式擬合和插值計算結(jié)果比較接近，這是因為升力多項式擬合結(jié)果和原來的數(shù)據(jù)誤差很小。

4 結(jié) 論

本文推導(dǎo)了有風(fēng)攻角的馳振臨界風(fēng)速和振幅計算公式，利用Den-Hartog公式求解臨界風(fēng)速和本文公式對比，并采用H形截面模型開展實驗對振幅加以驗證，討論了理論計算和試驗產(chǎn)生誤差的原因，有以下結(jié)論：

(1)Den-Hartog公式并不適用于來流風(fēng)和振動方向夾角較大的情況，通過和本文公式對比，認(rèn)為在較小的風(fēng)攻角之下，Den-Hartog公式誤差不大，但是較大的風(fēng)攻角，Den-Hartog公式計算結(jié)果會慢慢偏離本文所推導(dǎo)的公式計算結(jié)果，因此建議對于有風(fēng)攻角的馳振計算應(yīng)采用本文所推導(dǎo)公式。

(2)本文所推導(dǎo)的有風(fēng)攻角馳振振幅計算公式計算結(jié)果和試驗結(jié)果較為吻合，證實了本文所推導(dǎo)計算公式的正確性。