關(guān) 健，李自林

(1.天津城建大學(xué) 土木工程學(xué)院，天津市西青區(qū)津靜路26號 300384；2.河北水利電力學(xué)院 土木工程學(xué)院，河北省滄州市重慶路1號 061001)

近年來我國斜拉橋工程發(fā)展迅猛，隨著橋梁跨度不斷增加，斜拉索的剛度和阻尼越來越小,因此斜拉索極易發(fā)生振動的問題也日益突出，其中風(fēng)雨激振振幅最大，危害最為嚴(yán)重。所謂風(fēng)雨激振是指在一定天氣條件下，斜拉索受到風(fēng)雨的共同作用，之后發(fā)生的低頻率、大幅度的振動現(xiàn)象。自從1984年，Hikami等[1]首次觀測到該現(xiàn)象以來，世界多地區(qū)都出現(xiàn)了類似報道。隨后很多學(xué)者，如Yamaguchi[2]、Seidel和Dinkler[3]、Xu和Wang[4]、顧明[5]、李永樂[6]等紛紛提出多種理論對風(fēng)雨激振產(chǎn)生機(jī)理以及相關(guān)參數(shù)進(jìn)行研究。研究結(jié)果表明，上水線的形成和振蕩對風(fēng)雨激振的產(chǎn)生和發(fā)展有重要影響。由于實際水線在斜拉索表面的形態(tài)多樣，變化快，尺度小，難以準(zhǔn)確測量，單純依靠現(xiàn)場觀測和風(fēng)洞試驗的方法不足以全方位掌握其變化規(guī)律。隨著計算機(jī)技術(shù)的快速發(fā)展，計算流體力學(xué)(CFD)數(shù)值模擬技術(shù)逐漸成為研究復(fù)雜風(fēng)雨激振問題的有力工具之一。李壽英和顧明等[7-8]通過流體計算軟件CFX5.5建立了斜拉索的數(shù)值模型。敬海泉等[9]根據(jù)邊界層狀態(tài)提出一個新的數(shù)值模型來描述拉索表面的風(fēng)場，系統(tǒng)分析了多個因素對風(fēng)雨激振特性的影響。陳文禮等[10-11]應(yīng)用自編譯程序，提出氣泡破裂的理論并對風(fēng)雨激振現(xiàn)象進(jìn)行模擬。Lemaitre等[12]首次基于滑移理論，初步建立了二維水平靜止斜拉索表面的水膜運(yùn)動方程，并初步模擬了水線的形成。由于該理論能準(zhǔn)確反映斜拉索表面水膜變化形態(tài)，多名學(xué)者[13-18]在此基礎(chǔ)上先后考慮了重力、表面張力、斜拉索傾斜角度、斜拉索振動、隨時間變化的氣動力等因素的影響，現(xiàn)已逐步發(fā)展為研究風(fēng)雨激振問題的重要方法。

文中應(yīng)用基于滑移理論的斜拉索風(fēng)雨激振理論模型，考慮斜拉索振動與斜拉索表面水膜形態(tài)變化兩者之間的關(guān)系，并應(yīng)用COMSOL與MATLAB相結(jié)合的方法計算隨時間變化的氣動力系數(shù)，大幅提高了計算效率。通過42組數(shù)值模型，研究了不同風(fēng)速以及風(fēng)偏角條件下，斜拉索振動響應(yīng)、斜拉索表面水膜形態(tài)變化以及隨時間變化的空氣動力三者之間的關(guān)系。

1 模型

1.1 水膜運(yùn)動方程

如圖1(a)所示：風(fēng)速大小為U0，風(fēng)偏角為β，其中β的變化范圍為0≤β≤90°；斜拉索半徑為Rc、水平傾角為α，其中α的變化范圍為0<α≤90°；受到重力和空氣動力的共同作用，風(fēng)攻角為0°。在斜拉索上取A-A斷面，以拉索圓心為原點(diǎn)建立極坐標(biāo)系(er,eθ)。假設(shè)拉索表面有一層連續(xù)的水膜，在任意θ角處，水膜的厚度為h(θ)，該截面上拉索表面水膜受力如圖1(b)所示。

(a)斜拉索空間位置

(b)斜拉索表面水膜受力圖1 斜拉索模型圖

忽略軸向重力分量，在該斷面內(nèi)的重力分量：

