CFD方法的大跨度橋梁抖振荷載

劉婷婷，張文首

（大連理工大學 工業(yè)裝備結構分析國家重點實驗室，遼寧 大連116024）

風荷載是大跨度橋梁的主要環(huán)境荷載，其中，抖振荷載作為一種長期作用的隨機荷載，會使橋梁結構局部構件和連接部位處于交變應力作用，影響橋梁疲勞壽命，是大跨度橋梁抗風設計的重點問題。

目前的抖振荷載計算主要依靠風洞試驗方法。在一定縮比的橋梁節(jié)段模型上進行定常流風洞試驗，測定抖振力三分力氣動系數(shù)，然后帶入有限元模型（通常為魚骨頭模型）進行分析［1－5］。但是，風洞試驗方法的周期長、費用高、結果的可視性較差，并且，無法針對橋梁的精細有限元模型進行直接分析，制約了大跨度橋梁抖振分析方法的發(fā)展。計算流體力學（Computational Fluid Dynamic，CFD）的發(fā)展使數(shù)值模擬方法可以替代物理的風洞試驗，能夠節(jié)約試驗成本、進行重復試驗，彌補風洞試驗所測數(shù)據(jù)不足［6－8］。

采用CFD數(shù)值模擬替代傳統(tǒng)的風洞試驗，改進了大跨度橋梁抖振荷載計算。在香港青馬大橋精細有限元模型上，模擬橋板表面壓力分布，根據(jù)CFD模擬結果確定節(jié)點抖振力，由此求解橋板局部應力響應時程。CFD模擬能夠獲得橋板表面任何位置壓力，彌補了風洞試驗安裝測壓傳感器有限的缺點。而且，可以針對橋梁精細有限元模型，直接確定節(jié)點抖振力，將風致抖振分析范圍推進到局部構件，為橋梁風振疲勞研究提供了有力的工具，促進了橋梁健康監(jiān)測及局部抗風設計的發(fā)展。

1 抖振力

1．1 傳統(tǒng)抖振力

根據(jù)擬定常假定和片條原理，長度為L橋板節(jié)段彈性中心的阻力Dbfe，升力Lbfe及彎矩Mbfe可以表達為（如圖1）［9］：

式中，ρ是空氣密度，U是平均風速，B為橋板寬度，u（t）和w（t）為水平及豎直方向脈動風速時程。CD、CL和C M是與攻角α相關的無量綱三分力系數(shù)，主要通過風洞試驗在橋板彈性中心測得，C′D＝dCD／dα、C′L＝ dCL／dα、C′M＝ dCM／dα。χDbu，χDbw，χMbu，χMbw，χLbu，χLbw是氣動導納函數(shù)，當長度L較小時可看作常數(shù)1。

圖1 橋板彈性中心及節(jié)點抖振力

1．2 橋板壓力分布的CFD模擬

實質上，采用以上傳統(tǒng)方法計算抖振力時，就是將作用在橋板上的壓力合成為彈性中心的合力，因此，在風洞試驗中只關注彈性中心的氣動參數(shù)，而忽略了橋板的壓力分布及其產(chǎn)生的局部作用。另外，隨著有限元建模技術的發(fā)展，已經(jīng)可以建立更精細的橋梁有限元模型，以適應復雜分析的需要。例如，為了研究交通荷載引起的大跨度懸索橋疲勞問題，文獻10建立了青馬大橋的大型有限元模型，采用梁單元模擬橋板的桁架構型，并建立橋板節(jié)點的局部模型；針對青馬大橋健康監(jiān)測系統(tǒng)——WASHMS（Wind and Structure Health Monitoring System），更準確的青馬大橋有限元模型也被實現(xiàn)，并應用到風振引起的橋梁疲勞問題的預測［11－12］。與魚骨頭模型不同，采用精細有限元模型進行抖振計算時，需要的是單元節(jié)點的集中抖振力，這是單一采用風洞試驗方法無法獲得的。

