風(fēng)雨激振中斜拉索表面水線運(yùn)動(dòng)的三維數(shù)值模擬

王劍 畢繼紅 關(guān)健 喬浩玥 邵倩 周燕 管青海

摘要：基于滑移理論推導(dǎo)三維水膜運(yùn)動(dòng)方程，重點(diǎn)考慮斜拉索表面水膜形態(tài)變化對(duì)風(fēng)壓力系數(shù)和風(fēng)摩擦力系數(shù)的影響，建立模擬風(fēng)雨條件下三維節(jié)段斜拉索表面水膜形態(tài)變化的理論模型。在此基礎(chǔ)上推導(dǎo)出三維節(jié)段斜拉索氣動(dòng)力的計(jì)算公式，研究水線運(yùn)動(dòng)對(duì)拉索氣動(dòng)力及其振動(dòng)響應(yīng)的影響。結(jié)果表明：應(yīng)用三維滑移理論模型數(shù)值模擬得到的拉索表面上下水線的位置、形態(tài)與風(fēng)洞試驗(yàn)結(jié)果相近;發(fā)生風(fēng)雨激振時(shí)斜拉索的順風(fēng)向振幅明顯小于橫風(fēng)向振幅;水線沿拉索環(huán)向與軸向運(yùn)動(dòng)使得水膜形態(tài)發(fā)生低頻周期性變化，導(dǎo)致拉索氣動(dòng)升力與阻力的同頻周期性變化，從而引起拉索的大幅度振動(dòng)，產(chǎn)生風(fēng)雨激振現(xiàn)象。

關(guān)鍵詞：風(fēng)雨激振;斜拉索;三維模型;滑移理論;水線

中圖分類(lèi)號(hào)：TU312+。1;U443.38文獻(xiàn)標(biāo)志碼：A 文章編號(hào)：1004-4523（2020）03-0559-11

DOI：10.16385/i.cnki.issn.1004-4523.2020.03.015

引言

斜拉橋拉索由于具有質(zhì)量輕、剛度小和阻尼小的特點(diǎn)，在風(fēng)雨共同作用下，極易發(fā)生大幅度的低頻振動(dòng)，即風(fēng)雨激振現(xiàn)象。國(guó)內(nèi)外諸多學(xué)者通過(guò)一系列的現(xiàn)場(chǎng)觀測(cè)和風(fēng)洞試驗(yàn)研究，發(fā)現(xiàn)斜拉索表面上水線的形成與振蕩是風(fēng)雨激振現(xiàn)象的主要標(biāo)志和關(guān)鍵影響因素。為此，科研人員一方面采用人工水線風(fēng)洞試驗(yàn)測(cè)量水線對(duì)拉索氣動(dòng)性能的影響，另一方面通過(guò)人工降雨風(fēng)洞試驗(yàn)并采用各種測(cè)量方式觀測(cè)水線的運(yùn)動(dòng)特性。

由于水線尺寸較小，形態(tài)變化多種多樣，僅僅依靠風(fēng)洞試驗(yàn)研究很難全面掌握水線運(yùn)動(dòng)與斜拉索振動(dòng)間的內(nèi)在聯(lián)系，因此有必要采用數(shù)值模擬方法進(jìn)行相應(yīng)的研究。自2007年Lemaitre等首次應(yīng)用滑移理論研究水平靜止拉索表面上的水膜形態(tài)變化并模擬水線的形成以來(lái)，經(jīng)過(guò)Taylor等、許林汕等、畢繼紅和王劍等的不斷完善，滑移理論逐漸發(fā)展成為數(shù)值模擬水線運(yùn)動(dòng)的重要方法。

