斜拉索風雨激振數值模擬及其機理探究

畢繼紅， 武 驥， 關 健， 王 劍

(1. 天津大學 建筑工程學院， 天津 300072； 2. 濱海土木工程結構與安全教育部重點實驗室， 天津 300072;3. 天津城建大學 土木工程學院， 天津 300384； 4. 天津軟土特性與工程環境重點實驗室， 天津 300384)

斜拉索風雨激振數值模擬及其機理探究

畢繼紅1,2， 武 驥1， 關 健1， 王 劍3,4

(1. 天津大學 建筑工程學院， 天津 300072； 2. 濱海土木工程結構與安全教育部重點實驗室， 天津 300072;3. 天津城建大學 土木工程學院， 天津 300384； 4. 天津軟土特性與工程環境重點實驗室， 天津 300384)

將滑移理論與CFD技術相結合，考慮了風壓力系數和風摩擦因數隨位置和時間的變化，并引入了Spalart-Allmaras(S-A)湍流模型，運用有限單元法，實現了對斜拉索的風雨激振現象的數值模擬；考察了不同風速下的水線形成及形態，并與已有實驗成果進行比對，驗證了方法的合理性和準確性；研究了水膜形態、氣動力變化、斜拉索振動響應三者關系，得到了水線和拉索之間的振動耦合規律，探討了風雨激振機理。結果表明，水線、氣流和斜拉索之間的共振是引發風雨激振的主要因素之一。

風雨激振； 水線； 氣動力； 斜拉索振動

風雨激振是指在風雨共同作用下，斜拉索發生的低頻、大幅度振動，危害極大。對此，研究人員運用現場實測、風洞試驗、理論分析方法對此進行了大量研究。然而，各國內外學者對其產生機理并未達成統一。顧明等[1]進行了斜拉索人工雨線測壓試驗, 建立了三維拉索,三維水線風雨激振模型,分析認為平均氣動力系數的突降是致振機理。劉慶寬等[2]認為，卡門渦的脫落受到軸向流的抑制,低頻渦的產生得到激勵,受其影響，拉索產生渦激發散振動。Cheng 等[3]認為，在臨界雷諾數狀態下,尾流渦脫規律性發生變化，產生氣動升力，致使拉索振動。Bi等[4]認為，影響風雨激振現象產生的重要因素是上水線的產生和其振蕩。

近來，由于現場觀測的隨機性和不確定性，風洞實驗再現風雨激振的困難性，使得運用CFD數值模擬技術分析風雨激振機理得到發展。然而，現有的CFD數值模擬方法有諸多不足，需要進一步改進。其一，個別數值模擬方法需調用多個計算軟件，且各軟件之間兼容性不佳，模型參數和計算結果需人工大范圍處理，以致工作量極大，極易出現錯誤；其二，風雨激振是氣、液、固三種相態相互耦合的復雜問題，因方法限制和計算精度要求，以致計算時間步長極短，計算耗時冗長。其三為簡化計算，個別方法采用固定風壓力和風摩擦因數，忽略了水線形態和位置變化對其的動態影響，不能完全反映水線形成的物理場景，這與實際不符。針對上述弊病，本文做如下改進：① 運用COMSOL軟件建模，MATLAB計算，在風雨激振的數值模擬中將二者無縫鏈接，減少計算和人為錯誤；② 使用差分法解水膜方程，保證精度的前提下，增大計算時間步長，采用S-A湍流模型，提升計算效率；③ 克服采用固定氣動力系數的弊病，采用隨時間和位置變化的風壓力和風摩擦因數即瞬態氣動力系數，貼近風雨激振物理場景，提高準確性。在此基礎上,建立了拉索風雨激振的理論模型并與已有實驗進行了對比驗證。分析了風雨激振條件下的水膜形態、氣動力變化、斜拉索振動響應三者關系，并探究了風雨激振的產生機理。

1 模型建立

1.1 幾何模型

圖1 斜拉索空間位置

取圖1中的斜拉索A-A斷面為研究對象，如圖1,2所示。忽略重力沿斜拉索軸向的分量，則作用在斜拉索斷面內的重力分量gN：

(1)

圖2 斜拉索表面水膜受力

忽略軸向流的影響，只考慮垂直于斜拉索的氣流作用，則斜拉索斷面內的風速UN：

