基于CFD數值計算的大跨懸索橋碳纖維矩形截面吊桿風致振動研究

張建國, 張 春

(廈門大學 建筑與土木工程學院,福建 廈門 361005)

當今,在超大跨徑橋梁的建設中,懸索橋毫無疑問是競爭力最強的一種橋型,其跨越能力是其他類型橋梁無可匹及的。吊索作為懸索橋的重要傳力構件,其性能對懸索橋的影響巨大。由于吊索具有長細比大、質量輕、阻尼小等特點,使得在風荷載作用下其發生振動的概率大大增加[1]。

華旭剛等[2]強調長吊索和特長吊索的風致振動必須重點關注。馬存明等[3]通過節段的靜力和彈簧懸掛動力試驗,對H型吊桿的馳振和渦激振動性能進行分析。陳政清等[4]對H型吊桿的大攻角風致振動進行探究并提出直立桿件抗風設計的建議。岳明[5]使用數值模擬的方法研究了懸索橋并列吊桿的渦激振動,該研究的重點為研究不同長度吊桿的渦激振動。研究發現吊桿越長,吊桿剛度越大,越容易發生渦激振動且發生渦激振動的風速較小,易造成吊桿的疲勞損傷。徐昕宇等[6]采用CFD(computational fluid dynamics)模擬,研究扁平板式吊桿風致渦激振動特性,發現開槽后的吊桿斷面在各風速下振幅接近0,吊桿開槽能夠有效抑制渦振的發生。

以上對懸索橋吊桿的研究都是基于傳統的鋼絞線或H型鋼材料制作的吊桿進行風振探討,本文提出利用新型碳纖維材料制作吊桿,探究吊桿風振特性。與傳統鋼絞線相比,碳纖維材料具有質量輕、強度高、耐腐蝕的特點,為橋梁的安全提供了有力的保障[7]。由于碳纖維材料特性,吊桿做成矩形截面的情況較多,然而矩形截面容易發生馳振現象,故而矩形截面吊桿僅在歐洲為數不多的幾座橋梁上有所應用,國內橋梁應用的實例少,本文提出在滿足力學性能要求的基礎上,對矩形碳纖維吊索采用CFD方法研究不同風攻角、不同截面比對吊桿渦激振動以及馳振可能性的影響。

1 工程概況及自振頻率

1.1 工程概況

國內某大跨度懸索橋吊桿選用碳纖維材料后的鉸接式吊桿布置如圖1所示。

圖1 吊桿編號Fig.1 Hanger number

在圖1中選擇具有代表性的5根吊桿進行研究,分別為左端部1號,1/4跨20號,跨中39號,3/4跨58號以及右端部77號,如圖1圓圈所示。橋梁設計人員進行了設置碳纖維吊桿工況下大跨度懸索橋梁的受力分析和結構設計,求得更換后的典型碳纖維索截面積均為0.021 m2,密度為1.6 g/cm3,彈性模量為1.3×1011Pa,泊松比為0.307。表1為吊桿在重力荷載代表值下計算得到的索力值。

表1 替換后碳纖維索的索力

在進行吊桿風致振動研究時,可將吊桿簡化為質量-彈簧-阻尼系統,那么要研究吊桿的渦激振動問題,最重要的就是確定吊桿的剛度K,而吊桿的剛度與吊桿的自振頻率相關,因而確定吊桿的剛度問題也就轉化為吊桿自振頻率的求解問題。吊桿自振頻率的求解有兩種方法:有限元模態分析和理論公式計算法。

1.2 有限元模態分析

在ANSYS有限元分析軟件中,采用LINK180單元用來模擬吊桿,計算選擇Block Lanczos法計算前10階頻率,所得計算結果如表2所示。

表2 吊桿自振頻率

1.3 自振頻率理論計算

為了驗證吊桿模態分析所得到的自振頻率的正確性,采用考慮預應力的吊桿自振頻率理論計算公式進行吊桿的一階自振頻率計算。

