祝志文

(湖南大學土木工程學院， 湖南 長沙 410082)

引 言

目前斜拉橋的主跨跨度已達千米量級，隨著斜拉橋跨度的不斷增大，斜拉索數(shù)量越來越多，長度也越來越大。由此導致作用在斜拉索上的氣動力也越來越大，甚至可能超過作用在主梁上的氣動力?，F(xiàn)有斜拉橋顫振穩(wěn)定性分析，基本上局限在主梁自激力的模擬[1]，不考慮斜拉索可能的自激力作用，而在橋梁顫振臨界風速附近的高風速條件下，是否需要補充因拉索振動而產(chǎn)生的氣動力，以合理評價拉索氣動力和自激作用對超大跨度斜拉橋顫振穩(wěn)定性的影響，目前并不清楚。

顫振導數(shù)的獲得，通常是借助風洞試驗或計算流體動力學數(shù)值模擬。由于拉索外形為圓柱形，其繞流流態(tài)和氣動力特性對流動雷諾數(shù)的變化非常敏感，比如一般認為雷諾數(shù)在1.2×105～4×105為圓柱的臨界雷諾數(shù)區(qū)，在該區(qū)域，隨著雷諾數(shù)的增大，阻力系數(shù)快速減小[2]，而在雷諾數(shù)高于4×105的超臨界區(qū)，阻力系數(shù)又隨雷諾數(shù)的增大而逐漸增大。因此，顫振導數(shù)的識別必須保證拉索流動雷諾數(shù)的基本一致性。

在顫振臨界風速附近，拉索的流動雷諾數(shù)一般在5×105以上，由于拉索一般以低階頻率振動，且超長拉索的基頻一般非常低，比如蘇通大橋主跨所有拉索的基頻均小于1 Hz[3]。對實橋拉索在高風速下的振動，其振動對應的折算風速將非常高，可能超過無量綱值1 000。另外，在高風速下，作用在拉索上的平均氣動力可能顯著大于自激氣動力，加之風洞試驗模型慣性力的影響，這樣，對高風速下拉索顫振導數(shù)的識別，因需同時滿足雷諾數(shù)和折算風速的一致，這無論對自由振動風洞試驗還是強迫振動風速試驗均提出了極其困難而難于實現(xiàn)的要求。

采用CFD數(shù)值模擬顫振導數(shù)，一般是根據(jù)折算風速的要求，確定計算的來流風速和模型強迫振動頻率，通過確定合理的計算參數(shù)，可獲得作用在振動模型上的氣動力，進而確定對應折算風速下的顫振導數(shù)[4]。由于模型采用強迫振動且獲得的氣動力中沒有模型運動的慣性力，因而借助CFD的數(shù)值方法幾乎成為識別超臨界雷諾數(shù)下拉索顫振導數(shù)，以及評價氣動自激力特性的唯一途徑。本文以國內某大跨度斜拉橋某根拉索為研究對象，通過數(shù)值計算包含強迫振動拉索的計算域，研究了拉索順風向振動對應的大范圍折算風速內的拉索自激力特性，并識別與拉索阻尼項有關的顫振導數(shù)。

1 拉索顫振導數(shù)的準定常解

在超臨界雷諾數(shù)下斜拉索的振動實際可能為高折算風速振動，其特征是結構運動速度相對來流風速非常小，因而在拉索運動過程中其運動速度所產(chǎn)生的相對攻角變化非常小，因而可采用準定常假設來描述其運動過程中的非定常氣動力。

(1)

圖1 振動拉索上的非定常氣動力

上式已假設拉索振動速度遠小于來流風速，此時風軸坐標系下作用在拉索上的非定常氣動阻力和升力可分別表示為：

(2a)

(2b)

式中ρ為空氣密度；CL(t)，CD(t)分別為風軸坐標系下非定常升力和阻力系數(shù)，與拉索運動產(chǎn)生的相對攻角有關。在高風速下，拉索振動速度遠小于來流風速，設α(t)為拉索合速度與來流風的攻角，因此相對攻角可表示為

(3)

且有sinα(t)≈α(t)；cosα(t)≈1。

如將式(2)轉化到以拉索中心為原點的體軸坐標系OXY下，則體軸坐標系下的升力和阻力可表示為：

FV=FDsinα(t)+FLcosα(t)

(4a)

FH=FDcosα(t)-FLsinα(t)