目前，滑移理論的研究還主要集中于二維斷面模型研究，而斜拉索的風(fēng)雨激振是復(fù)雜的三維問(wèn)題，應(yīng)考慮水的軸向流動(dòng)和氣流的軸向流動(dòng)。對(duì)此，Bi等基于滑移理論推導(dǎo)出了考慮拉索振動(dòng)和氣流作用的三維節(jié)段拉索表面的水膜運(yùn)動(dòng)方程，將拉索振動(dòng)響應(yīng)作為已知條件帶入水膜運(yùn)動(dòng)方程，重點(diǎn)研究拉索振動(dòng)對(duì)水線運(yùn)動(dòng)的影響規(guī)律。眾所周知，水線運(yùn)動(dòng)受重力、氣流作用和拉索振動(dòng)三個(gè)方面共同作用，而水膜形態(tài)變化又會(huì)對(duì)氣流作用產(chǎn)生影響。對(duì)此，本文應(yīng)用已建立的三維水膜運(yùn)動(dòng)方程，忽略拉索振動(dòng)對(duì)水膜形態(tài)的影響，而重點(diǎn)考慮斜拉索表面水膜形態(tài)變化對(duì)風(fēng)壓力系數(shù)和風(fēng)摩擦力系數(shù)的影響，建立模擬風(fēng)雨條件下三維節(jié)段斜拉索表面水膜形態(tài)變化的理論模型;在此基礎(chǔ)上推導(dǎo)出三維節(jié)段斜拉索氣動(dòng)力的計(jì)算公式，研究水線運(yùn)動(dòng)對(duì)拉索升力、阻力的影響;而后將拉索氣動(dòng)力施加到斜拉索模型上得到斜拉索的振動(dòng)響應(yīng);與已有的風(fēng)洞試驗(yàn)及數(shù)值模擬結(jié)果進(jìn)行比對(duì)驗(yàn)證，揭示風(fēng)雨激振現(xiàn)象的產(chǎn)生機(jī)理。相比于以往的基于滑移理論的研究，本文突破了二維斷面模型的限制，建立的三維模型能夠考慮水線沿斜拉索的軸向流動(dòng)，更加接近實(shí)際，可靠性更高。

1 模型

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

參考文獻(xiàn)推導(dǎo)水膜運(yùn)動(dòng)方程，半徑為R、傾角為a（0°≤a≤90°）的斜拉索節(jié)段，受水平方向氣流和重力的共同作用，如圖1所示，風(fēng)速為U，風(fēng)偏角為β（0°≤β≤90°）。

將重力分解，則作用在斜拉索斷面內(nèi)的重力分量g_N和沿斜拉索軸向的重力分量g_z分別為

采用圖2所示的柱坐標(biāo)系（e_r，e_θ，e_z），根據(jù)滑移理論，假設(shè)斜拉索表面存在一層連續(xù)的水膜，斜拉索節(jié)段任一斷面（A-A）上的水膜受力如圖3所示。水膜內(nèi)任一點(diǎn)的坐標(biāo)為（r，θ，z），R≤r≤R+h，h為水膜厚度;速度u表示為分量形式u=u_rr_r+u_θe_θ十u_ze_r，則三維Navier-Stokcs公式可寫(xiě)為

方程的邊界條件包括：

（1）位移邊界條件：水膜在底面（r-R）處相對(duì)于斜拉索表面靜止，即

式中 t為時(shí)間，u為水的動(dòng)力黏度系數(shù)，r為水在空氣中的表面張力系數(shù)（r，u與θ無(wú)關(guān)），δ和δ_g分別為水膜和空氣的應(yīng)力張量，I為單位向量，n為水膜與空氣交界處的法向向量，K為水膜表面的曲率（K=-▽·n），p為水膜內(nèi)的壓強(qiáng);p_g為水膜表面所受空氣壓力，τ_g為空氣黏滯力張量。

（3）自由邊界條件：

水膜在與空氣交界面F（r，θ，Z，t）=R+h（θ，z，t）-r=0處滿足

將邊界條件帶人式（4）并進(jìn)行無(wú)量綱處理，得到無(wú)量綱的三維水膜運(yùn)動(dòng)方程

1.2 拉索氣動(dòng)力與拉索振動(dòng)方程

斜拉索任一斷面（A-A）的受力如圖4所示，圖中F_r（θ，z）和F_θ（θ，Z）分別為水膜底面（r=R）處的法向力和切向力。水膜底面（r=R）處的應(yīng)力張量為：

整個(gè)三維剛性斜拉索節(jié)段的氣動(dòng)阻力和升力分別為

2 數(shù)值求解