gN=gcosα

(1)

若只考慮二維斷面內(nèi)的氣流作用，則斜拉索斷面內(nèi)的風(fēng)速：

(2)

δ=arctan(sinα·tanβ)

(3)

水膜的運(yùn)動方程由二維Navier-Stokes公式推導(dǎo)。建立極坐標(biāo)系(er,eθ)，如圖1(b)所示，水膜內(nèi)點(diǎn)的坐標(biāo)為(r,θ)，其中R≤r≤R+h，h為水膜厚度。將速度u表示為分量形式u=urer+uθeθ，則N-S方程可寫為：

(4a)

(4b)

(4c)

根據(jù)Lemaitre等[12]的計算結(jié)果可知，水膜運(yùn)動的雷諾數(shù)Reh=hu/v≈1(v為水的運(yùn)動粘性系數(shù))，式(4a)和式(4b)可寫為：

ρ[gNsin(θ-δ)-sinθ]

(5a)

(5b)

方程的邊界條件包括以下方面。

(1)水膜底面采用無滑移邊界條件，即

ur|r=Rc=uθ|r=Rc=0

(6)

(2)水膜在與空氣交界面F(r,θ,t)=Rc+

h(θ,t)-r=0處所受的應(yīng)力與表面張力相平衡，即

(σg-σ)n=κγn

(7a)

σ=-pI+[μu+(u)T]

(7b)

σg=-pg(θ)I+τg

(7c)

式中：γ為水在空氣中的表面張力系數(shù)，σ和σg分別為水膜和空氣的應(yīng)力張量，I為單位向量，pg為水膜表面所受空氣壓力，τg為空氣粘滯力張量；n為水膜與空氣交界處的法向向量，κ為水膜表面的曲率。

除了滿足以上2種邊界條件之外，拉索表面水膜與空氣的交界面還滿足自由界面F(r,θ,t)=Rc+h(θ,t)-r=0的連續(xù)條件，即水膜在與空氣交界面處滿足

(8)

將邊界條件和式(8)帶入方程(5a)、(5b)和(4c)中可得水膜運(yùn)動方程：

(9)

1.2 斜拉索振動方程

圖2為斜拉索A-A斷面受力圖。設(shè)拉索質(zhì)量為M，在二維模型中為單位長度上的質(zhì)量，阻尼為c，剛度為K，阻尼比為ξ0，自振圓頻率為ω0，根據(jù)單自由度振動理論在橫風(fēng)向(y軸方向)建立拉索運(yùn)動方程：

(10)

式中Fy為橫風(fēng)向的拉索升力，按下式計算：

(11)

圖2 斜拉索斷面受力

水膜底面的應(yīng)力張量為：

σ0={-pI+μ[u+(u)T]}|r=R

(12)

式中I為單位向量，則

Fr(θ)=er·σ0·er=-p+2μ?rur|r=R

(13a)

(13b)

將式(13a)和式(13b)分別代入式(11)，即可得到升力Fy。

1.3 隨時間變化的空氣動力的計算

滑移理論假設(shè)拉索外表面有一層連續(xù)的水膜，因此水膜外表面與空氣接觸，內(nèi)表面與拉索接觸。水膜作為空氣和拉索的媒介，在受到兩者共同影響的同時也反作用于兩者。已有研究結(jié)果表明[15-17]，考慮隨時間變化的空氣動力系數(shù)能更加準(zhǔn)確地反映風(fēng)雨激振現(xiàn)象。文中沿用畢繼紅等[18]的研究方法，采用COMSOL計算隨時間變化的風(fēng)壓力系數(shù)Cp和風(fēng)摩擦力系數(shù)Cf。在每一個時間步長內(nèi)，依據(jù)該時刻水膜形態(tài)建立“類圓柱”外表面，并認(rèn)為該時間步長內(nèi)水膜形態(tài)采用定常流，采用SA湍流模型進(jìn)行計算，計算區(qū)域設(shè)置、網(wǎng)格劃分和邊界條件設(shè)置都與文獻(xiàn)[18]一致。

2 數(shù)值求解

2.1 水膜運(yùn)動方程的求解