(4b)

如考慮拉索順風向振動及對應的阻力，且考慮拉索在高風速下的振動為其平衡位置附近的小幅振動，則可將(4b)式在零攻角附近作泰勒展開，考慮圓形截面阻力、升力及其導數(shù)特性，并忽略高階小量，有

(5)

(6)

式中K=ωD/U∞為拉索運動的折算頻率；ω=2πf為拉索振動圓頻率；f為振動頻率。

2 數(shù)值方法

2.1 控制方程

繞圓柱形拉索斷面的非定常二維不可壓流動可用下面的雷諾時均Navier-Stokes方程來描述：

(7)

(8)

上述雷諾應力的引入使得控制方程不封閉，需要引入湍流模型得以求解。如果基于渦粘假設，可將雷諾應力表示為

(9)

式中μt=ρCμk2/ε為湍流粘性;Cμ為經(jīng)驗常數(shù);k和ε分別為湍動能及耗散率，需要通過求解湍流模型方程來確定。由于標準k-ε模型往往過高估計了流動滯點區(qū)域的湍動能，一般認為，它不能用于風工程問題的數(shù)值模擬[6]。本文采用綜合了標準k-ε模型和k-ω湍流模型的SSTk-ω湍流模型。SSTk-ω湍流模型利用了標準k-ε模型適合剪切層模擬而k-ω模型適合近壁區(qū)模擬的優(yōu)點，從而通過設定一個混合函數(shù)，使得k-ω模型能在邊界層內靠近壁面使用，而邊界層外使用k-ε模型求解。相關研究認為，對于分離點附近邊界層內的非平衡流動，SSTk-ω能給出明顯優(yōu)于標準k-ε模型的模擬結果[6]。

2.2 計算域和計算參數(shù)描述

圖2 計算域分區(qū)

圖3 拉索周圍的網(wǎng)格

Z1區(qū)和Z2區(qū)均采用結構化貼體正交四邊形網(wǎng)格，Z3區(qū)為非結構四邊形網(wǎng)格。拉索表面等分為140個網(wǎng)格，為進行網(wǎng)格無關性檢查，貼近該表面的第一個網(wǎng)格點到物面的距離h0分別為5×10-6，1×10-5和2×10-5m，分別稱之為最細網(wǎng)格G1、細網(wǎng)格G2和粗網(wǎng)格G3。沿物面外法向，3套網(wǎng)格均采用1.06的網(wǎng)格生長率(文獻[6]建議的網(wǎng)格增長率不大于1.15)，以保證在物面附近流動變量變化梯度大的位置獲得高的網(wǎng)格分辨率。Z2和Z3區(qū)的公共邊兩側網(wǎng)格尺度也保持平緩變化。對3套網(wǎng)格系統(tǒng)，Z3區(qū)的網(wǎng)格劃分完全相同，不同的只是Z1和Z2區(qū)的網(wǎng)格數(shù)量和網(wǎng)格分辨率。這樣處理，在拉索周圍和尾跡區(qū)的大范圍內能獲得高質量的正交網(wǎng)格。3套網(wǎng)格系統(tǒng)各自的總單元數(shù)N見表1。

表1 網(wǎng)格劃分參數(shù)

對計算域外邊界，入口處采用自由流速度條件，水平向速度等于來流速度，垂直水平向的速度等于零，來流湍流度取為0.25%。在計算域出口采用流動出口條件，即在出口邊界上沿垂直于該邊界的法線方向，速度梯度等于零。計算域的上、下邊界均采用對稱邊界條件，即垂直于該邊界的速度為零，其它流動變量以該邊界內外分別對稱。在拉索表面，采用無滑移邊界條件，拉索表面粗糙高度與拉索直徑D的比值取3×10-5[2]。

2.3 網(wǎng)格無關檢查