2.1 風(fēng)壓力系數(shù)C_p與風(fēng)摩擦力系數(shù)C_f

在滑移理論中，氣流對(duì)水膜形態(tài)變化的影響主要體現(xiàn)在水膜運(yùn)動(dòng)方程中的風(fēng)壓力系數(shù)C_p和風(fēng)摩擦力系數(shù)C_f;而這兩個(gè)參數(shù)又會(huì)隨著水膜形態(tài)的變化而改變，如圖5所示。Bi等在應(yīng)用三維模型研究風(fēng)雨激振時(shí)忽略了水膜形態(tài)變化對(duì)風(fēng)壓力系數(shù)C_p和風(fēng)摩擦力系數(shù)C_f的影響，采用固定參數(shù)，具有一定的局限性，無(wú)法體現(xiàn)氣流與水膜運(yùn)動(dòng)之間的相互作用。由于風(fēng)雨激振中的水線位置和形狀多種多樣，很難通過(guò)風(fēng)洞試驗(yàn)來(lái)確定每一時(shí)刻的風(fēng)壓力系數(shù)和風(fēng)摩擦力系數(shù)。因此，本文仍采用二維斷面模型的相關(guān)研究方法，暫時(shí)忽略軸向氣流的作用，應(yīng)用COMSOL軟件計(jì)算每一斷面上隨時(shí)間變化的風(fēng)壓力系數(shù)和風(fēng)摩擦力系數(shù)。具體參數(shù)設(shè)置詳見(jiàn)文獻(xiàn)。

2.2 數(shù)值計(jì)算流程

圖6顯示了采用MATLAB軟件數(shù)值求解三維水膜運(yùn)動(dòng)方程（式（8））及計(jì)算拉索氣動(dòng)力（式（11））和拉索振動(dòng)響應(yīng)（式（12））的基本流程。

計(jì)算C_p和C_f時(shí)在COMSOL軟件中劃分流場(chǎng)網(wǎng)格，此過(guò)程需要用到上一時(shí)刻的水膜形態(tài)，因此在劃分網(wǎng)格時(shí)需進(jìn)行參數(shù)化處理;水膜運(yùn)動(dòng)方程（式（8））為四階非線性偏微分方程，采用預(yù)測(cè)-校正差分格式數(shù)值求解;依據(jù)式（11a）和（11b），在環(huán)向和軸向采用二重?cái)?shù)值積分方法分別計(jì)算拉索的氣動(dòng)阻力F_x和升力F_v;根據(jù)拉索振動(dòng)方程（式（12）），采用四階Runge-Kutta法求解順風(fēng)向與橫風(fēng)向振動(dòng)響應(yīng)。

2.3 基本參數(shù)

Li等采用超聲波測(cè)厚系統(tǒng)在人工降雨風(fēng)洞試驗(yàn)中測(cè)量了風(fēng)雨激振時(shí)斜拉索表面的水膜厚度分布，為此后數(shù)值研究水膜形態(tài)變化提供了實(shí)驗(yàn)基礎(chǔ)。參照該試驗(yàn)，本文研究選取的基本參數(shù)如下：斜拉索半徑R=0.05m，自振頻率f₀=0.952Hz，風(fēng)偏角β=22.5°，水膜初始厚度h₀=0.2mm，重力加速度g=9.8m/s²，水密度p=1.0×10³kg/m³，水的運(yùn)動(dòng)黏性系數(shù)v=1.0×10^-6m²/s，水在空氣中的表面張力系數(shù)γ=7.2×10^-2N/m，空氣密度p_g=1.225kg/m³，空氣的運(yùn)動(dòng)黏性系數(shù)v_g=1.51×10^-5m²/s，斜拉索的線密度p_s=8.57kg/m，阻尼比ξ₀=0.17%。

L1等在試驗(yàn)中發(fā)現(xiàn)，傾角α=30°、風(fēng)偏角β=22.5‘的斜拉索在風(fēng)速U=6.76-8.04m/s時(shí)會(huì)發(fā)生風(fēng)雨激振現(xiàn)象，在風(fēng)速U=7.72m/s時(shí)拉索振動(dòng)尤為明顯。故本文在數(shù)值計(jì)算過(guò)程中選擇傾角（α=30°、風(fēng)偏角β=22.5°、風(fēng)速U=7.72m/s作為計(jì)算參數(shù)。