橋板風壓分布是計算單元節(jié)點抖振力關鍵，常用方法是制作布置了壓力傳感器的橋板節(jié)段模型，通過風洞試驗獲得。CFD技術的發(fā)展使數(shù)值模擬方法可以替代物理的風洞試驗，并能夠節(jié)約試驗成本、彌補風洞試驗測量數(shù)據(jù)的不足。計算流體力學CFD的思想是：把時間及空間上連續(xù)的物理量（速度或壓力），用有限個離散點上的變量集合代替，通過流動基本方程建立關于這些離散點上變量關系的代數(shù)方程組，然后進行求解獲得近似值。通過這種對流體流動的數(shù)值模擬，可以得到復雜問題的流場內各個位置上的基本物理量（速度、壓力、溫度等）的分布，甚至這些物理量隨時間變化的情況。

對于橋板壓力分布，CFD數(shù)值模擬的對象是橋板周圍流場，可以采用二維定常不可壓縮控制方程［13－14］，考慮湍流效應，引入標準k－ε雙方程湍流模型。

其中，u和v分別是流場內沿水平軸和豎直軸的速度，ρ是空氣密度，p為壓力，k和ε分別是是湍流動能和湍流動能耗散率；有效粘性系數(shù)υeff＝υ＋υ1，其中，υ為運動粘性系數(shù)，是湍流粘性系數(shù)；系數(shù)為：Cμ＝0．09、C1＝1．44、C2＝1．92、σk＝1．0以及σε＝1．3。

1．3 節(jié)點抖振力

若通過CFD數(shù)值模擬已知橋板壓力分布，可將壓力轉換為作用在橋板表面梁單元的節(jié)點力：

其中，?Fl1和?Fl2是第l個單元兩端節(jié)點力，nl是CFD

式中，p k1和p k2為計算網(wǎng)格節(jié)點處壓力，l k是網(wǎng)格長度。

將式3計算的單元節(jié)點力由單元坐標轉換為p－h(huán)－α坐標后，疊加各相鄰單元節(jié)點力，就可將橋板壓力等效為各單元節(jié)點力，其中第j個節(jié)點力為：是作用的節(jié)點力總數(shù)。則節(jié)點氣動參數(shù)可以表示為：

式中，是CFD模擬中的平均風速。

最后，式1可以改寫為節(jié)點抖振力形式：

2 精細有限元模型的結構運動方程

在全局坐標系x－y－z中，在抖振力作用下，橋梁運動方程為：

其中，a（t）＝ ｛u1，w1，…，um，wm｝T是2m維脈動風速向量，m為橋板節(jié)段總數(shù)；Pbf為6N×2m氣動參數(shù)矩陣，由每個橋板節(jié)段的氣動參數(shù)矩陣Q ibf組成。

大跨度橋梁的有限元模型的自由度較大，計算量十分龐大，通常采用模態(tài)疊加法進行分析。節(jié)點位移X（t）可表示為：模擬計算中將此單元劃分的網(wǎng)格數(shù)量；和是在第k個網(wǎng)格末端的壓力合力（如圖1），其表達為：

式中，Φ＝ ［Φ1，Φ2，…，ΦNΦ］是6N×NΦ模態(tài)矩陣，q（t）＝ ｛q1（t），q2（t），…，q NΦ（t）｝T是廣義位移向量，NΦ是計算中采用的模態(tài)總數(shù)。式7的模態(tài)運動方程為：

模態(tài)運動方程（10）可采用直接積分法求解廣義位移時程q（t），然后由式9確定節(jié)點位移。如果已知某單元模態(tài)應力Γr，此單元應力為：

3 算例研究

3．1 青馬大橋精細有限元模型

香港青馬大橋是連接市區(qū)和大嶼山機場的主要通路，是世界目前最長的公鐵兩用懸索橋之一，主跨達1 377 m。青馬大橋三維精細有限元模型采用商業(yè)軟件MSC NASTRAN建立（圖2）。有限元模型的幾何構型與真實橋梁一致，采用三維2節(jié)點梁單元模擬加勁桁架，主索及吊桿用圓截面梁單元模擬，橋面采用正交各向異性板單元，用等效剛度方法考慮加勁肋對橋面剛度影響。2個橋塔是多室空間結構，采用空間梁單元進行劃分，并對其剛度進行換算。模型邊界條件完全符合實際情況。對該模型進行了動力特性分析，與實測結果符合很好［6］。