對拉索繞流的非定常計算，控制方程的時間離散采用二階隱式格式，空間離散采用二階迎風格式。壓力方程和動量方程的解耦采用SIMPLEC算法和欠松弛迭代。在網(wǎng)格無關和時間步長無關檢查中，來流風速在所有工況中維持為63 m/s，對應一個恒定的拉索繞流流動雷諾數(shù)5.18×105。為了有效地捕捉流動的非定常特性，非定常計算的時間步長一般要求為圓柱繞流漩渦脫落周期的1/300～1/500[7]，為此，本文通過試算大致確定了漩渦脫落周期，并設定在所有計算工況下的時間步長均為0.000 02。對每一類時間步的數(shù)值計算，每一子步進行50次迭代，當動量方程和湍流方程的殘差小于10-5時，認為這一時間步的迭代計算已經(jīng)收斂。氣動力系數(shù)和其它參數(shù)的獲取，是在足夠多的時間步進數(shù)值計算，充分剔除初始計算影響，即氣動力系數(shù)時程和趨勢基本穩(wěn)定后開始記錄的。

定義拉索截面的非定常阻力系數(shù)CD、升力系數(shù)CL和扭矩系數(shù)CM分別為：

(10)

式中FD，F(xiàn)L和M分別對應作用在拉索截面上的阻力、升力和扭矩，其中阻力順來流流向、升力垂直來流方向向上、扭矩以拉索順時針轉動為正。

圖4中0.3 s以前的時程是固定拉索采用G1網(wǎng)格系統(tǒng)計算得到的氣動力系數(shù)時程。因受初始計算的影響，要經(jīng)過0.1 s大概5 000個時間步計算，氣動力時程數(shù)據(jù)才表現(xiàn)為有規(guī)律的周期數(shù)據(jù)。升力系數(shù)的主頻率為阻力系數(shù)的二分之一，也即對應拉索的漩渦脫落頻率。為保證固定拉索氣動力系數(shù)完全擺脫初始計算的影響，本文對上述3套網(wǎng)格系統(tǒng)的繞流計算，均模擬了0.3 s計15 000個時間步。由上述3套網(wǎng)格系統(tǒng)計算得到的阻力和升力系數(shù)平均值和均方根值見圖5，這些值是根據(jù)0.1～0.3 s的時程數(shù)據(jù)統(tǒng)計得到。

圖4 拉索氣動力系數(shù)時程

圖5 拉索阻力和升力系數(shù)平均值和均方根值

從圖5可見，固定拉索的升力系數(shù)平均值非常小，這是因為拉索外形上下對稱(其流態(tài)上下并不對稱)，其RMS值顯著大于其平均值。與此相反，阻力系數(shù)的平均值顯著大于其RMS值。對3套網(wǎng)格系統(tǒng)，隨著物面網(wǎng)格尺度的減小，氣動力系數(shù)值有稍微的變化，表現(xiàn)為阻力系數(shù)平均值和升力系數(shù)RMS值有少量的增大，但從G2到G1網(wǎng)格的變化不到1%。阻力系數(shù)RMS和升力系數(shù)平均值在3套網(wǎng)格系統(tǒng)上基本沒有變化。因此可以認為，從網(wǎng)格系統(tǒng)G2到G1，數(shù)值模擬結果沒有明顯的變化，可認為已獲得了與網(wǎng)格無關的解。因從網(wǎng)格系統(tǒng)G2到G1，網(wǎng)格數(shù)量沒有顯著增加，因此，在后續(xù)運動拉索數(shù)值模擬中，本文均采用網(wǎng)格系統(tǒng)G1開展研究。另外，由G1網(wǎng)格系統(tǒng)得到的阻力系數(shù)平均值為CD=0.81，可作為準定常理論需要的阻力系數(shù)值，這與超臨界雷諾數(shù)下拉索阻力系數(shù)建議值0.8基本吻合[2]，表明了數(shù)值方法的有效性，因而可將G1網(wǎng)格和所用計算參數(shù)來開展后續(xù)拉索自激力特性研究。

3 結果和分析

為獲得拉索在來流中運動時，作用在拉索上的氣動力，CFD模擬采用單自由度強迫振動法，即強迫拉索作順風向單自由度、單頻率諧振動，拉索強迫運動位移按下式給定

將甘薯淀粉（SPS，sweet potato starch）與魔芋膠（KGM，konjac gum）按以下比例混合（質量比10：0，9.5：0.5，9.0：1.0，8.5：1.5，8.0：2.0），準確稱取各配比下的甘薯淀粉、魔芋膠于燒杯中加入去離子水充分混合，配制成質量分數(shù)為8%的均一懸浮液（以干基計），于沸水浴中攪拌、加熱糊化15 min。除糊化特性外，老化特性、流變學特性的測定均采用該方法制備樣品。

p(t)=p0sinωt

(11)

