王 劍，畢繼紅，何旭輝，關(guān) 健，喬浩玥

(1. 天津城建大學 天津市土木建筑結(jié)構(gòu)防護與加固重點實驗室，天津 300384；2.天津大學 建筑工程學院, 天津 300072；3.中南大學 土木工程學院橋梁工程系，長沙 410000)

隨著現(xiàn)代斜拉橋跨徑的不斷增大，斜拉索越來越長、剛度和阻尼越來越小，導致其在特定風速及降雨條件下產(chǎn)生風雨激振現(xiàn)象。自1988年Hikami等[1]首次觀測到風雨激振現(xiàn)象以來，研究人員對風雨激振的機理進行了大量的研究，并認識到拉索表面水線的形成和振蕩是風雨激振的主要原因。對此Lemaitre等[2-3]等將滑移理論應用到風雨激振研究中，用風壓力和風摩擦力來表示氣流對水膜的作用，模擬水平靜止拉索表面上的水線形成。許林汕等[4]考慮拉索運動對水膜形態(tài)變化的影響，根據(jù)已有研究成果假設拉索的運動規(guī)律，建立液固耦合模型研究水線的運動規(guī)律。Taylor等[5]則運用Lemaitre模型，通過數(shù)值求解得到隨時間變化的風壓力系數(shù)和風摩擦力因數(shù)，建立氣液耦合模型分析水平靜止的拉索表面上的水線運動規(guī)律。畢繼紅等[6-7]等將滑移理論和單質(zhì)點運動理論結(jié)合，建立了氣液固三相耦合的風雨激振模型。

國內(nèi)外研究人員通過一系列的現(xiàn)場觀測和風洞試驗研究，發(fā)現(xiàn)在降雨強度為小到中雨的條件下，斜拉索僅在一定的風速范圍內(nèi)發(fā)生風雨激振現(xiàn)象，而風速過大或過小均不會出現(xiàn)，即風雨激振有“限速”的特點。現(xiàn)階段關(guān)于風雨激振的起振風速區(qū)間的研究大多以風洞試驗為主[8-9]，理論分析方面的研究較少，特別是從水線運動角度的研究就更少。對此，本文對氣液固三相耦合模型的數(shù)值計算過程進行改進，采用有限差分法求解水膜運動方程，并且應用有限元軟件COMSOL求解隨時間變化的風壓力系數(shù)和風摩擦力因數(shù)，在保證計算精度的同時大幅度提高計算效率；在此基礎上，研究風速變化對水線運動及斜拉索振動的影響，分析風雨激振的產(chǎn)生機理。

1 模 型