方程(9)為水膜運(yùn)動方程，該方程為四階偏微分方程，采用有限差分法進(jìn)行求解。根據(jù)Robertson等人的研究[13]，將水膜在圓周空間上離散為128點(diǎn)可以保證計算精度，文中也照例設(shè)置，即相鄰兩個離散點(diǎn)與圓柱中心的連線夾角為Δθ=0.711π。時間步長取為0.001 s。水膜運(yùn)動方程離散后的差分格式和差分方程可參考文獻(xiàn)[18]。

2.2 數(shù)值求解流程

數(shù)值求解流程如圖3所示。計算伊始，根據(jù)拉索和水膜的初始狀態(tài)輸入初始值，參數(shù)無量綱化后進(jìn)入迭代過程。在迭代過程中，首先將計算結(jié)果中的水膜厚度數(shù)據(jù)讀取出來；然后以此為依據(jù)在COMSOL中建立幾何模型，在模型中劃分網(wǎng)格，并根據(jù)上一步的拉索速度，設(shè)定邊界條件從而建立數(shù)值模型；然后計算該時刻的水膜外表面各點(diǎn)處的氣動力系數(shù)；再依據(jù)上一步求出的拉索加速度，求解水膜運(yùn)動方程，即可獲得下一時刻各點(diǎn)的水膜厚度；將水膜厚度變化結(jié)果代入方程(13)求解Fr(θ)和Fθ(θ)，進(jìn)而求解方程(11)獲得拉索升力Fy；然后將升力代入方程(10)，采用四階Runge-Kutta法求解，可獲得拉索新的位移、速度以及加速度。接下來判斷是否進(jìn)入最終時刻，直到計算全部結(jié)束。

2.3 基本參數(shù)選取

根據(jù)陳文禮等[10]的試驗，斜拉索表面各個位置處的水膜初始厚度h0=0.2 mm，應(yīng)用有限差分法求解水膜運(yùn)動方程需要水膜在斜拉索表面連續(xù)分布，因此在水膜形態(tài)演變過程中，最大水膜厚度設(shè)為1 mm，最小水膜厚度設(shè)為0.02 mm。斜拉索斷面上的水膜可以近似認(rèn)為滿足質(zhì)量守恒。文中選取1×10-3s作為數(shù)值計算的時間步長。其余基本參數(shù)見表1，計算工況見表2。

圖3 計算流程圖

表1基本參數(shù)

Tab.1Parametersincalculations

重力加速度g/(m·s-2)斜拉索半徑R/m斜拉索傾角α/(°)斜拉索自振頻率f0/Hz空氣的運(yùn)動粘性系數(shù)vg空氣密度ρg/(kg·m-3)水的運(yùn)動粘性系數(shù)v/(m2·s-1)水密度ρ/(kg·m-3)水的表面張力系數(shù)γ/(N·m-1)9.80.05300.9521.51×10-51.2251×10-610007.2×10-2

表2 數(shù)值計算工況

3 數(shù)值模擬結(jié)果

3.1 斜拉索振動響應(yīng)

陳文禮等[10]在哈爾濱工業(yè)大學(xué)風(fēng)洞與浪槽聯(lián)合實驗室中進(jìn)行了人工降雨風(fēng)洞試驗。將風(fēng)偏角為22.5°的計算結(jié)果提取出來并與試驗結(jié)果相對比，結(jié)果如圖4所示。其中深色線為試驗結(jié)果，淺色線為計算結(jié)果?？梢园l(fā)現(xiàn)，風(fēng)速過大或者過小時，斜拉索振幅都不大，在風(fēng)速為7.4～7.72 m/s范圍內(nèi)時，斜拉索振幅大于0.08 m。這體現(xiàn)出斜拉索風(fēng)雨激振“限速”的特點(diǎn)?？傮w上斜拉索振幅與試驗結(jié)果非常接近。驗證了本文中提出的斜拉索風(fēng)雨激振理論模型和求解方法的準(zhǔn)確性和可靠性。

圖4 各個風(fēng)速下的拉索振幅與試驗結(jié)果對比Fig.4 The cable amplitude comparison of calculated results and test results at every wind speed