式中p0為順風向振動位移幅值，通常幅值越大，自激氣動力也增大，但如振動幅值很大時，顫振導數(shù)識別假定的線性小擾動前提將不滿足，可能會出現(xiàn)較大的流動非線性；另外，在CFD模擬中，因在每一個時間步后要重新劃分網(wǎng)格，因此振動幅度越大，計算網(wǎng)格的變形幅度也越大，網(wǎng)格的畸變程度越高，這會影響CFD的計算精度?？紤]來流風速較大，為增大低振動頻率下自激氣動力對拉索渦脫力的比值，并權衡模型運動導致的氣動力非線性，本文取p0=0.1D，并在Fluent中確認了該振動幅值下網(wǎng)格的畸變不大。

折算風速定義為Vr=U∞/(fD)。因此拉索不同折算風速下的流動模擬僅需通過改變強迫拉索振動的頻率，從而得到不同的折算風速[8]。由于來流風速和拉索尺寸不變，因此不同的折算風速模擬仍然保持了相同計算雷諾數(shù)，從而使得所有的數(shù)值計算條件，如網(wǎng)格、時間步長等參數(shù)，在不同折算風速下的模擬保持一致，保證了通過網(wǎng)格無關檢驗確定的參數(shù)能在所有折算風速下完全一致性。

對每一個折算風速下的模擬，為使拉索強迫振動模擬不受初始計算的影響，首先對固定拉索均進行0.3 s共計15 000個時間步的數(shù)值計算，大量時間步計算使得拉索繞流能形成穩(wěn)定的漩渦脫落狀態(tài)，這可從氣動力系數(shù)的規(guī)則振蕩看出，這也是所有折算風速模擬的共同初始條件。然后強迫拉索按給定的位移模式作單自由度運動，振動拉索繞流通過Fluent的動網(wǎng)格實現(xiàn)，在每一時間步計算完成后更新計算域網(wǎng)格。

圖4顯示了順風向振動時作用在拉索上的氣動力系數(shù)時程，0.3 s后為拉索按10 Hz頻率強迫振動的2個周期計算結果，對應時間步長的采樣頻率為50 kHz。從固定拉索狀態(tài)啟動強迫運動，CFD求解在0.3 s后幾個時間步上有較劇烈的數(shù)值振蕩，這主要表現(xiàn)在阻力系數(shù)時程上有突躍，這是模型突然運動，導致速度導數(shù)不連續(xù)導致。隨后該現(xiàn)象消失，氣動力呈現(xiàn)有規(guī)律振蕩，表現(xiàn)為拉索渦脫氣動力信號受到低頻強迫振動信號的調質，且強迫振動的頻率越低，氣動力的調幅幅度越小，而振動前后阻力和升力系數(shù)時程的幅值變化并不大。從渦脫力和自激氣動力的幅值來看，阻力和升力系數(shù)的渦脫力明顯大于自激力，即渦脫阻力和升力主導氣動力；從扭矩系數(shù)來看，因固定拉索繞流的扭矩系數(shù)非常小，因而強迫振動顯著增大了拉索的扭矩系數(shù)。從振動拉索氣動力信號的頻域特征分析來看，拉索振動并沒有改變拉索繞流自身的漩渦脫落頻率。

為了獲得拉索強迫振動的自激氣動力，需要將總氣動力中的渦脫力剔除出來，本文假設小幅振動下氣動力滿足可疊加性。由于自激力頻率明顯低于渦脫力，因而可采用低通濾波的方法，將高頻渦脫力過濾掉。圖6是圖5振動拉索在0.3～0.4 Hz的一個完整振動周期內的總氣動三分力，以及采用截止頻率為15 Hz的低通濾波器過濾后的氣動力，對應的折算風速為53。因為過濾所得信號的頻率仍為10 Hz，因此本文將其視為由強迫拉索振動產(chǎn)生的自激氣動力。從圖6的3個氣動力系數(shù)時程幅值，可見氣動力耦合仍然存在。另與總氣動力值相比，自激氣動力均非常小，但相比而言，沿拉索運動方向的自激阻力系數(shù)幅值最大，因而自激力效應將主要體現(xiàn)在與運動一致的自由度方向。

圖6 強迫振動總氣動力與自激氣動力系數(shù)時程(Vr=53)