(2)

(3)

1.2 基于滑移理論的水膜運動方程

根據滑移理論，認為斜拉索表面存在水膜[5-6]，采用極坐標系(er,e0)下的采用極坐標系下Navier-Stokes公式(4)～(6)描述水膜在拉索表面的運動，初始水膜厚度為h0，厚度為h(θ,t)，水膜內任意一點坐標為(r，θ),其速度可表示為u=urer+uθe0。

(4)

(5)

(6)

經代入邊界條件，結合滑移理論進行簡化，并經過無量綱化處理后，在重力、表面張力、風壓力和風摩擦力下的單向耦合水膜運動方程為(7)，并通過差分法進行求解;

(7)

(8)

(9)

式中：u為水膜內的流場速度;ρg和ρ分別為空氣和水的密度;p為水膜內的壓強;ν為水的動力黏度系數;γ為水在空氣中的表面張力系數;σg為空氣的應力張量;Pg為水膜表面所受的空氣壓力;n和t分別為水膜與空氣交界處的法向向量和切向向量。其具體的無量綱化見文獻[7]。

1.3 應用COMSOL技術獲得風壓力系數Cp和風摩擦因數Cf

此前的計算，Cp和Cf多為靜態固定不變，不能準確反映拉索表面氣動力的變化，以致不能準確模擬水膜水線的形成和分析風雨激振的真正機理。為彌補此缺陷，運用COMSOL軟件計算隨時間和空間位置動態變化的Cp和Cf，并利用MATLAB編程自動調用瞬態Cp和Cf進行計算。

選取的計算區域如圖3所示，設拉索直徑為D=2R，計算區域取邊長為27.5D×20D的矩形，右側為入口，左側為出口，將拉索中心置于坐標原點，拉索中心距離入口10D，距離出口17.5D，距離上下邊界10D。計算流場網格劃分如圖4。采用RANS(Reynolds Average Navier-Stokes)湍流數值模擬方法，進行穩態計算。采用Spalart-Allmaras湍流模型。邊界條件設置如下：拉索為無滑移固壁邊界；上部和下部均為滑移邊界；右側為速度入口邊界，設置風場風速為UN；左側為出口邊界，邊界條件為壓力P=0，并且抑制回流。在保證計算結果準確性的前提下，較其他方法增大時間步長，取為dt=0.001 s,提升計算效率。

圖3 流場計算區域

圖4 流場網格示意圖

1.4 計算參數

參照LI等[8-10]的實驗工況選取模型基本參數、模型驗證參數分別如表1、表2；根據斜拉橋設計規范[11]建議的斜拉索設計宜用傾角區間和常見自然環境風向，選取拉索傾角和風偏角。

表1 模型基本參數

表2 模型驗證參數

2 模型驗證

本文將滑移理論與CFD技術相結合，建立有限元模型模擬水線形成。為驗證此方法的可靠性和準確性，本文考察不同風速下上水線的形成，考察水線初始形成形態，并與實驗數據進行比對。如圖5，水線形成時間隨風速增大逐漸變短，在風速小于8.04 m/s時，水線形成時間隨風速增大而急劇變短；當風速大于8.04 m/s時，水線形成時間隨風速增大而平緩變化，上述規律與Bi[12]所述規律相同。

如圖6所示，穩定的上水線形成位置處在50.63°～ 81.56°內，與實測測數據52.48°～ 82.51°符合很好。如圖7所示，上水線平均高度0.52 mm，與實驗觀測到的0.509 mm[10]非常接近。如圖8所示，上水線平均寬度7.34 mm，亦非常接近實驗觀測到的水線寬度7.96 mm。說明，其一，本文模擬方法準確性高、穩定性強，由此方法推導的結論有很強的可靠性。其二，風速是上水線振蕩和水膜形態變化的顯著影響因素。

圖5 不同風速下的水線初始形成時間

圖6 不同風速下的水線初始形成位置

圖7 不同風速下的水線高度

3 水線模擬結果分析

3.1 斜拉索表面水膜形態變化