Clough等[8]提出理論:當“弦”兩端固定,不考慮抗彎剛度的影響時張力與頻率之間的關系為

(1)

式中:T為吊桿張力;m為吊桿單位長度質量;l為吊桿長度;fn為吊桿第n階橫向振動頻率。

參考表1給出的吊桿最大索力,計算轉換為最大初張力,結合已知的碳纖維吊桿的物理參量,采用上述理論方法計算得到吊桿的一階自振頻率,并與有限元模態分析計算得到的吊桿一階自振頻率表2進行對比得到表3。

表3 吊桿一階自振頻率理論法和有限元法計算比較

由表3給出的吊桿一階自振頻率理論公式求解和有限元模態分析求解結果比較可知,理論公式法的結果與有限元模態分析求解的結果基本相同,說明本文的有限元模型有效,因此,可運用本文有限元模型求得的結果進行后續風致振動研究。另外,端部1號和端部77號兩根吊桿的頻率極其接近,而1/4跨20號和3/4跨58號兩根吊桿的自振頻率也很接近,因此,在后續的風致振動分析時,沒有針對3/4跨58號吊桿和右端部77號吊桿進行研究,僅選取了1號、20號和39號3根吊桿進行詳細分析。

2 CFD計算模型

2.1 控制方程

計算風工程又稱數值風洞方法,與直接風洞試驗相比具有模擬真實風環境的能力,可以構建原型尺度的計算模型,還具有周期短、費用低的特點。在直角坐標系下,對二維不可壓縮黏性流體,其運動規律用納維-斯托克斯方程(N-S方程)[9]描述

(2)

(3)

式中:ρ為空氣密度;ν為運動黏性系數;u為速度矢量,f為作用在單位質量流體微團的質量力;p為流體壓強。

2.2 數值模擬

對數值計算而言,計算域選取的大小對于結果的影響很大。從理論上來講,計算域越大,與實際情形越接近,結果就會越精確。但由于計算機內存及時間的限制,不可能選取太大的計算域[10]。文獻[11]認為計算域的進口、上下邊界距柱體的距離至少為10D,出口距柱體的距離應大于20D,才能保證結果的獨立性。本文計算域取為42D×24D,進口處距吊桿12D,出口距吊桿30D,上下邊界距吊桿12D。此外,當流體流過鈍體的時候,由于流體本身黏性的作用,使得物體表面上的流體速度為零,而離開物體表面一定距離的流體速度則不受黏性影響,此處的流動可以按照無黏來處理。在物面和流體之間的可以按無黏處理的這一部分流體就是邊界層。邊界層通常是很薄的一層,但在這一薄層中各項流體參數卻發生著劇烈的變化,存在著較大的速度梯度。因此,在繞流問題中精確地模擬邊界層對于結果至關重要,這就要求邊界層網格要足夠密[12],本文所劃網格如圖2所示。

圖2 截面比為1 ∶1時計算域及網格示意圖Fig.2 Schematic diagram of the calculated domain and mesh with a cross-section ratio of 1 ∶1

SSTk-ω模型合并了來源于ω方程中的交叉擴散,湍流黏度考慮到了湍流剪應力的傳播,而且模型常量不同,這些特點使得該模型在近壁區域的自由流、逆壓梯度流動、翼型、跨音速激波等方面有著更高的精度[13]。本數值模擬選擇SSTk-ω湍流模型,具體數值模擬中,采用SIMPLE格式求解壓力速度耦合方程組,空間離散采用二階迎風格式,時間離散采用二階全隱格式。已有的研究數據表明,方柱繞流的無量綱頻率約為0.12,在一個渦脫落周期內至少需要20～25個時間步長[14],經計算后時間步長取Δt=0.02 s。

計算域邊界條件設置如下:

(1)左側,速度入口,U=U0,V=0。

(2)右側,出口邊界條件,平均靜壓P為零。

(3)沿展向兩邊界采用對稱邊界條件,各變量沿法向的分量為零。

(4)其余邊界采用無滑移固壁邊界條件,U=0,V=0。