圖5給出了6個典型風(fēng)速下的斜拉索振動響應(yīng)時程曲線，可以發(fā)現(xiàn)，當(dāng)風(fēng)速較小時，拉索振動迅速趨于穩(wěn)定，振幅很小。當(dāng)風(fēng)速進(jìn)入激振區(qū)范圍內(nèi)時，斜拉索振幅初期也是比較小，經(jīng)過大約30 s的能量積累之后，振幅基本穩(wěn)定在最大值附近。表3是對斜拉索振動時程曲線進(jìn)行頻譜分析的結(jié)果，該結(jié)果表明，在風(fēng)速較小以及在風(fēng)雨激振范圍內(nèi)時，斜拉索振動主頻接近斜拉索自振頻率0.952 Hz，其頻譜分析中的幅值逐漸增加；在風(fēng)速超過激振區(qū)范圍后，主頻開始發(fā)生偏移，且幅值減小，說明曲線頻譜特性不明顯。

(a)UN=6.00 m/s (b)UN=6.76 m/s (c)UN=7.40 m/s

(d)UN=7.72 m/s (e)UN=8.04 m/s (f)UN=10.00 m/s圖5 六個風(fēng)速下的拉索振幅時程曲線Fig.5 The time history of cable vibration at six wind speeds

表3斜拉索振動主頻(UN=7.72 m/s)

Tab.3Thedominantfrequencyofcablevibration(UN=7.72 m/s)

風(fēng)速/(m·s-1)6.06.767.47.57.67.728.048.38.59.010主頻/Hz0.9890.9890.9890.9890.9890.9890.9890.9891.033無無幅值/(°·Hz-1)0.0120.0150.0740.0740.0740.0750.0530.0280.003無無

圖6給出的是風(fēng)速在6.00 m/s到9.00 m/s的范圍內(nèi)，斜拉索振幅隨風(fēng)偏角的變化趨勢，可以發(fā)現(xiàn)風(fēng)偏角為22.5°，風(fēng)速為7.72 m/s時的振幅最大。當(dāng)風(fēng)速為6 m/s和9 m/s時，風(fēng)偏角在10°～50°變化時，斜拉索振幅都很小。這說明相同風(fēng)偏角的情況下，隨著風(fēng)速改變都是先增大后減小，風(fēng)偏角的改變對風(fēng)雨激振區(qū)間影響不大。圖7顯示的是風(fēng)速為7.72 m/s時，斜拉索振幅隨風(fēng)偏角改變的變化規(guī)律。其先增大后減小的特征明顯，在風(fēng)偏角為22.5°時，斜拉索振幅最大，接近0.1 m。

圖6 不同風(fēng)偏角下的斜拉索振幅Fig.6 The amplitudes of the cable vibrations at different wind yaw angles

圖7 風(fēng)速為7.72m/s時，斜拉索振幅隨風(fēng)偏角的變化規(guī)律Fig.7 The variation of the amplitude of cable vibration with the wind deflection at the wind speed of 7.72m/s

3.2 上水線形態(tài)變化

圖8顯示的是不同風(fēng)速以及風(fēng)偏角條件下的上水線厚度。由圖8可以發(fā)現(xiàn)，在風(fēng)偏角相同的條件下，上水線厚度隨風(fēng)速增加而增大；在風(fēng)速不變的情況下，上水線厚度隨著風(fēng)偏見增加而減小。這個規(guī)律與斜拉索振動響應(yīng)規(guī)律有較大不同。為分析其原因，對上水線厚度變化時程進(jìn)行頻譜分析，并將風(fēng)偏角在10°～50°、風(fēng)速在6.00～9.00 m/s范圍內(nèi)的主頻總結(jié)于表4中。由表4可以發(fā)現(xiàn)，在風(fēng)速小于等于7.72 m/s范圍內(nèi)時，上水線厚度變化主頻均為0.989 Hz，接近斜拉索自振頻率，而當(dāng)風(fēng)速進(jìn)一步增加以后，上水線主頻逐漸偏離斜拉索自振頻率。圖9顯示的是風(fēng)偏角為22.5°時，3個典型風(fēng)速下的斜拉索表面水膜形態(tài)演變規(guī)律。由圖9可以發(fā)現(xiàn)，在風(fēng)速為6.00 m/s時，斜拉索表面并未形成明顯上水線。當(dāng)風(fēng)速增加到7.72 m/s，斜拉索振動穩(wěn)定后，上水線在斜拉索表面周期性地形成和消失，其演變頻率與斜拉索自振周期相接近。當(dāng)風(fēng)速增大到10.00 m/s后，上水線依然會形成，且最大水膜厚度顯然比風(fēng)速為7.72 m/s時更大，但同時可以發(fā)現(xiàn)，上水線形態(tài)雜亂，其周期性并不明顯。