根據(jù)文獻(xiàn)中三維水膜方程的求解精度研究，本文選擇0.5m長(zhǎng)的三維剛性節(jié)段斜拉索作為研究對(duì)象;沿斜拉索軸向劃分為25個(gè)斷面，每個(gè)斷面環(huán)向離散為120個(gè)點(diǎn)，時(shí)間步長(zhǎng)為10^-3s;水膜的初始厚度0.1mm。忽略雨滴的影響，參考文獻(xiàn)假設(shè)三維節(jié)段拉索表面的水量守恒，總水量（體積y）不隨時(shí)間變化，即滿足

應(yīng)用差分法求解水膜運(yùn)動(dòng)方程需要水膜在斜拉索表面連續(xù)分布，為此設(shè)定水膜的最小厚度h_min=0.01mm。

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

以往的研究表明，風(fēng)雨激振的形成需要數(shù)10s以上的時(shí)間，故本文進(jìn)行了150s的數(shù)值計(jì)算，并選取120-140s內(nèi)的計(jì)算結(jié)果進(jìn)行研究。

3.1水線的形成

圖7顯示t=0.9s時(shí)拉索表面的水膜厚度分布情況。為便于觀察，將三維斜拉索圓柱體展開(kāi)，橫軸表示斜拉索軸向方向，縱軸表示環(huán)向方向，圖中顏色表示水膜厚度。t=0.9s時(shí)，水膜在拉索下方聚集形成下水線，范圍在274.5°-280°之間，寬度約為8.8mm，厚度約為0.86mm，此時(shí)還未形成明顯的上水線。至t=1.4s時(shí)，在74°-94°范圍內(nèi)出現(xiàn)了明顯的上水線，寬度約為13.9mm，厚度約為0.33mm，如圖8所示。

3.2 三維節(jié)段斜拉索表面的水膜形態(tài)變化

圖9-14顯示了t=136.0-137.0s時(shí)間段內(nèi)斜拉索表面的水膜厚度分布。下水線位于275°-280°之間，最大厚度約為1.48mm。其位置、形態(tài)變化很小，只有厚度在斜拉索軸向方向有明顯變化，體現(xiàn)出了氣動(dòng)力和重力共同作用下水的軸向流動(dòng)。上水線則在55°-90°之間運(yùn)動(dòng)，最大厚度約為0.46mm，與風(fēng)洞試驗(yàn)的觀測(cè)結(jié)果相近。與試驗(yàn)結(jié)果不同的是，數(shù)值計(jì)算結(jié)果中上水線的環(huán)向運(yùn)動(dòng)幅值小于試驗(yàn)結(jié)果，其原因主要是計(jì)算中未考慮拉索振動(dòng)對(duì)水膜形態(tài)變化的影響;此外，數(shù)值計(jì)算結(jié)果中還出現(xiàn)了水滴自上水線位置沿迎風(fēng)側(cè)下滑并匯合到下水線的現(xiàn)象。產(chǎn)生這一現(xiàn)象的因素除計(jì)算中未考慮拉索振動(dòng)對(duì)水膜形態(tài)變化的影響外;另一因素可能是應(yīng)用滑移理論數(shù)值模擬風(fēng)雨激振時(shí)無(wú)法模擬降雨過(guò)程中雨滴下落至斜拉索表面的現(xiàn)象。

3.3 水線的環(huán)向運(yùn)動(dòng)

圖15一18分別表示t=120-140s內(nèi)Z=0.1m，Z=0.2m，Z=0.3m和Z=0.4m斷面處的水膜厚度時(shí)程變化。圖像橫軸表示時(shí)間，縱軸表示拉索環(huán)向位置，顏色表示水膜厚度?？梢钥吹竭@四幅圖像區(qū)別不大。首先，與試驗(yàn)結(jié)果相一致，下水線的形態(tài)、厚度幾乎一致，均穩(wěn)定在相同范圍內(nèi)（273°-281°）;其次，上水線寬度明顯大于下水線，而厚度明顯小于下水線，與試驗(yàn)觀測(cè)結(jié)論一致;再次，上水線的運(yùn)動(dòng)范圍在40°-90°之間，中心位置與試驗(yàn)結(jié)果一致而范圍略大;最后，這4個(gè)斷面處均有水從上水線沿迎風(fēng)側(cè)下滑匯聚至下水線，與二維斷面研究的結(jié)果一致。

