箱梁斷面氣動導數的CFD模擬研究

劉 鑰，張志田，陳政清

(1.中國市政工程西北設計研究院有限公司武漢分院，武漢 430056;2.湖南大學 風工程試驗研究中心，長沙 410082)

按照風洞試驗模型的振動情況來分類，目前所采用的氣動導數試驗識別方法可分為自由振動法與強迫振動法兩大類。近年來，隨著計算機技術及計算流體動力學的發展，采用數值仿真的方法識別橋梁氣動導數已成為可能。雖然與風洞試驗相比，數值方法還存在一定局限性，但數值方法已成為識別橋梁氣動導數的另一種方法。Walther借助離散渦法，在氣動彈性數值分析方面做了開創性的工作［1］;Larsen基于離散渦法開發了DVMFLOW軟件，數值計算了從流線型至H型逐漸變化的5種典型橋梁斷面的氣動導數［2］;曹豐產用有限單元法求解二維不可壓粘性流體的N-S方程計算運動物體的氣動力，識別了幾種典型橋梁斷面的氣動導數［3］;周志勇采用離散渦法，提取了南京二橋閉口箱梁斷面的氣動導數［4］;祝志文通過有限體積法對橋梁斷面的氣動導數及顫振臨界風速進行了數值模擬［5］，張倩采用計算軟件FLUENT，在橋梁氣動導數的計算方面提出了剛性邊界層網格+三角形動網格+靜止網格的網格分區方法［6］;辛大波分別采用了Fluent中的RNGk-ε和Realizablek-ε兩種湍流模型識別了氣動導數［7］;劉慕廣采用動網格強迫結構簡諧振動，識別了薄平板、大帶東橋及腹板開孔H型截面的氣動導數［8］。

本文在對某大橋跨中箱梁斷面0°攻角下的氣動導數進行數值識別前，先通過商用CFD軟件，采用動網格法強迫模型運動的方式對薄平板的氣動導數進行了數值識別并與試驗結果進行了對比，驗證CFD數值識別的可行性，隨后以跨中箱梁斷面模型為原型，提取了8個氣動導數，并與強迫振動試驗法得到的氣動導數結果進行了對比。

1 薄平板氣動導數

圖1 薄平板模型(mm)Fig.1 Model of thin flat plate(mm)

薄平板計算模型尺寸如圖1所示，寬高比為22.5，模型計算域為二維矩形，入口距平板模型中心為5B(B為平板模型寬)，出口距模型中心為10B，上下距模型中心各為5B。不同雷諾數時數值試驗研究所需要的網格精細程度也不同，主梁壁面的第一層網格需要滿足y+≈1。可以用下面的公示估算y+，并確定第一層網格控制點離壁面的距離 Δy(Schlichting，1979):y+=0.172Δy/LRe0.9。由于主梁斷面在風洞試驗時很難知道雷諾數的大小，于是采用 Δy試算法得出Δy已滿足y+≈1這一條件。采用ANSYS ICEM CFD 10.0生成計算網格，構造網格時僅采用結構網格，靠近物面的第一條網格線距物面距離為Δy=0.09 mm，即0.000 2B，以1.3的生長率徑向逐漸放大網格(圖2)，共生成了2.82萬左右的結構網格。

采用ANSYS CFX 10.0進行邊界條件給定及數值計算。入口為速度邊界條件，即Vx=6 m/s，Vy=0 m/s;出口為壓力邊界條件，參考壓力為零;上下面采用自由滑移壁面條件，即Vx=6 m/s，Vy=0 m/s;主梁斷面結構表面采用無滑移壁面條件，即Vx=0 m/s，Vy=0 m/s;側壁采用對稱邊界條件。

模型作豎向振動的振幅為0.02B，作扭轉運動的振幅為2°。首先固定平板模型做定常計算，然后以此為初場使平板強迫振動，進行非定常計算，計算中采用SST湍流模型。通過變化模型振動頻率的方式來提高無量綱風速，計算采用的時間步長為0.002 s，共計算得到了10 s內的氣動力時程曲線。

圖2 計算網格Fig.2 Meshing of thin flat plate