學界普遍認為風雨激振主要受上水線的形成和振蕩影響，為揭示風雨激振機理，本文選取了激振風速U0=7.72 m/s的工況，重點對上水線作了研究分析。如圖9所示，t=0.27 s左右開始形成上水線，初始形成位置約為67.5°，其振蕩區間大致位于55°～85°，高度約為0.49 mm，寬度約為9.64 mm。上水線寬度、高度發展較大時，由于重力，上水線產生向下滑落現象，引起水膜、水線形態變化。

圖8 不同風速下的水線寬度

(a) 2D水膜形態

(b) 3D水膜形態

圖10顯示了在0°位置處，風速在U0=7.72 m/s時水膜的高度時程變化和其頻譜分析。如圖10(a)所示：水膜形態出現了兩階段變化：t=0 s～20 s為劇烈變化階段，由于水線的形成、滑落劇烈，使得水線震蕩、水膜厚度起伏變化幅值較大；t=20 s以后為平穩變化階段，氣流產生的風壓力和風摩擦力的變化趨于平穩，水線振蕩和水膜厚度變化逐漸顯示了明顯的規律性和穩定性。如圖10所示，可以發現水膜厚度其振蕩頻率f≈1 Hz，與試驗中與斜拉索的自振頻率f0≈0.952 Hz相近。

(a) 時程曲線

(b) 頻譜分析

3.2 斜拉索表面氣動力

t=60 s處于平穩變化階段，并以此時刻為研究對象，如圖11所示，風壓力系數Cp、風摩擦因數Cf、表面風速、風壓力等高線在50°～85°范圍內發生了顯著變化，這與圖9顯示的上水線的形成和振蕩位置相吻合。這表明，氣動力的突變區域位于上水線的振蕩區域，分析得知，上水線的振蕩對Cp、Cf有重大影響，受到影響的氣動力又反作用于水膜形態，同時影響拉索振動，各因素相互作用，密切聯系。進一步說明考慮Cp、Cf隨時間和位置的變化對準確模擬風雨激振的必要性和重要性。

圖12(a), 13(a)顯示了U0=7.72 m/s時采用固定Cp、Cf時的氣動力變化：氣動升力波動區間為-0.9 N～-0.4 N，氣動阻力波動區間為5.0 N～5.2 N,與采用瞬態Cp、Cf時的氣動力相比,振幅明顯偏小。這表明是否采用瞬態Cp、Cf進行計算，對氣動力的計算結果影響很大，進而很大程度上影響風雨激振數值模擬的準確性。

(a) 風壓力系數Cp

(b) 風摩擦因數Cf

(c) 表面風速隨位置變化圖示(m/s)

(d) 氣動力等高線隨位置變化圖示(Pa)

圖12(b),13(b)顯示了U0=7.72 m/s時的采用瞬態Cp、Cf斜拉索氣動力變化。如圖所示，氣動力時程顯示了與水膜形態變化一致的階段性：t=0 s～20 s為劇烈變化階段，氣動升力從-0.4 N～0.8 N快速變大為-0.6 N～1.0 N，氣動阻力從4.7 N～6.5 N快速降為5 N～6.3 N。t=20 s以后為平穩變化階段，氣動力逐步穩定，并且升力幅值處在-0.6 N～1.2 N之間，阻力氣動力幅值處在5.2 N～5.8 N之間，變化穩定，周期性明顯。分析認為，最初，水膜形態在初始風速的作用下劇烈變化，導致氣動力幅值劇烈變化；此后，水膜形態變化和氣動力變化相互耦合，使得較初始狀態變化更為緩和；最后，水線的振蕩和水膜形態逐漸趨于穩定，導致氣動力變化逐漸趨于穩定，穩定變化的氣動力又反過來作用于水線，保證水膜形態的穩定變化，各因素相互作用、相互協調，達到平衡，協同保證了此現象的穩定性和規律性。

如圖12(c)和圖13(c)所示：升力和阻力的振蕩頻率分別為fy≈1 Hz和fx≈1 Hz，與LI等[8-10]的實驗結果相近。氣動力振蕩頻率與0°處的水膜厚度振蕩頻率相同再一次證明，氣動力的周期性振蕩與水膜形態周期性變化密切相關。

(a) 升力時程曲線(固定Cp、Cf)

(b) 升力時程曲線(瞬態Cp、Cf)

(c) 升力頻譜分析(瞬態Cp、Cf)