為研究上水線位置處水膜形態(tài)變化的運(yùn)動(dòng)特性，分別研究Z=0.1m，Z=0.2m，Z=0.3m和Z=0.4m四個(gè)斷面θ=63°位置處的水膜厚度，如圖19-22所示。這四個(gè)斷面上水線位置處的水膜厚度變化范圍均在0.1-0.4mm范圍內(nèi);對(duì)其分別進(jìn)行頻譜分析，發(fā)現(xiàn)水膜厚度變化的主頻主要有三個(gè)：0.253，0.558和0.802Hz，均小于拉索自振頻率（0.952Hz），與試驗(yàn)結(jié)果有一定區(qū)別。對(duì)比文獻(xiàn)僅考慮拉索振動(dòng)影響的研究，可以發(fā)現(xiàn)產(chǎn)生這一差別的主要原因是本文的計(jì)算中未考慮拉索振動(dòng)的影響，說(shuō)明拉索振動(dòng)是導(dǎo)致水線周期性環(huán)向振蕩的主要因素之一。

由于在計(jì)算中出現(xiàn)了水從迎風(fēng)側(cè)自上水線下滑至下水線的現(xiàn)象，故在上述四個(gè)斷面上的迎風(fēng)側(cè)各選取同一位置（θ=9°）的水膜厚度作為研究對(duì)象，如圖23-26所示。這4個(gè)位置處的水膜厚度變化范圍在0.1-0.5mm內(nèi)，略大于上水線位置處的水膜厚度變化幅度。頻譜分析顯示迎風(fēng)側(cè)的水膜厚度變化主頻與上水線位置處一致，仍為0.253，0.558和0.802Hz。比較這4個(gè)斷面8個(gè)位置的水膜厚度變化，可以發(fā)現(xiàn)均有0.802Hz這個(gè)頻率。

3.4 水線的軸向運(yùn)動(dòng)

為研究水線的軸向運(yùn)動(dòng)，分別選取t=127-132s內(nèi)Z=0.30m和Z=0.32m兩個(gè)相鄰斷面作為研究對(duì)象。由于上水線在40°-90°之間周期性振蕩，同時(shí)有水周期性的沿迎風(fēng)側(cè)從上水線下滑匯聚到下水線，故研究上水線（θ=60°和θ=63°）及迎風(fēng)側(cè)（θ=9°和（θ=12°）的水膜厚度時(shí)程變化，如圖27所示。Z=0.30m和x=0.32m是數(shù)值計(jì)算中的兩個(gè)相鄰斷面，圖27中的對(duì)比顯示這兩個(gè)斷面的相同環(huán)向位置處的水膜厚度變化的幅值基本一致，且迎風(fēng)側(cè)的變化幅度明顯大于上水線位置處。

在同一角度位置，兩條水膜厚度時(shí)程曲線相近，但水膜最大厚度出現(xiàn)的時(shí)間明顯不同，有明顯的時(shí)間差存在，說(shuō)明水在沿著斜拉索軸向流動(dòng);比較同一斷面內(nèi)相鄰角度的時(shí)程曲線，可以發(fā)現(xiàn)水膜厚度變化亦非完全同步。因此，通過(guò)三維模型的數(shù)值計(jì)算可以看出，水線在拉索表面同時(shí)沿環(huán)向和軸向運(yùn)動(dòng)。事實(shí)上，水線沿拉索環(huán)向和軸向運(yùn)動(dòng)是在重力、氣流和拉索振動(dòng)共同作用下形成的，但在本文的研究中，忽略了拉索振動(dòng)的影響，僅研究重力和氣流作用下的水線運(yùn)動(dòng)。

3.5 拉索氣動(dòng)力