2.3 模型驗證

首先進行了網格獨立性測試。對1 ∶1截面時0°攻角計算共選用了4種網格數量(時間步長均為0.02 s)進行了測試,網格獨立性測試結果如表4所示。表5為已有文獻單方柱0°攻角時繞流數據對比。

表4 網格獨立性驗證

表5 文獻單方柱0°攻角時繞流數據

根據表4和表5結果,選擇40 000網格數進行后續計算既能滿足計算精度要求,又能達到計算效率要求。

3 數值模擬結果

3.1 馳振判斷

馳振是由于負氣動阻尼引起的橫風向結構振動,常出現在具有特殊截面形狀的細長結構中。Den Hartog在研究覆冰導線時提出了橫風向馳振判據式(4)。

(4)

圖3和圖4為不同截面比且攻角為0°時的三分力系數時程曲線圖;圖5～圖7為平均升力系數、平均阻力系數數值模擬結果以及不同截面馳振系數隨攻角變化情況。

圖3 攻角為0°時升力時程曲線Fig.3 Lift force time curve of 0° attack angle

圖4 攻角為0°時阻力時程曲線Fig.4 Damp force time curve of 0° attack angle

圖5 平均升力系數隨攻角變化曲線Fig.5 The mean lift coefficient varies with the angle of attack curve

圖6 平均阻力系數隨攻角變化曲線 Fig.6 The mean damp coefficient varies with the angle of attack curve

圖7 不同截面馳振系數隨攻角變化圖 Fig.7 The galloping coefficient of different sections with the angle of attack

按準定常馳振理論[20],馳振臨界風速如式(5)所示。

(5)

式中,B為相應的迎風面寬度,同時根據陳政清等的觀點,認為維持穩定的馳振狀態的區間應在一定攻角范圍以內。按一般設計思想取最小馳振系數進行計算,雖然最為保守但實際上卻是明顯不合理的。因為在實際的紊流風作用下桿件風攻角(即風場的方向角)是變化的,攻角大于馳振度數后桿件會迅速退出馳振狀態,故而認為取處于穩定的馳振區間的平均值較為合理。表6為不同截面比矩形吊桿的臨界馳振風速計算結果。

表6 臨界馳振風速計算結果

按照Parkinson等于1964 年提出的非線性準定常模型可知,此理論中橫風向脈動氣動力系數可以表示成

(6)

CFY可以由柱體靜風試驗直接測量得到的平均升力系數和平均阻力系數推算求得

CFY(β)=-[CL(β)+CD(β)·tan(β)]·sec(β)

(7)

利用該數學模型預測所得柱體氣動失穩(馳振)臨界風速

(8)

對1 ∶1截面所得數據進行處理后,曲線擬合所得系數A=1.27,B=29.54,C=72.51,D=55.92,代入式(8)中計算可得該模型預測1 ∶1截面比時馳振臨界風速為端部吊桿臨界馳振風速為66.08 m/s,1/4跨處吊桿臨界馳振風速為413.80 m/s,跨中吊桿臨界馳振風速為1 934.18 m/s。表6中計算數據與此模型計算結果相差不大,數據可信。

根據表6結果,可知對于跨中及1/4跨處吊桿,計算出的臨界馳振風速均較大,大于懸索橋的設計風速,即不用考慮吊桿馳振。對于端部吊桿而言,當截面比為1 ∶1和1 ∶2時不用考慮馳振情況,但是當截面比為1:3時,端部吊桿相應的臨界馳振風速較低,需要考慮相應的抗馳振措施。

3.2 渦振判斷

當大氣邊界層中近地風繞過橋梁時,可能會產生流動分離和周期性的旋渦脫落,使結構上下或左右兩側表面出現交替變化的正負壓力和力矩,稱為渦激力。渦激力可能會引起結構橫風向或扭轉方向有限振幅的振動,稱為渦激振動,簡稱渦振。渦振雖然不會像馳振、顫振一樣導致橋梁產生災難性的破壞,但是它發生時的風速低、頻率高,有可能導致桿件裂紋或疲勞破壞,或影響行車的舒適性和安全性[21]。