圖1所示的半徑為R、傾角為α(0°≤α≤90°的斜拉索，受重力和水平方向風的共同作用，風速大小為U0，風偏角為β(0°≤β≤90°)。

圖1 斜拉索空間位置Fig.1 Spatial position of stay cable

取斜拉索任一橫截面A-A為研究對象，如圖2所示。對拉索的二維模型進行研究，忽略重力在拉索軸向方向的影響。

圖2 斜拉索表面水膜受力Fig. 2 Force of water film around cable

根據(jù)滑移理論，假設拉索表面存在一層連續(xù)的水膜，兩自由度無量綱化的水膜運動方程可表示為

(1)

其中，

(2)

式中：t為時間；μ為水的動力黏度系數(shù)；γ為水在空氣中的表面張力系數(shù)；h為水膜厚度；Rc為拉索半徑；M為單位長度拉索的質(zhì)量；ρg和ρ分別為空氣和水的密度；κ為水膜表面的曲率；n和t分別為水膜與空氣交界處的法向向量和切向向量；σg分別為空氣的應力張量；pg為水膜表面所受的空氣壓力。

斜拉索在橫風向和順風向的振動方程為

(3a)

(3b)

其中，

(4a)

(4b)

2 數(shù)值求解

2.1 有限差分格式

(5)

其中，

(6a)

(6b)

(6c)

2.2 風壓力系數(shù)Cp與風摩擦力因數(shù)Cf

水膜運動方程中的風壓力系數(shù)Cp和風摩擦力因數(shù)Cf是影響水膜運動的關(guān)鍵參數(shù)。在已有的研究中，畢繼紅等人采用Fluent軟件獲取表面帶有不同形狀水膜的斜拉索的Cp和Cf，計算時間很長，效率極低，僅能獲得前20 s的運動狀態(tài)；而Gu等[10]的試驗數(shù)據(jù)表明，從試驗開始到風雨激振真正形成，需要數(shù)十秒甚至上百秒的時間。對此，本文采用有限元軟件COMSOL計算隨時間變化的風壓力系數(shù)和風摩擦力因數(shù)，大幅提高了計算效率和計算穩(wěn)定性。

本文應用COMSOL軟件計算時仍采用與原先相同的假設，即假設水膜所受的氣流作用與干燥拉索表面所受的氣流作用相同。每一時間步內(nèi)均根據(jù)上一步計算得到的水膜形狀，將水膜考慮為固態(tài)，建立流場中的繞流障礙物，應用Spalart-Allmaras湍流模型，采用穩(wěn)態(tài)計算此時刻的Cp和Cf。

設干燥拉索直徑為d，計算區(qū)域及網(wǎng)格劃分如圖3所示。邊界條件設置如下：

左側(cè)進口處采用速度邊界條件(u=UN,v=0)；右側(cè)出口處采用壓強邊界條件(p=0)；上下邊界處采用完全滑移邊界條件；圓柱表面處采用無滑移邊界條件(u=0,v=0)。

圖3 流場網(wǎng)格劃分Fig.3 Numerical grids of the flow field

將數(shù)值求解得到的風壓力系數(shù)Cp和風摩擦力因數(shù)Cf與Fage等[11]的試驗數(shù)據(jù)、Celika等[12]及Bi等的數(shù)值計算結(jié)果進行對比，考慮圓柱直徑為0.16 m,風速為10 m/s的試驗工況，如圖4所示，應用COMSOL軟件求解得出的結(jié)果與前人的研究結(jié)果非常接近，驗證了數(shù)值方法的準確性。

圖4 圓柱表面風壓力系數(shù)Cp和風摩擦力因數(shù)Cf的分布圖Fig.4 Distribution of wind pressure coefficient Cp and friction coefficient Cf around cylinder

2.3 基本參數(shù)及數(shù)值計算工況基本

借鑒文獻[13-14]中的風洞試驗數(shù)據(jù)，基本參數(shù)選取如下：斜拉索半徑0.05 m、傾角30°、自振頻率0.952 Hz、線密度8.57 kg/m、阻尼比0.17%；拉索的初始位移、速度和加速度均為零；水膜的初始厚度為0.2 mm；水密度1.0×103kg/m3、水的運動黏性系數(shù)1.0×10-6m2/s、水在空氣中的表面張力系數(shù)7.2×10-2N/m、空氣密度1.225 kg/m3、空氣的運動黏性系數(shù)1.51×10-5m2/s；風偏角22.5°，風速6.0 m/s，6.76 m/s，7.4 m/s，7.5 m/s，7.6 m/s，7.72 m/s，8.2 m/s，8.5 m/s，9.0 m/s，9.5 m/s，10.5 m/s，12.0 m/s，13.5 m/s。

參照已有的研究成果，在采用有限差分法求解水膜運動方程時，時間步長dt=1×10-5s，為保證拉索表面水膜的質(zhì)量守恒，根據(jù)試驗觀測數(shù)據(jù)，假設水膜的最大厚度為1 mm；為使水膜在拉索表面連續(xù)分布，假設水膜的最小厚度為0.02 mm。

2.4 數(shù)值求解流程

水膜運動方程(見式(1))為四階偏微分方程，拉索振動方程(見式(3))為二階微分方程，采用差分法求解方程組，應用MATLAB軟件進行數(shù)值計算并利用有限元軟件COMSOL求解各個時間步不同水膜形態(tài)下的風壓力系數(shù)Cp和風摩擦力因數(shù)Cf，基本流程如圖5所示。

圖5 數(shù)值計算流程圖Fig.5 Flowchart of numerical calculation

3 數(shù)值計算結(jié)果

3.1 拉索振動響應

圖6顯示了各個風速工況下斜拉索的橫風向及順風向振幅變化。數(shù)值計算得到的起振風速區(qū)間是6.76 m/s

圖6 風速對拉索振動的影響Fig.6 Influence of wind velocity on cable vibrations

圖7為風速7.5 m/s下拉索在0～150 s內(nèi)的振動響應。在0～29 s以內(nèi)，拉索橫風向振幅迅速增大，振幅達到0.19 m左右，如圖6(a)所示；而后振幅迅速衰減，至50 s附近，振幅減至0.13 m；50～100 s內(nèi)拉索振幅繼續(xù)小幅減小，在100 s以后逐步趨于穩(wěn)定，保持在0.10 m左右，與實驗結(jié)果(0.09 m附近)相近，如圖6所示。

7.5 m/s時拉索的順風向振幅如圖7(b)所示。初始時刻，拉索的順風向振幅為0.14 m；前45 s內(nèi)，拉索振幅迅速減小，最小值為0.03 m；之后小幅增大，72 s后逐漸穩(wěn)定在0.05 m左右，近似等于拉索橫風向振幅的一半，與現(xiàn)場觀測結(jié)果相同。與前人的研究結(jié)論相一致，拉索在兩個方向的振動均需要幾十秒(50 s左右)的時間達到穩(wěn)定振動狀態(tài)，因此本文對各風速工況均進行了150 s的計算，選用100～150 s的計算結(jié)果進行分析研究。

圖7 風速7.5 m/s時的拉索振動時程曲線Fig.7 Vibration responses of cable with wind velocities 7.5 m/s

3.2 拉索氣動升力和阻力

各個風速工況下的拉索氣動升力及阻力變化幅度曲線，如圖8所示。隨著風速的增加，升力和阻力在總體上均呈現(xiàn)出不斷增大的趨勢。但當風速處于6.76 m/s

對每一風速下拉索的氣動力變化時程進行頻譜分析，如圖9(a)所示。絕大多數(shù)風速下升力和阻力的主頻相一致。當風速較低(U0≤6.76 m/s)時，升力和阻力的變化均無明顯的主頻；當風速位于風雨激振的風速區(qū)間(6.76 m/s

圖8 風速對氣動力變化最大幅度的影響Fig.8 Influences of wind velocity on aerodynamic forces

圖9(b)為拉索氣動力在其主頻附近的變化幅值與風速變化間的關(guān)系曲線。6.76 m/s

綜上所述，決定拉索振動強弱的主要因素是氣動力變化在拉索自振頻率附近幅值，而非其整體的最大變化幅度。

3.3 水線形態(tài)及其運動特征

隨著風速從6 m/s增大至9.5 m/s，上水線的最大厚度不斷增加，從0.22 mm逐漸增加至0.81 mm，且其增加的幅度越來越大；至10.5 m/s以上時，上水線的最大厚度達到計算假設的1.0 mm，如圖10所示。上述趨勢與拉索氣動力最大幅度的變化趨勢相近，說明上水線厚度對升力和阻力的變化幅度有顯著影響。

從圖11所示的各風速工況下的上水線振蕩平衡位置可以看出，風速為6 m/s時的上水線中心位置大致在40°附近；隨著風速的增大，上水線平衡位置不斷向背風側(cè)偏移；到13.5 m/s風速時移至約70°。文獻[13]和文獻[14]中的坐標系與本文數(shù)值計算時選取的坐標系相差90°，經(jīng)坐標變換后進行對比，發(fā)現(xiàn)數(shù)值計算結(jié)果與風洞試驗結(jié)果基本一致，如圖11所示。

圖9 不同風速拉索的氣動力變化頻譜分析Fig.9 Frequency spectrum analysis of aerodynamic forces with various wind velocities

圖10 各風速工況下的上水線厚度Fig. 10 Upper rivulet thickness with various wind velocities

圖11 各風速工況下的上水線振蕩平衡位置Fig.11 Equilibrium positions of upper rivulet oscillations with various wind velocities

圖12所示顯示了各風速工況下的上水線振蕩幅度。由于風速小于7.4 m/s時，上水線振蕩幅度很小，故在此只研究風速大于7.4 m/s時的上水線振蕩幅度。當風速位于風雨激振的風速區(qū)間(6.76 m/s

圖12 各風速工況下的上水線振蕩幅度Fig.12 Amplitudes of upper rivulet oscillations with various wind velocities

通過對各風速工況下的上水線位置處水膜厚度變化時程進行頻譜分析(見圖13)，可以進一步研究水線運動與拉索振動間的相互影響。

圖13(a)為上水線振蕩頻率與風速的關(guān)系，當風速較低(U0≤6.76 m/s)時，上水線厚度較小且基本無運動；當風速位于風雨激振的風速區(qū)間(6.76 m/s

圖13(b)顯示了各風速工況下上水線位置處的水膜厚度變化在其主頻附近的幅值。當風速位于風雨激振的風速區(qū)間(6.76 m/s

對比圖6所示的拉索振動與圖10、圖13所示的水線振蕩，可以看出影響拉索振動強弱的主要因素是水線的振動頻率而非其厚度，只有當水線振蕩頻率與拉索自振頻率相近時，才會使得拉索氣動力的頻率接近拉索自振頻率，產(chǎn)生風雨激振現(xiàn)象；此時，拉索的振動幅度與拉索自振頻率附近上水線振蕩的強弱有密切聯(lián)系。

分別選取風速6.0 m/s、7.72 m/s、9.0 m/s和13.5 m/s的水膜厚度時程變化用來分析研究水線運動與拉索振動間的相互影響，如圖14所示。

圖13 不同風速拉索的上水線頻譜分析Fig.13 Frequency spectrum analyses of upper rivulets with various wind velocities

6.0 m/s風速下，上、下水線均保持穩(wěn)定狀態(tài)，振蕩幅度很?。?.72 m/s風速下，上、下水線均做周期性的環(huán)向振蕩，周期大致與拉索的自振周期一致，且上水線的振蕩幅度遠大于下水線；9.0 m/s風速下，上、下水線亦在振蕩，但其周期顯著減小，遠小于拉索自振周期；至13.5 m/s風速時，上水線振蕩的周期性已不再明顯。

圖15顯示了風速6.0 m/s，7.72 m/s，9.0 m/s和13.5 m/s時上水線位置處水膜厚度變化的頻譜分析。

風速6.0 m/s時，由于上水線基本保持穩(wěn)定狀態(tài)而無明顯的振蕩頻率，導致拉索氣動力變化幅度很小且無明顯主頻，因此拉索振動很弱，不會產(chǎn)生風雨激振現(xiàn)象；風速7.72 m/s時，上水線的環(huán)向振蕩顯現(xiàn)出明顯的周期性特征，主頻為0.98 Hz，接近拉索的自振頻率，導致拉索氣動力亦發(fā)生同頻率的周期性變化，從而引發(fā)拉索的大幅度振動；當風速為9.0 m/s時，上水線的振蕩幅度明顯增大，導致氣動力變化的幅度亦明顯增大，但水線振蕩主頻為2.82 Hz，遠大于拉索自振頻率，致使拉索氣動力變化頻率亦遠大于拉索自振頻率，拉索振幅明顯減?。划旓L速增大至13.5 m/s時，上水線的振蕩頻率已增大至49.98 Hz，遠遠大于拉索自振頻率，沒有出現(xiàn)風雨激振現(xiàn)象。因此，只有當一定風速范圍內(nèi)，上水線振蕩頻率與拉索自振頻率相近時，才會導致拉索發(fā)生顯著的風雨激振現(xiàn)象。

圖14 拉索表面水膜厚度時程變化Fig.14 Temporal evolution of water film thickness

圖15 上水線頻譜分析Fig. 15 Frequency spectrum analysis of upper rivulet

4 結(jié) 論

本文考慮拉索橫風向和順風向的振動對水膜形態(tài)變化的影響，將滑移理論和兩自由度振動理論相結(jié)合，建立風雨激振時拉索和水膜耦合的理論模型；通過數(shù)值計算得到了不同風速工況下150 s的拉索振動響應、氣動力變化和水線運動情況，計算結(jié)果與人工降雨條件下拉索節(jié)段模型的風洞試驗數(shù)據(jù)吻合；在此基礎上分析風速對水線運動及拉索振動的影響，得出以下結(jié)論：

(1) 由兩自由度模型得到的拉索振動響應、氣動力及水膜形態(tài)變化規(guī)律均與風洞試驗相近，驗證了模型及數(shù)值模擬方法的準確性。

(2) 150 s的拉索振動時程表明拉索在兩個方向的振動均需要經(jīng)過幾十秒的不穩(wěn)定期后才會產(chǎn)生風雨激振現(xiàn)象，與已有的現(xiàn)場觀測結(jié)論相一致，顯示出提高數(shù)值計算效率的必要性。

(3) 隨著風速的增大，水線厚度逐漸增加，上水線環(huán)向振蕩的中心位置逐漸向背風側(cè)偏移，但拉索只在一定風速范圍內(nèi)發(fā)生大幅度振動，且其橫風向振幅接近順風向振幅的兩倍。

(4) 只有當風速處于特定范圍內(nèi)時，上水線在拉索表面環(huán)向振蕩的幅度最大，且其周期性振蕩的頻率接近拉索自振頻率，導致升力、阻力發(fā)生相近頻率的周期性變化，引發(fā)拉索在橫風向和順風向發(fā)生大幅度振動；而當風速過大或過小時，水線運動及升力、阻力變化的頻率與拉索自振頻率相去甚遠，導致拉索振動明顯減弱。因此，上水線與拉索之間的共振是導致拉索發(fā)生風雨激振的主要原因之一。