圖28和29分別為t=120-140s內(nèi)拉索氣動(dòng)升力與阻力的時(shí)程變化曲線和頻譜分析。兩者的變化趨勢(shì)一致，但升力的變化幅值較大，升力的變化范圍是—0.65-2.45N，阻力的變化范圍是4.50-6.90n.頻譜分析顯示升力與阻力的主頻均為0.253，0.558和0.802Hz，與上水線位置處的水膜厚度變化頻率一致，亦與水線下滑導(dǎo)致的迎風(fēng)側(cè)水膜厚度變化頻率一致。這一現(xiàn)象說(shuō)明水膜形態(tài)變化是導(dǎo)致拉索氣動(dòng)力變化的直接重要因素。

3.6 拉索振動(dòng)響應(yīng)

圖30所示的t=120-140s內(nèi)拉索橫風(fēng)向振動(dòng)響應(yīng)顯示橫風(fēng)向振幅約為0.114m，與風(fēng)洞試驗(yàn)觀測(cè)數(shù)據(jù)基本一致。對(duì)比圖31顯示的順風(fēng)向振動(dòng)的時(shí)程曲線，可以發(fā)現(xiàn)拉索發(fā)生風(fēng)雨激振時(shí)的順風(fēng)向振幅約為0.054m，明顯小于橫風(fēng)向振幅（0.114m），與Ni等的現(xiàn)場(chǎng)觀測(cè)結(jié)果及以往的二維斷面模型研究結(jié)論相一致，產(chǎn)生這一現(xiàn)象的主要原因是升力的變化幅值明顯大于阻力。由于氣動(dòng)升力有正有負(fù)，故拉索橫風(fēng)向振動(dòng)的平衡位置大致位于原點(diǎn)附近的0.012m處;而阻力一直為正，導(dǎo)致順風(fēng)向振動(dòng)的平衡位置遠(yuǎn)離原點(diǎn)，大致位于0.142m處。

結(jié)合水膜形態(tài)研究，可以看出風(fēng)雨激振是由于水膜形態(tài)（厚度）的低頻周期性變化導(dǎo)致拉索氣動(dòng)升力與阻力的同頻周期性變化，進(jìn)而引起拉索的大幅度振動(dòng)。

4 結(jié)論

本文基于滑移理論，忽略拉索振動(dòng)對(duì)水膜形態(tài)的影響，應(yīng)用三維水膜運(yùn)動(dòng)方程研究重力和氣流共同作用下三維斜拉索表面的水膜形態(tài)變化;并推導(dǎo)出三維節(jié)段斜拉索氣動(dòng)力的計(jì)算公式，研究水線運(yùn)動(dòng)對(duì)拉索升力、阻力及拉索振動(dòng)的影響，得到下列結(jié)論：

（1）應(yīng)用三維滑移理論模型數(shù)值模擬得到的拉索表面上下水線的位置、形態(tài)與實(shí)驗(yàn)結(jié)果相近。

（2）在三維節(jié)段模型中考慮水膜形態(tài)變化對(duì)風(fēng)壓力系數(shù)和風(fēng)摩擦力系數(shù)影響，可以較好地模擬水線沿拉索環(huán)向與軸向運(yùn)動(dòng)。上水線的環(huán)向振蕩及水線下滑導(dǎo)致上水線附近及迎風(fēng)側(cè)水膜厚度的周期性變化，二者主頻一致，均小于拉索自振頻率。

（3）阻力的變化幅度小于升力，導(dǎo)致斜拉索發(fā)生風(fēng)雨激振時(shí)的順風(fēng)向振幅明顯小于橫風(fēng)向振幅，與已有的研究結(jié)果相一致。

（4）拉索氣動(dòng)升力與阻力的變化頻率與上水線及迎風(fēng)側(cè)位置處水膜厚度變化頻率一致。水膜形態(tài)變化是導(dǎo)致拉索氣動(dòng)力變化的直接重要因素。水膜形態(tài)周期性變化導(dǎo)致拉索氣動(dòng)力的周期性變化，引發(fā)拉索橫風(fēng)向與順風(fēng)向的大幅振動(dòng)，形成風(fēng)雨激振現(xiàn)象。