(a) 阻力時程曲線(固定Cp、Cf)

(b) 阻力時程曲線(瞬態Cp、Cf)

(c) 阻力頻譜分析(瞬態Cp、Cf)

3.3 斜拉索振動響應

圖14(a)顯示了采用固定Cp、Cf的拉索振幅時程：其平衡位置在-0.25 m處，波動區間為-0.10 m～0.06 m，與實驗結果[8]差別較大(平衡位置大約在0 m處，波動區間在-0.12 m～0.13 m左右)。如圖14(b)所示，采用瞬態Cp、Cf的振動響應顯示了與水膜形態、氣動力變化一致的階段性：t=0 s～20 s為劇烈變化階段，在劇烈變化的氣動力的作用下，拉索振幅迅速變大到-0.14 m～0.15 m；t=20 s以后為平穩變化階段，周期性明顯，平衡位置在0 m處，振幅區間為-0.13 m～0.15 m，與實驗結果[8]符合良好，這表明，采用瞬態Cp、Cf更貼近風雨激振物理環境，使得計算結果更為準確，數值模擬更為合理。

如圖14(c)所示：拉索振動響應振蕩頻率為fy≈1 Hz。綜合結果得知，水線的振蕩頻率、氣動力頻率、拉索的振動響應頻率，斜拉索的自振頻率四者相近。另外根據陳文禮[9]的實驗得知，U0=7.72 m/s位于風洞試驗中發生風雨激振現象的顯著風速范圍內，結合水膜形態變化和氣動力變化說明了風雨激振現象的氣、液、固耦合機制：在激振風速范圍內，包裹在拉索表面的水膜主要受氣流和重力的作用形成上水線，并呈現周期振蕩和滑落；水膜水線的周期性變化又導致了氣動力的周期性變化；周期性變化的氣動力和水線振蕩作用在斜拉索上，又由于其變化頻率與拉索自振頻率相近，從而引發了拉索的受迫振動，發生共振；反過來，斜拉索的振動又影響水膜形態和氣動力的變化，三者相互作用、相互影響。

(a) 拉索振幅時程曲線(固定Cp、Cf)

(b) 拉索振幅時程曲線(瞬態Cp、Cf)

(c) 拉索振幅頻譜分析(瞬態Cp、Cf)

4 結 論

本文將滑移理論與CFD技術相結合, 并基于S-A湍流模型，考慮隨時間和位置變化的風壓力和風摩擦因數，對斜拉索的風雨激振現象進行了數值模擬。通過考察水膜形態、氣動力和拉索振動響應之間的關系，探究了風雨激振的產生機理，得到以下結論。

(1) 風速是上水線振蕩和水膜形態變化的顯著影響因素:風速在6 m/s～10 m/s范圍內，隨風速增大，水線形成時間變短，形成位置向背風側移動；水線高度、寬度發生不同程度的增大和波動。

(2) 風雨激振狀態下，上水線及水膜形態存在周期性的形成、滑落和振蕩，并呈現了階段性的發展過程：劇烈變化階段、平穩變化階段。

(3) 氣動力的突變區域與上水線的振蕩區域基本重合，上水線的振蕩對風壓力系數Cp和風摩擦因數Cf有重大影響，Cp和Cf受到影響的氣動力又反作用于水膜形態，同時影響拉索振動。考慮Cp和Cf隨時間和位置的變化和與水膜形態拉索振動的相互作用，對風雨激振的準確計算極其重要。

(4) 水線振蕩、氣動力變化、拉索振動相互作用，由于頻率相近引發共振，是產生風雨激振現象重要因素之一。

[1] 顧明,李壽英，杜曉慶.斜拉橋拉索風雨激振理論模型和機理研究[J].空氣動力學學報，2007,25(2)169-174.

GU Ming,LI Shouying, DU Xiaoqing. Testing study on wind pressure distributions of stayed cables with a fixed artificial rivulet[J]Acta Aerodynamica Sinica,2007,25(2)169-174.

[2] 劉慶寬, 張峰, 喬富貴. 軸向流對斜拉索氣動穩定性影響的試驗研究[J]. 石家莊鐵道學院學報(自然科學版), 2008, 21(4): 16-19.