圖3 數值模擬時的氣動力時程曲線(U/fB=10)Fig.3 Aerodynamic forces time-history curves in CFD(U/fB=10)

圖3為無量綱風速U/fB=10，薄平板分別數值模擬時做豎向與扭轉運動前10 s時氣動力時程，圖4和圖5為平板豎向運動和扭轉運動時壓力等高線及流線圖。由圖3可見，薄平板強迫運動產生的氣動自激力諧波特性極好。由圖4中的流線圖看見，平板模型在豎向運動過程中無明顯尾流分離區，尾流中也沒有旋渦的脫落與運動，這跟平板寬高比為22.5，較好的流線性有關。

從圖4和圖5中的薄平板時程升力曲線最后一周期內渦脫壓力等高線圖可以看出，在t=0時，薄平板上頂面和下底面的壓力等高線相互對稱，上下壓力相互抵消，此時薄平板對應的升力為平均值，其值接近于0;在t=1/4T時刻，下底板的絕對值大于上頂板，此時升力達到最大值;在t=1/2T時刻跟t=0時刻一樣，此時升力為平均值;在t=3/4T時刻，下底板的絕對值小于上頂板，升力達到最小值。從圖中我們可以看出t=0時刻跟1/2T時刻、t=1/4T時刻跟t=3/4T時刻對應的壓力等高線恰好成鏡像關系。

圖6中的試驗結果為采用強迫振動法識別出的平板氣動導數值，由圖識別結果可以看出:

(4)由以上平板數值方法識別結果可以看出，采用CFD動網格法識別顫振導數的合理性及采用數值方法識別橋梁顫振導數可行。

2 橋梁箱形斷面的氣動導數

某大橋是一混凝土剛構橋，主梁跨中斷面模型寬B=0.3 m，數值計算時沒有考慮欄桿等附屬物。計算域同樣為二維矩形，入口距模型剪心為5B(B為橋面寬)，出口距剪心為10B，上下距剪心各為5B。靠近物面的第一條網格線距物面距離為Δy=0.000 5 m，即0.001 7B，滿足y+≈1這一條件，并以1.3的生長率徑向逐漸放大網格(圖6)，共生成了3.02萬左右的結構網格。

計算域入口為邊界條件與前面平板計算模型設定一致。模型做豎向振動的振幅為0.015B，做扭轉運動的振幅為2°。首先固定模型做定常計算，然后以此為初場，進行非定常強迫振動計算，湍流模型為SST模型。通過變化模型振動頻率的方式來提高無量綱風速U/fB，計算采用的時間步長為0.002 s，共計算得到了前10 s的時程曲線。

由圖8可看出，氣動力波形沒有薄平板那樣有規則，由于跨中斷面分離流較為復雜，氣動力信號中含有較多的干擾成份。寬高比為3.43的跨中斷面與薄平板模型明顯的區別在于，跨中斷面存在一個尾流分離區，且分離區的大小及位置隨斷面運動而變化，故跨中斷面的壓力等高線分布比平板也要復雜一些。圖9和圖10為斷面豎向和扭轉運動時對應一個周期內壓力等高線及流線變化。

從圖9中的跨中斷面升力時程曲線最后一周期內渦脫壓力等高線圖和流線圖可以看出，在t=0時，箱梁上頂板和下底板附近的尾流出現速度旋渦，此時對應地在下底板后緣的尾流出現壓力旋渦;在t=1/4T時刻，箱梁上頂板出現速度旋渦，尾流的速度旋渦向下游運動，此時對應地在離下底板較遠的地方出現壓力旋渦，且在上頂板后邊緣有形成壓力旋渦的趨勢;在t=1/2T時刻，箱梁上頂板和下底板附近的尾流出現旋渦，此時對應地在上頂板和下底板的后緣出現壓力旋渦;在t=3/4T時刻，箱梁上頂板和下底板附近的尾流出現旋渦，下底板后緣的旋渦逐漸增大并有向下游運動的趨勢，此時對應地在上頂板和下底板的后緣出現壓力旋渦，兩個旋渦都向下游運動。