對FLUENT計算得到的1 ∶1吊桿的力矩系數Cm時程數據進行FFT(fast Fourier transformation)后可以得到其頻譜圖如圖8所示。可見力矩系數頻譜存在明顯的峰值頻率。

圖8 截面為1 ∶1 0°攻角時力矩頻譜圖Fig.8 FFT of the moment coefficient at 0°attack angle and 1 ∶1 cross section

匯總所有風向角下3種不同截面吊桿力矩系數的峰值頻率數據,給出斯特勞哈爾數St如式(9)所示。

St=fSD/U

(9)

式中:fs為旋渦脫落頻率;U為本次數值模擬中的平均風速,3 m/s。圖9即為主頻對應的斯特勞哈爾數隨風攻角變化圖。

圖9 主頻對應的斯特勞哈爾數隨風攻角變化圖Fig.9 The main frequency corresponding to the Strouhal number changes with the attack angle

對于吊桿,如果旋渦脫落的卓越頻率fs與在橫向風作用下吊桿振動特征頻率f相同,就會導致渦激共振。因此,由fs=f定義的臨界風速Ucrit為式(10)。

(10)

Sohankar[22]的研究發現無論邊界層的性質如何,尖銳邊緣的物體都會導致流動分離,但尤其是在較高的雷諾數下對雷諾數不敏感。由于本文研究的吊桿為矩形截面,且研究為高雷諾數情況,因此可不考慮雷諾數的影響。而Re=UD/ν,即在此次研究中不考慮空氣運動黏性系數以及特征長度的改變,也就是說Re僅與風速有關,即本次研究中風速的改變基本上對St無影響。由此知D/St越小,同時自振頻率越小,則渦振風速越小,故而對不同截面比的吊桿取相應的D/St最小值,用對應吊桿自振頻率進行計算,可以推算其最小渦振風速見表7。由表7可知,端部吊桿在3種截面比情況下,均容易發生多階的渦激振動,最小渦振風速都較小,而1/4跨和跨中吊桿,在3種截面比時,僅可能發生低階的渦振,高階渦振對應的風速均較大。

表7 渦振風速

3.3 流場分析

基于FLUENT后處理軟件繪制本次數值模擬中二維矩形斷面的風速云圖及風壓云圖,本文通過這些云圖來分析吊桿周圍的流場分布情況,進一步說明矩形吊桿的渦激共振的機理,具體云圖如圖10所示。需要說明的是,由于本文沒有進行流固耦合分析,流場分析中未出現“頻率鎖定”現象。

圖10 風速云圖(左)與風壓云圖(右)Fig.10 Wind speed cloud graph (left) and wind pressure cloud graph (right)

從云圖中可以看出,氣流在吊桿的上下側出現明顯的來流分離和旋渦脫落,并在吊桿后方產生空腔。吊桿上下側旋渦交替脫落以及吊桿后方的空腔是引起吊桿渦振的重要因素。吊桿后方的尾流區上下兩側旋渦交替反向,總有一側的尾渦占主導地位。隨著截面比變小,流體附著點越多,出現再附著的可能性更大。隨著攻角變化,分離點發生變化,同時附著點也發生變化,對于非正方形截面,當吊桿處于大風攻角時,尾流相較而言紊亂得多。對于吊桿,迎風面所受風壓大,在尾部形成負壓區,同時尾流區散布負壓區。

4 結 論

文章對吊桿采用碳纖維材料后對其采用不同截面比的矩形斷面進行了相關風振情況研究,通過數值模擬方法,探究其產生渦激共振與馳振的可能性以及相應的臨界風速計算,得出以下結論:

(1)1 ∶1方形截面和1 ∶2矩形截面吊桿在抗馳振方面表現良好。

(2)無論截面比為多少,端部的吊桿比跨中和1/4跨處吊桿更容易發生渦振。

(3)截面比越小,其渦振風速越小,且更容易發生多階渦振。