圖7是強迫振動頻率為0.5 Hz，拉索在0.3～2.3 s的一個完整振動周期內的總氣動三分力，以及采用截止頻率為1 Hz的低通濾波器過濾后得到的氣動力時程，對應的折算風速為1 050。與上圖對比可見，在超高折算風速下，拉索自激氣動力變化幅值明顯小于低折算風速值，特別是沿運動方向的自激阻力系數(shù)，這說明隨著折算風速的增大，拉索運動產(chǎn)生的自激氣動力在總氣動力總的比重將減小，本文計算的折算風速為1 050時自激氣動阻力系數(shù)最大值與阻力系數(shù)平均值的比值β約千分之一，因此非常小，見表2。

因此可以認為，在超臨界雷諾數(shù)下，拉索自激力顯著小于其渦脫力，且由于氣動力變化頻率特征仍主要表現(xiàn)為渦脫頻率特征，因此自激力可忽略，風荷載計算只需考慮渦脫力。

圖7 強迫振動總氣動力與自激氣動力系數(shù)時程(Vr=1 050)

圖8 計算的與理論解的對比

表2 CFD計算工況和相關計算值

圖9 拉索一個強迫振動周期內關鍵時刻的漩渦脫落圖

圖9是拉索一個強迫振動周期內的在4個關鍵時刻的漩渦脫落圖，分別對應拉索離開平衡位置到達最大振幅的時刻T/4、回到平衡位置時刻T/2、反向離開平衡位置到達最大振幅的時刻3T/4，以及完成一個周期振動回到平衡位置時刻，其中T為拉索振動周期。雖然可以看到拉索尾跡漩渦的產(chǎn)生、拉長、脫落和隨尾跡漂移。

4 結 論

本文基于CFD數(shù)值方法研究拉索順風向振動的自激力特性，得到下述結論：

3)在超高折算風速下，拉索自激氣動阻力顯著小于其阻力平均值，自激升力明顯小于升力幅值，即拉索渦脫力顯著主導氣動力，拉索振動沒有改變拉索繞流自身的漩渦脫落頻率。

4)在超臨界雷諾數(shù)和拉索橫風向振動情況下，拉索運動產(chǎn)生的自激氣動力在總氣動力總的比重非常小，對頻率特性的改變也非常小，因而自激氣動力可忽略，拉索風荷載只需考慮渦脫力。

參考文獻：

[1] Ge Y J， Tanaka H. Aerodynamics flutter analysis of cable-supported bridges by multi-mode and full-mode approaches[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2000，86：123—153.

[2] Poulin Sanne, Larsen Allan, Drag loading of circular cylinders inclined in the along-wind direction[J]. Journal of Wind Engineering and Industrial Aerodynamics,2007,95:1 350—1 363.

[3] 楊素哲, 陳艾榮, 超長斜拉索的參數(shù)振動[J]. 同濟大學學報(自然科學版), 2005,33(10):1 303—1 308.Yang Su-zhe , Chen Ai-rong, Parametric oscillation of super long stay cables[J]. Journal of Tongji University(Natural Science), 2005, 33(10): 1 303—1 308.

[4] Zhi-Wen Zhu, Ming Gu, Zheng-qing Chen. Wind tunnel and CFD study on identification of flutter derivatives of a long-span self-anchored suspension bridge[J]. Computer-Aided Civil and Infrastructure Engineering,2007,22:541—554.

[5] Scanlan R H, Tomko J J. Airfoil and bridge deck flutter derivatives [J]. Journal of Engineering Mechanics, ASCE, 1971, 97(6): 1 717—1 737.

[6] Casey M, Wintergerste T. ERCOFTAC Special interest group on “Quality and trust in industrial CFD”, Best practice guidelines[R]. Fluid Dynamics Laboratory, Sulzer Innotec, 2000.1.

[7] Claudio Mannini, Ante Soda, Ralph Voβ, et al. Unsteady RANS simulations of flow around a bridge section[J]. Journal of Wind Engineering and Industrial Aerodynamics,2010，98：742—753.

[8] 祝志文，陳政清. 數(shù)值模擬橋梁斷面的顫振導數(shù)和顫振臨界風速[J].中國公路學報，2004，15(4):41—50．Zhiwen Zhu, Zhengqing Chen. Numerical simulations for aerodynamic derivatives and critical flutter velocity of bridge deck[J]. China Journal of Highway and Transport,2004,17(3):41—50.