圖2 青馬大橋三維精細有限元模型和主梁3種典型截面

3．2 CFD數(shù)值模擬

青馬橋的橋板壓力分布計算采用商用軟件ANSYS FLOTRAN進行計算。為了模擬定常流風洞試驗，分析域和邊界條件設定如圖3所示。其中，流體分析域采用二維流體單元FLUID141。在入口邊界上，來流均勻分布，湍流強度小于1%。青馬橋主梁主要有3種截面形狀，對應有限元模型分別是主跨12節(jié)點截面、青衣塔10節(jié)點截面和8節(jié)點截面（圖2）。其中，12節(jié)點橋板截面為主梁的主要構型，計算網(wǎng)格劃分如圖5所示，另外2個截面網(wǎng)格劃分與其相似。

圖3 計算域與邊界條件

圖4 12節(jié)點橋板截面CFD計算模型

3．3 計算參數(shù)

將主梁橋板沿長度方向分為120個18 m長的節(jié)段，根據(jù)CFD數(shù)值模擬獲得的橋板壓力分布，由式3－6計算橋梁精細有限元模型節(jié)點抖振力。其中，橋板水平和豎直方向脈動風速時程采用諧波合成法進行模擬［15］。風速譜分別采用Simiu順向風譜和Lumley－Pnofsky豎向風譜，即：

圖5 青馬大橋主跨中點脈動風速時程（平均風速U＝18 m／s）

式中，u＊為摩擦速度，該文為1．15 m／s；n為頻率，單位是Hz，f（z）＝nz／U為折減頻率。主梁離海平面平均高度z＝60 m，模擬的時程采樣頻率為S50 Hz，頻率上限為25πHz，頻率段數(shù)為214，時長為10 min；Davenport形式的相關函數(shù)中的衰減因子取為16。在本算例中，分別對強風和一般風速下的橋梁應力響應進行分析，取主梁的平均風速分別為40 m／s、18 m／s和8 m／s，圖5是平均風速為18 m／s時，主跨中點順風方向u（t）和豎直方向w（t）的風速時程。

對青馬橋局部應力進行數(shù)值分析時，采用模態(tài)疊加法，用MSC．Nastran對橋梁精細有限元模型進行動力分析，主要考慮前80階位移模態(tài)和應力模態(tài)，模態(tài)阻尼比為1%。采用Newmark－β法求解廣義位移q（t），其中β＝0．25。

3．4 應力響應結果及分析

在青馬大橋精細有限元模型中，共有15 904個梁單元模擬主梁桁架。通過對所有梁單元應力計算結果分析比較后，發(fā)現(xiàn)青馬大橋局部應力較大的位置主要集中在：橋墩 M2（139．5 m）、馬灣側橋塔（495 m）、以及青衣側橋塔（1 872 m）（如圖6所示）。圖7所示為馬灣側橋塔的應力較大的桿件單元，主要是主梁桁架的下部順橋向、主跨一側的水平桿件，以及下部橫向桿件。其他兩個位置的重點單元部位與馬灣側橋塔相似，但桿件應力值與重點桿件的數(shù)量均不及馬灣側橋塔位置。

圖6 青馬大橋危險部位

圖7 馬灣側橋塔重點單元位置

表1是馬灣側橋塔的主梁桁架重點桿件，單元中部的應力標準差計算值，包括了梁單元截面上部和下部。從應力標準差結果可以看出，不同的平均風速下，單元應力隨平均風速增加而增大，并呈現(xiàn)平方增長趨勢，因此，在強風或臺風時，更加需要關注局部桿件的應力變化。

外側桁架下部的順橋向水平桿件單元（單元34111）的單元應力最大。在平均風速40 m／s，即通過青馬橋監(jiān)測數(shù)據(jù)推算的120 a一遇最大風速，可達6．53 MPa，因此，在抖振力的長期作用下，更容易產(chǎn)生疲勞損傷，在監(jiān)測維護時，需要對這類桿件進行重點監(jiān)護，圖8為單元34111在平均風速18 m／s時，截面下部單元應力時程。