從圖10中跨中斷面升力時程曲線最后一周期內渦脫壓力等高線圖和流線圖可以看出，在t=0時，離箱梁較遠的尾流出現速度旋渦，此時對應地在后緣的尾流出現壓力旋渦;在t=1/4T時刻，箱梁上頂板后緣出現速度旋渦，尾流的速度旋渦向下游運動，此時對應地在上頂板后緣和離下底板較遠的地方出現壓力旋渦;在t=1/2T時刻，箱梁上頂板和下底板附近的尾流出現旋渦，此時對應地在上頂板和下底板的后緣出現壓力旋渦;在t=3/4T時刻，箱梁下底板后緣和離上頂板較遠的尾流出現旋渦，此時對應地在箱梁下底板后緣和離上頂板較遠的尾流出現壓力旋渦。

圖9 箱形斷面豎向運動時的壓力等高線與流線Fig.9 Streamlines and pressure contour lines around box girder within a cycle as the vertical movement of mesh

圖10 箱形斷面扭轉運動時的壓力等高線與流線Fig.10 Streamlines and pressure contour lines around box girder within a cycle as the torsion movement of mesh

表1 某大橋箱形斷面的氣動導數Tab.1 Aerostatic derivatives?of one long －span bridge’s box girder

表2 某大橋箱形斷面的氣動導數Tab.2 Aerostatic derivatives of one long － span bridge’s box girder

跨中斷面氣動導數的試驗結果如表1所示。跨中斷面氣動導數的數值模擬結果如表2所示。數值模擬方法識別到的跨中斷面氣動導數如圖11所示，圖中試驗值為強迫振動法得到的氣動導數結果。

圖11 箱形斷面顫振導數Fig.11 Aerostatic derivatives of box girder

3 結論

本文通過CFD數值模擬結合強迫振動風洞試驗識別了薄平板、主梁跨中斷面的氣動導數。結果表明:

(1)對流線性好的薄平板，CFD數值計算氣動導數值與理論解曲線吻合良好。

(2)鈍體箱形斷面數值識別結果與試驗值存在一定誤差。由于分離流與旋渦的復雜性，數值模擬結果在低折減風速更接近試驗值，而在高折減風速，差別較為明顯。識別精度與數值計算中網格密度及時間步長緊密相關。

(3)跨中斷面識別的氣動導數與試驗值吻合程度較好，結合薄平板的識別結果可見，數值方法強迫模型運動提取氣動導數的過程受模型本身流線程度和離散網格質量的影響較大，斷面流線型越好，網格質量越高，數值模擬識別得到的氣動導數精度就越高。

(4)從氣動導數的趨勢來看，CFD計算結果與風洞試驗有較好的一致性。但對于鈍體特性十分明顯的橋梁箱形斷面，目前采用CFD方法識別氣動導數存在明顯不足之處，尤其是算法的穩定性及效率問題。因此，對于鈍體特性十分明顯的斷面，CFD方法目前尚不能完全取代物理風洞研究其流固耦合問題。

［1］ Walther J H，Larsen A.Two dimensional discrete vortex method for application to bluff body aerodynamics［J］.Journal of Wind Engineering and Industrial Aerodynamics，1997，67-68:183-193.

［2］Larsen A，Walther J H.Aeroelastic analysis of bridge girder sections based on discrete vortex simulations［J］.Journal of Wind Engineering and Industrial Aerodynamics，1997，67-68:253-265.

［3］曹豐產，項海帆，陳艾榮.橋梁斷面氣動導數和顫振臨界風速的數值計算［J］.空氣動力學學報，2000，18(1):26-33.

［4］周志勇，陳艾榮，項海帆.渦方法用于橋梁斷面氣動導數和顫振臨界風速的數值計算［J］.振動工程學報，2002，15(3):327-331.

［5］祝志文，陳政清.數值模擬橋梁斷面氣動導數和顫振臨界風速［J］.中國公路學報，2004，17(3):41-45.

［6］張 倩，廖海黎.橋梁主梁斷面氣動力特性的數值模擬研究［J］.工程結構，2008，28(4):83-86.

［7］辛大波.基于數值方法的大跨橋梁顫振分析［D］.哈爾濱:哈爾濱工業大學，2004，1-82.

［8］劉慕廣.兩類大長細比橋梁構件的風振特性研究［D］.長沙:湖南大學，2009，1-169.