有風攻角的棱柱體馳振計算方法研究

2018-09-27 12:58:20謝蘭博廖海黎

振動與沖擊 2018年17期

謝蘭博，廖海黎

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

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

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

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

當來流風和振動方向垂直時，將升力系數表達為風攻角的奇數次多項式形式，簡記為

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

(1)

那么振動方程可以寫為

(2)

式中：ξ為結構阻尼比；ω為結構自振圓頻率；ρ為空氣密度；B為特征尺寸；m為結構質量；vr為等效風速；v為來流風速。

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

(3)

(4)

令

(5)

則式(4)可以進行參數無量綱化

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

(6)

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

(7)

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

(8)

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

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

1 馳振臨界風速的計算

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

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

(9)

等效風速vr滿足

(10)

根據式(10)可得

(11)

結構振動方程可以寫成

(12)

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

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

(13)

(14)

則可以得到體軸下的升力系數為

(15)

(16)

觀察式(16)，當等式兩邊阻尼相等時是馳振發散的臨界點，即可得臨界風速為

(17)

上述推導都是基于體軸坐標系，而風軸坐標系下的推導思路和上述一致，即升力為

Fv=

(18)

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

(19)

根據式(19)，可以得到馳振臨界風速為

v0=

(20)

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

2 馳振振幅的計算

根據式(13)，可以得到當風攻角為α0時振動方程為

(21)

觀察等效風攻角和振動速度的關系，牽涉反正切計算，因此本文將升力系數表達為風攻角正切值的多項式形式，以方便公式推導。和Parkinson一致，本文也取到七次多項式，即

(22)

結合式(21)，可得升力表達式為

(23)

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

(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)

由于常數項只是產生一個不變動的位移，所以可以刪去，令

(26)

則振動方程可以寫為

(27)

假定結構馳振時振動為簡諧振動，振幅保持不變，且振動頻率和0風速時保持一致，即

Y=asinτ,Y′=acosτ

(28)

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

(29)

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

T(x)=mh2ω2η×

(30)

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

B1x+B3x2+B5x3+B7x4=0，

(31)

令

(32)

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

z3-αz-β=0

(33)

參考恒等式

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

(34)

可得

α=3mn，

β=n3+m3

(35)

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

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

(36)

其中

(37)

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

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

馳振達到穩定的振幅時一個周期內外力做功必然是0，但是外力做功是0并不一定能達到穩定的馳振狀態。只有當振幅稍微增大外力做功為負，振幅稍微減小而外力做功為正時，這個振幅才是穩定的振幅，否則這個振幅只能說理論上存在而已。這個條件可以表達為

(38)

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

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

在振動過程中系統的總能量為

(39)

能量對時間的導數為外力的功率，因此可以得到

(40)

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

(41)

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

(42)

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

dτ=

(43)

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

3 工程實例

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

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

實驗模型截面為H形，高寬比為0.845∶1，模型長度為1 m，豎向為y方向，其尺寸如圖2所示。

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

升力和阻力系數定義如下

(44)

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

圖3 升力和阻力系數Fig.3 Lift and drag coefficient

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

定義馳振力系數如下所示

Fdsin2α0+Fd

(45)

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

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

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

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

如圖5所示，風攻角設置和實際一致，實驗風攻角為0°、6°、8°和20°四個工況。試驗模型質量為6.5 kg，y方向振動頻率為2.82 Hz，扭轉頻率為5.4 Hz，y方向振動阻尼比約為0.05%，扭轉阻尼比約為2.5%。所測y方向位移曲線如圖6所示。

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

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

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

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

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

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

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

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

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

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

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

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

圖8 網格劃分Fig.8 Mesh partition

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

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

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