該文計算的是梁單元中部應力值，在進行疲勞分析時，還需要考慮桿件連接部位應力集中問題。通過桿件單元應力計算結果，可確定危險桿件位置，然后參考鋼結構疲勞規(guī)范獲得應力集中系數(shù)，或者對連接部位建立更具體的局部有限元模型，進行熱點應力分析。

圖8 單元34111截面下部單元應力時程

4 結論

采用CFD數(shù)值模擬方法替代傳統(tǒng)風洞試驗，獲得大跨度橋梁主梁表面壓力分布，確定節(jié)點抖振力，在香港青馬大橋精細有限元模型上，求解了橋板局部應力響應。在不同風速下，對計算整個主梁桁架單元應力，確定了青馬橋應力較大、需要關注的重點位置和重點桿件單元，這些單元均處于橋塔與主梁連接部位。這些關鍵桿件的具體位置及單元應力響應的確定，為橋梁抖振引起的疲勞分析以及結構健康監(jiān)測設計奠定了基礎。

［1］項海帆．現(xiàn)代橋梁抗風理論與實踐［M］．北京：人民交通出版社，2005．

［2］SUN D K，XU Y L，KO J M．Fully coupled buffeting analysis of long－span cable－supported bridges：formulation［J］．Journal of Sound and Vibration，1999，228：569－588．

［3］XU Y L，KO J M，ZHANG W S．Vibration studies of Tsing Ma suspension bridge［J］．Journal of Bridge Engineering，1997，2（4）：149－156．

［4］ZHU L D，XU Y L．Buffeting response of long－span cable－supported bridges under skew winds．Part 1：theory［J］．Journal of Sound and Vibrations，2005，281：647－673．

［5］ZHU L D，XU Y L．Buffeting response of long－span cable－supported bridges under skew winds．Part 2：case study［J］．Journal of Sound and Vibrations，2005，281：675－697．

［6］SHINICHI KURODA．Numerical simulation of flow around a box girder of a long span suspension bridge［J］．Journal of Wind Engineering and Industrial Aerodynamics，1997（67／68）：239－252．

［7］SHUJI SHIRAI， TOSHIO UEDA． Aerodynamic simulation by CFD on flat box girder of super－long－span suspension bridge［J］．Journal of Wind Engineering and Industrial Aerodynamics，2003，91：279－290．

［8］SHIGERU WATANABEA，KOICHIRO FUMOTO．Aerodynamic study of slotted box girder using computational fluid dynamics［J］．Journal of Wind Engineering and Industrial Aerodynamics，2008，96：1885－1894．

［9］XU Y L，GUO W W，CHEN J，et al．Shum and H．Xia，Dynamic response of suspension bridge to typhoon and trains．II：numerical results［J］．Journal of Structural Engineering，2007，133（1）：12－21．

［10］LI Z X，ZHOU T Q，CHAN T H T，et al，Multi－scale numerical analysis on dynamic response and local damage in long－span bridges ［J］． Engineering Structures，2007（29）：1507－1524．

［11］LIU T T，XU Y L，ZHANG W S，et al．Buffetinginduced stresses in a long suspension bridge：structural health monitoring oriented stress analysis［J］．Wind and Structure，2009，12（6）：479－504．

［12］XU Y L，LIU T T，ZHANG W S．Buffeting－induced fatigue damage assessment of a long suspension bridge［J］．International Journal of Fatigue，2009，31：575－586．

［13］約翰D．安德森．計算流體力學基礎及其應用［M］．吳頌平，劉趙淼．譯．北京：機械工業(yè)出版社，2007．

［14］BLAZEK J．Computational Fluid Dynamics Principles and Applications［M］．Oxford：Elsevier Science Ltd．

［15］CAO Y H，XIANG H F，ZHOU Y．Simulation of stochastic wind velocity field on long－span bridges［J］．Journal of Engineering Mechanics，2000，126（1）：1－6．