圖8 不同風(fēng)速以及風(fēng)偏角條件下的上水線厚度Fig.8 Thickness of water film at different wind speeds and wind deflection angles

表4 上水線附近水膜厚度變化主頻

(a)UN=6.00 m/s (b)UN=7.72 m/s (c)UN=10.00 m/s圖9 3個典型水膜形態(tài)變化Fig.9 Three typical water film evolutions

3.3 斜拉索表面升力變化

圖10是各個工況下的斜拉索表面升力幅值。由此可知：斜拉索表面升力隨風(fēng)速增加呈增大趨勢；在風(fēng)速不變的前提下，斜拉索升力隨風(fēng)偏角增大而減小。這是由于風(fēng)偏角增加，斷面內(nèi)風(fēng)速會減小，導(dǎo)致氣動力減小，進(jìn)而引起斜拉索升力下降。升力主頻規(guī)律在表5中給出，由表5可知：風(fēng)速較小以及風(fēng)雨激振發(fā)生時，各個工況中升力變化主頻均接近斜拉索自振頻率；風(fēng)速較大時，升力主頻逐漸偏離斜拉索自振頻率，這也是斜拉索振動主頻偏移自振頻率的主要原因，因為此時斜拉索受迫振動效果更明顯。而斜拉索表面升力偏移斜拉索自振頻率的原因是上水線運(yùn)動主頻逐漸消失。

圖10 各個工況下的斜拉索表面升力幅值Fig.10 Lift on the cable surface at every case

表5 不同風(fēng)偏角時的升力主頻

3.4 斜拉索和水線的雙向耦合作用分析

由式(9)可以看出，斜拉索表面水膜運(yùn)動同時受到斜拉索振動和空氣動力的影響。當(dāng)風(fēng)速較小時，斜拉索振動對水膜運(yùn)動的影響是主要的，而水膜運(yùn)動對于斜拉索振動反饋作用相對較小，因此斜拉索表面受到的升力(由于水膜運(yùn)動產(chǎn)生)較小，此時斜拉索振動類似自由振動。這也是斜拉索表面水膜厚度變化主頻接近斜拉索自振頻率的原因。隨著風(fēng)速的增大，水膜運(yùn)動逐漸劇烈，水膜運(yùn)動對斜拉索振動的作用也逐漸增加，但由于此時水膜運(yùn)動主頻依然接近斜拉索自振頻率，因此此時水膜與斜拉索發(fā)生了同頻諧振，導(dǎo)致斜拉索振幅大幅增加，風(fēng)雨激振由此產(chǎn)生。當(dāng)風(fēng)速增大到激振區(qū)范圍外時，水膜運(yùn)動對于斜拉索振動的影響成為主要方面，由于此時斜拉索表面水膜在空氣動力作用下變化劇烈，水膜形態(tài)雜亂，其周期性已不再明顯。雖然此時上水線厚度比之前更大，但由于水膜與斜拉索之間的同頻諧振現(xiàn)象已經(jīng)消失，斜拉索振幅并未隨著升力幅值的進(jìn)一步增加而增大，反而逐漸減小，斜拉索逐漸退出風(fēng)雨激振狀態(tài)。

4 結(jié)論

(1)不同風(fēng)偏角的條件下，隨風(fēng)速增加，斜拉索振幅都是先增加后減小的趨勢。風(fēng)偏角的改變對風(fēng)雨激振的激振區(qū)間范圍影響不大。

(2)在激振區(qū)間內(nèi)，斜拉索振幅隨風(fēng)偏角的增加呈現(xiàn)先增加后減小的趨勢，在風(fēng)偏角為22.5°左右時，斜拉索振幅最大。

(3)相同風(fēng)速下，隨著風(fēng)偏角的增加，上水線厚度以及升力幅值呈減小趨勢。風(fēng)偏角不變時，上水線厚度以及升力幅值均隨風(fēng)速增加而增大，風(fēng)速超過激振范圍外后，兩者主頻均逐漸偏離斜拉索自振頻率。

(4)斜拉索與上水線之間的同頻諧振是斜拉索風(fēng)雨激振產(chǎn)生的主要原因之一。