LIU Qingkuan, ZHANG Feng, QIAO Fugui. Effect of axial flow on rain-wind induced vibration of stay-cables[J].Journal of Shijiazhuang Railway Institute (NaturalScience), 2008, 21(4): 16-19.

[3] CHENG S, IRWIN P A, JAKOBSENL B, et al. Divergent motion of cables exposed to skewed wind[C] //Proc. of the 5th In.t Symposium on Cable Dynamics. Italy: Santa Margherita Ligure, 2003: 271-278

[4] BI J H，WANG J，SHAO Q，et al. 2D numerical analysis on evolution of water film and cable vibration response subject to wind and rain[J]. Journal of Wind Engineering and Industrial Aerodynamics，2013,121：49-59.

[5] LEMAITRE C，HéMON P，DE LANGRE E. Thin water film around a cable subject to wind[J]. Journal of Wind Engineering and Industrial Aerodynamics，2007， 95(9/10/11)：1259-1271.

[6] LEMAITRE C，DE LANGRE E，HéMON P. Rainwater rivulets running on a stay cable subject to wind[J].European Journal of Mechanics-B/Fluids,2010，29(4)：251-258.

[7] 王劍.斜拉索表面水膜形態及斜拉索氣動力變化規律[J].天津大學學報(自然科學與工程技術版), 2014,47(4)：1-12.

WANG Jian. Variation of water film morphology and aerodynamic force of stay cable[J]. Journal of Tianjin University(Science and Technology),2014，47(4)：1-12.

[8] LI F C，CHEN W L，LI H, et al. An ultrasonic transmission thickness measurement system for study of water rivulets characteristics of stay cables suffering from wind-rain-induced vibration[J]. Sensors and Actuators A：Physical，2010，159(1)：12-23.

[9] 陳文禮. 斜拉索風雨激振的試驗研究與數值模擬[D]. 哈爾濱: 哈爾濱工業大學土木工程學院, 2009.

[10] LI H，CHEN W L，XU F. A numerical and experimental hybrid approach for the investigation of aerodynamic forces on stay cables suffering from rain-wind induced vibration[J]. Journal of Fluids and Structures，2010，26(7/8)：1195-1215.

[11] 中華人民共和國交通部.公路斜拉橋設計規范(試行):JTJ 027—96[S].北京：人民交通出版社，1996.

[12] 畢繼紅，王劍，逯鵬，等.不同風速下拉索表面水線的形成[J].天津大學學報(自然科學與工程技術版), 2014, 47(7)：577-582.

BI Jihong, WANG Jian, LU Peng, et al. Formation of rivulets on cable surface under different wind speeds[J].Journal of Tianjin University(Science and Technology, 2014, 47(7)：577-582.

Numerical simulation of rain-wind induced vibration of stay cables and its mechanism study

BI Jihong1, 2, WU Ji1, GUAN Jian1, WANG Jian3,4

(1. School of Civil Engineering，Tianjin University，Tianjin 300072，China；2. Key Laboratory of Coast Civil Structure Safety，Ministry of Education，Tianjin University，Tianjin 300072，China:3. School of Civil Engineering，Tianjin Chengjian University，Tianjin 300384，China；4. Key laboratory of Soft Soil Characteristic and Engineering Environment of Tianjin，Tianjin 300384，China)

Using the lubrication theory combined with the computational fluid dynamics(CFD) method, considering wind pressure coefficient and wind friction one changes with variation of position and time, introducing Spalart-Allmaras(S-A) turbulence flow model and adopting the finite element method, the numerical simulation for rain-wind induced vibration (RWIV) phenomena of stay cables was realized. To verify the reasonableness and correctness of this method，the formation and morphology of rivulets at different wind speeds were inspected，and they were compared with the existing test results. Through investigating the relations among water film morphology, lift variation and vibration responses of stay cables, vibration coupling laws between rivulets and stay cables were obtained, the mechanism of RWIV were explored. Similar to test results, the simulation ones showed that the resonance among upper rivulet, lift and stay cables is one of the main factors to cause RWIV.

rain-wind induced vibration；rivulet；lift；vibration of cable

國家自然科學基金(51408399)

2015-12-30 修改稿收到日期：2016-03-25

畢繼紅 女，博士，教授,1965年生

U443.38

A

10.13465/j.cnki.jvs.2017.11.017